# The Definition Method and Optimization of Atomic Strain Tensors for Nuclear Power Engineering Materials

Authors: Xiangguo Zeng; Ying Sheng; Huayan Chen; Tixin Han
Journal: Science and Technology of Nuclear Installations
 (2016)

---

## Abstract

A common measure of deformation between atomic scale simulations and the continuum framework is provided and the strain tensors for multiscale simulations are defined in this paper. In order to compute the deformation gradient of any atomm, the weight function is proposed to eliminate the different contributions within the neighbor atoms which have different distances to atom m, and the weighted least squares error optimization model is established to seek the optimal coefficients of the weight function and the optimal local deformation gradient of each atom. The optimization model involves more than 9 parameters. To guarantee the reliability of subsequent parameters identification result and lighten the calculation workload of parameters identification, an overall analysis method of parameter sensitivity and an advanced genetic algorithm are also developed.


---

## Body

## 1. Introduction

Titanium alloys have been largely used as nuclear power engineering materials [1], and it is important and significant to analyze the atomic-level strain distribution of these materials. The strain tensors are commonly defined by the local deformation of the continuum. Unlike displacement, strain is not a physical quantity that can be measured directly, and it is calculated from a definition that relies on the gradient of the continuous displacement field. At the microscale, it is difficult to define the local deformation according to the position of each atom which is obtained from the adjacent discrete time interval, so there is no universally accepted definition of strain tensors of atomic scale so far.Many engineering problems involving physical phenomena need to calculate strain tensors at atomic scale. Wang et al. [2] pointed that it was important to analyze the atomic-level strain distribution and get the atomic stress-strain curve while studying the mechanical behavior of Zr-based metallic glass under indentation. Hirth et al. [3] considered that the computation of the deformation gradient and strain tensors made the approach useful for evaluation of continuum models, development of microstructure and mechanical property relationships, and identification of dislocations and disclinations, as well as for quantification of plastic spin and strain gradients.In recent years, many researchers are challenging to provide a common measure of deformation between atomic scale simulations and the continuum framework and define the strain tensors for multiscale simulations. Zimmerman et al. [4] defined the slip vector according to the positions of atoms and successfully identified the lattice distortion and the formation of dislocation structures, but these measures could not be utilized in the continuum framework. Mott et al. [5] presented a definition of the local atomic strain increments in three dimensions and an algorithm for computing them. First, an arbitrary arrangement of atoms was tessellated into Delaunay tetrahedra, and then the deformation gradient increment tensor for interstitial space was obtained from the displacement increments of the corner atoms of Delaunay tetrahedra. However, it was complicated to establish the tetrahedral elements of atoms. Gullett et al. [6] proposed an atomic strain tensor that is based on the definition of a discrete equivalent to the continuum deformation gradient that accounts for the relative motion of an atom and its neighbors in a nonlocal fashion. This method was computationally efficient, because the deformation gradient arose from an optimization procedure that did not rely on a geometric decomposition, and the strain tensors were computed directly from the deformation gradient and were appropriate for general finite, multiaxial deformation states. When the deformation gradient at an atom was formed, a weight function should be built to eliminate the different contributions within the neighbor atoms which had different distances to the atom. However, Gullett et al. [6] did not study the establishment and optimization method of the weight function which played an important role in the formulation of the discrete deformation but just used the invariant weight function of the artificial assumption to calculate the discrete deformation gradient at the atom.By summarizing the shortcomings of the existing methods mentioned above, the work done by this paper can be categorized into three parts: first, the strain tensors for multiscale simulations are defined, and the weighted least squares error optimization model involving more than 9 parameters is established to seek the optimal coefficients of the weight function and the optimal local deformation gradient of each atom; next, to guarantee the reliability of subsequent parameters identification result and also to lighten the calculation workload of parameters identification, an overall analysis method of parameter sensitivity, based on Latin Hypercube Sampling method and Spearman rank correlation method, is proposed; furthermore, on the fundamentals of the result of parameter sensitivity and basic genetic algorithm (GA), an advanced genetic algorithm based on the advanced niche genetic algorithm, global peak value determination strategy, and local accurate searching techniques is developed. Finally, taking alpha titanium as an example, the strain tensors of atoms are computed by means of the method proposed in this paper, and the method is proved to be correct and feasible by comparing with the results got by other methods from the existing reference.


## 2. Modeling Approach

### 2.1. Deformation Gradient and Strain Tensors

In order to describe the positions of atoms at the initial timet
0 and at the current time t, we assume a fixed Cartesian coordinate system as shown in Figure 1.Figure 1
General motion in the neighborhood of a discrete atomic particle.In Figure1, X and x are, respectively, the coordinate vectors of atoms in the reference configuration Ω
0 at the initial time t
0 and in the current configuration Ω
1 at the current time t, and χ maps the atoms from X to x. Consider(1)
X
=
X
i
E
i
,
x
=
x
i
e
i
,
x
=
χ
X
.Assuming sufficient continuity, the local deformation atX is characterized as the gradient of the motion and can be defined as(2)
F
≡
∂
χ
∂
X
=
∂
x
∂
X
.The deformation of an infinitesimal segmentd
x at one point in the reference configuration can be expanded in a Taylor series. If the higher-order terms are omitted, d
x can be written as(3)
d
x
=
χ
X
+
d
X
-
χ
X
=
χ
X
+
∇
χ
X
·
d
X
+
o
d
X
-
χ
X
≈
∇
χ
X
·
d
X
≈
F
·
d
X
.Then the deformation gradientF can be got by (3) which is not only applicable to atomic scale but also available for the continuum framework. Consider an atomic system shown in Figure 1. The deformation in the neighborhood of atom m is characterized by the changes in the relative position of its neighbors. Atom m is located at the position X
m in the reference configuration Ω
0 and position x
m in the current configuration Ω
1. Then the relative position of neighboring atom n and the deformation gradient F
m at atom m are given as(4)
Δ
X
m
n
=
X
n
-
X
m
Δ
x
m
n
=
x
n
-
x
m
Δ
x
m
n
=
F
m
·
Δ
X
m
n
.The deformation gradient of atomm is related to its neighboring atoms, and the deformation gradient of each neighboring atom with respect to atom m is different, so the deformation gradient F
m of atom m cannot generally be got by a single atom n. Therefore, we should seek an optimal local deformation gradient F
^
m, which can make the weighted squares error W
m shown in the following equation minimized:(5)
W
m
=
∑
n
=
1
N
Δ
x
m
n
-
F
^
m
Δ
X
m
n
T
Δ
x
m
n
-
F
^
m
Δ
X
m
n
ω
n
,where N is the number of neighboring atoms and ω
n is a weight factor.The optimal deformation gradientF
^
m of the atom m is determined by the deformation gradients of all the neighboring atoms with respect to atom m, but the contribution of each atom n to F
^
m is different and related to the distance between atom n and atom m. Therefore, the weight function used to eliminate the different contributions within the neighbor atoms which had different distances to the atom. In the cutoff radius r
c
u
t that specifies the domain around atom m, ω
n of the nearest neighbor atom is equal to 1.0 and gradually reduced to 0.0 with the increase of the distance between atom n and atom m. The weight function ω
n plays an important role on calculating the atomic strain. For each material, the weight function curve must be uniquely determined, but it is difficult to get the exact analytic expression of the weight function corresponding to the weight function curve, so many forms of ω
n can be assumed and the coefficients of the function should be optimized to fit the weight function curve. The S-curve model of biological population growth is proved to be one of the appropriate forms to fit the weight function curve [7], so the appropriate weight function ω
n can be assumed as(6)
ω
n
=
2
+
k
1
1
+
k
1
-
1
1
+
k
1
e
-
k
2
r
r
=
R
m
n
-
r
g
1
r
c
u
t
,where k
1 and k
2 are undetermined coefficients, R
m
n is the distance between atom n and atom m, and r
g
1 is the distance between the closest neighbor atom and atom m.The purpose of parameters identification is to find a set of solutions of discrete design variablesX to satisfy the objective function W
m
(
X
) as shown in (7)
o
b
j
:
m
i
n
W
m
X
=
∑
n
=
1
N
Δ
x
m
n
-
F
^
m
Δ
X
m
n
T
Δ
x
m
n
-
F
^
m
Δ
X
m
n
ω
n
X
=
F
^
11
m
,
F
^
12
m
,
F
^
13
m
,
F
^
21
m
,
F
^
22
m
,
F
^
23
m
,
F
^
31
m
,
F
^
32
m
,
F
^
33
m
,
k
1
,
k
2
T
s
.
t
.
:
X
~
i
m
i
n
≤
X
i
≤
X
~
i
m
a
x
.When the objective functionW
m
(
X
) reaches the minimum, the parameters F
^
i
j
m (i
,
j
=
1,2
,
3) are the optimal components of the optimal deformation gradient matrix F
^
m.WhenF
^
m is determined by the optimization model as shown in (7), the Lagrangian Green strain tensor E is defined with respect to reference coordinates in terms of this quantity as(8)
E
=
1
2
F
T
F
-
I
.The strain tensorE is a 3 × 3 matrix, and the effective strain of a point ε
e
f
f can be described by the second invariants of E, shown as(9)
ε
e
f
f
=
3
J
2
=
3
E
11
E
12
E
21
E
22
+
E
22
E
23
E
32
E
33
+
E
11
E
13
E
31
E
33
.


### 2.2. Parameter Sensitivity Analysis

The optimization model (see (7)) involves 11 parameters. To guarantee the reliability of subsequent parameters identification result and also to lighten the calculation workload of parameters identification, an overall analysis method of parameter sensitivity, based on Latin Hypercube Sampling method and Spearman rank correlation method, is proposed.Sensitivity analysis means picking out those sensitive parameters, which could significantly affect the reliability of result, among wide range of uncertain factors. A traditional analysis method is called single factor analysis method. And the basic procedures are as follows [8].(1)
Build systematic model, meaning function relation between system feature and factors,P
=
W
(
x
1
,
x
2
,
…
,
x
m
), and get the benchmark sets of parameters.
(2)
When analyzing the effect of one particular parameterx
i on feature P, keep the rest of parameters constant in basic value, and only allow parameter x
i to vary within the proper range. If minor change of x
i leads to drastic changes of P, that means P is very sensitive to x
i, and the sensitivity value of x
i is large; if x
i varies within a rather wide range, while the change in P is not obvious, that means P is not sensitive to x
i, and the sensitivity value of x
i is small.The optimization model (see (7)) obtained in this paper is nonlinear, and parameters are related to each other. All the parameters collectively affect the fitting accuracy of the objective function W
m
(
X
) as shown in (7). Hence, overall analysis of parameter sensitivity is required, which means allowing every parameter simultaneously to change within assigned ranges, to observe the effect of each parameter on objective function.To thoroughly analyze all the changing parameters, first, we need to pick up samples from the whole parameters space. Latin Hypercube Sampling (LHS) is a uniformed space sampling method, proposed earliest by McKay et al. in 1979 [9], which possesses advantages of fewer sampling times and effectively avoiding repeat sampling, and applies to complicated multidimensional parameters space sampling. The basic procedures of LHS are as follows.(1)
Set a parameter set consisting ofm parameters, and set the value range of each parameter as n nonoverlap zones with equal probability.
(2)
Randomly pick a sample of every parameterx
i within its range, meaning each parameter x
i consists of n sample values.
(3)
Randomly arrange then sample values of each x
i, to form m random arrangements.
(4)
Pick up one sample value for each parameter from them arrangements, which could form a sample parameter set X
k and pick n times in sequence, and then one will get n sample parameter sets {
X
1
,
X
2
,
…
,
X
n
}.The realization process schematic drawing of LHS method is shown in Figure2.Figure 2
The realization process schematic drawing of LHS method.Since the optimization model (see (7)) obtained in this paper is nonlinear, the output result of random sampling sample sets is uncertain. Therefore, it is necessary to adopt nonparameter statistical method to carry out the relativity analysis among random parameters. Spearman rank correlation method is a nonparameter statistical analysis method, proposed by Spearman in 1904 [10]. This method applies to relativity analysis among multiple parameters, to analyze the effect of input parameters on the output result. This is a very practical method [11]. The basic procedures of Spearman rank correlation method are as follows.(1)
Definition of rank: by arranging the sample sets{
A
1
,
A
2
,
…
,
A
n
} from large to small, an ordered sequence A is formed. Supposing A
i is at the order of R
i in the sequence A, then R
i is the rank of A
i in {
A
1
,
A
2
,
…
,
A
n
}.
(2)
By substitution ofn sets of sample parameters {
X
1
,
X
2
,
…
,
X
n
} into (7), the solution sets of objective function {
W
1
,
W
2
,
…
,
W
n
} can be obtained.
(3)
Supposing the rank of the parameterx
i
k in {
x
i
1
,
x
i
2
,
…
,
x
i
n
} is α
k and the rank of the parameter W
k in {
W
1
,
W
2
,
…
,
W
n
} is β
k, then the sensitivity of parameter x
i equals the absolute value of Spearman rank correlation coefficient, as shown in(10)
r
i
=
n
∑
k
=
1
n
a
k
b
k
-
∑
k
=
1
n
a
k
∑
k
=
1
n
b
k
n
∑
k
=
1
n
a
k
2
-
∑
k
=
1
n
a
k
2
n
∑
k
=
1
n
b
k
2
-
∑
k
=
1
n
b
k
2
.Arrange the sensitivity degreer
i of all the parameters from large to small, the larger sensitivity degree means greater impact on output result, which requires close attention in subsequent parameter identification analysis.The result of parameter sensitivity analysis obtained through this method forms the basic theory for the purpose of determining the density of optimization model parameter discrete zone. This method not only reduces calculation workload of subsequent parameter identification but also guarantees the reliability of the result of parameter identification.


### 2.3. Parameter Identification Method

Based on the result of parameter sensitivity analysis, in order to achieve identification of the undetermined parameters in the optimization model (see (7)), one possible method is genetic algorithm (GA). GA is a self-adapting overall optimization probability searching algorithm based on biological genetic and evolution process in the nature. Compared to traditional optimization algorithms, GA has better overall searching ability. The basic feature of GA is using overall evolution, which means to find the best solution through constant evolution of species.The controlling condition of basic GA includes individual coding method, fitness evaluation function, initial group, group size, the selection operator, the crossover operator, the mutation operator, and termination condition of genetic algorithm.However, there are many drawbacks in basic GA; for example, the accuracy, reliability, and calculating time cannot be satisfied simultaneously, and it is likely to encounter disadvantages such as earliness and poor local searching ability, so basic GA requires further update.This paper made a lot of effective update to basic GA, adopting the advanced niche genetic algorithm, suspicious peak value determination strategy, and local accurate searching techniques, so that the overall searching ability of genetic algorithm is improved, which can be shown as follows.(
1) Improvement to Traditional Niche Genetic Algorithm. First, build optimal individual preservation strategy. Each niche evolves individually and saves the current optimal individual after evolution is done for each generation. When evolution of this generation is done before the evolution of each generation starts, if the optimal individual is within the group, then leave it as it is; otherwise, use the copy of current optimal individual to replace the worst individual in current generation. This method guarantees that the optimal individual will not be eliminated during the process of evolution and will speed up the converge process.Otherwise, build data exchange mechanism among niches. Perform optimal individual exchange when each evolution is completed, meaning use the optimal individual of the first niche to replace the worst individual in the second niche. Under the premise of ensuring the diversity of the population, this method can improve the proportion of the good individuals and speed up the convergence.(
2) Building Suspicious Peak Value Determination Strategy. When independent evolution of each niche is done, it will converge to a peak value. When the amount of niche is large, then all the peak values must contain all the global optimization and local optimization points. Before the final global optimization point is determined, all the obtained peak values are considered as “suspicious peak value.” The optimization point of A amount of all the suspicious peak value point must be the global optimization point, which is referred to as objective function F
∗, and the corresponding variables are {
x
1
∗
,
x
2
∗
,
…
,
x
m
∗
}. Set the objective function of a suspicious peak value point B as F, and the corresponding variables are {
x
1
,
x
2
,
…
,
x
m
}. The condition that point B is another global optimization point different from A is(11a)
F
-
F
∗
F
∗
<
α
(11b)
d
i
=
∑
k
=
1
n
x
i
-
x
i
∗
2
>
β
i
i
=
1,2
,
…
,
m
,where α
(
0
<
α
<
1
) is a constant, meaning the relative searching range around the global optimization value; β
i is a constant, which is used to judge whether two points are neighboring; and d
i is the Hamming distance between two points. If d
i
≤
β
i, it means that points A and B are neighboring, and they are essentially the same global optimization point with different accuracy. Equation (11a) is used to determine that B is also a global optimization point, and (11b) is used to determine that B is another global optimization point different from A.(
3) Local Accurate Searching for All the Different Global Minimum. Supposing the variable sets of one global optimization point are {
x
1
∗
,
x
2
∗
,
…
,
x
m
∗
} and the range of each variable is x
i
l
≤
x
i
≤
x
i
u, the new searching range of each variable will be changed to(12)
x
i
∈
x
i
l
,
x
i
u
∩
x
i
∗
-
γ
i
,
x
i
∗
+
γ
i
i
=
1,2
,
…
,
m
,where γ
i is the searching range around the variable x
i.Finally, the genetic evolution operation will be performed to the variables within the new searching range. Until the termination condition is satisfied, the new solution got must have higher accuracy.The program flow diagram of advanced GA is shown in Figure3.Figure 3
The program flow diagram of advanced GA.


## 2.1. Deformation Gradient and Strain Tensors

In order to describe the positions of atoms at the initial timet
0 and at the current time t, we assume a fixed Cartesian coordinate system as shown in Figure 1.Figure 1
General motion in the neighborhood of a discrete atomic particle.In Figure1, X and x are, respectively, the coordinate vectors of atoms in the reference configuration Ω
0 at the initial time t
0 and in the current configuration Ω
1 at the current time t, and χ maps the atoms from X to x. Consider(1)
X
=
X
i
E
i
,
x
=
x
i
e
i
,
x
=
χ
X
.Assuming sufficient continuity, the local deformation atX is characterized as the gradient of the motion and can be defined as(2)
F
≡
∂
χ
∂
X
=
∂
x
∂
X
.The deformation of an infinitesimal segmentd
x at one point in the reference configuration can be expanded in a Taylor series. If the higher-order terms are omitted, d
x can be written as(3)
d
x
=
χ
X
+
d
X
-
χ
X
=
χ
X
+
∇
χ
X
·
d
X
+
o
d
X
-
χ
X
≈
∇
χ
X
·
d
X
≈
F
·
d
X
.Then the deformation gradientF can be got by (3) which is not only applicable to atomic scale but also available for the continuum framework. Consider an atomic system shown in Figure 1. The deformation in the neighborhood of atom m is characterized by the changes in the relative position of its neighbors. Atom m is located at the position X
m in the reference configuration Ω
0 and position x
m in the current configuration Ω
1. Then the relative position of neighboring atom n and the deformation gradient F
m at atom m are given as(4)
Δ
X
m
n
=
X
n
-
X
m
Δ
x
m
n
=
x
n
-
x
m
Δ
x
m
n
=
F
m
·
Δ
X
m
n
.The deformation gradient of atomm is related to its neighboring atoms, and the deformation gradient of each neighboring atom with respect to atom m is different, so the deformation gradient F
m of atom m cannot generally be got by a single atom n. Therefore, we should seek an optimal local deformation gradient F
^
m, which can make the weighted squares error W
m shown in the following equation minimized:(5)
W
m
=
∑
n
=
1
N
Δ
x
m
n
-
F
^
m
Δ
X
m
n
T
Δ
x
m
n
-
F
^
m
Δ
X
m
n
ω
n
,where N is the number of neighboring atoms and ω
n is a weight factor.The optimal deformation gradientF
^
m of the atom m is determined by the deformation gradients of all the neighboring atoms with respect to atom m, but the contribution of each atom n to F
^
m is different and related to the distance between atom n and atom m. Therefore, the weight function used to eliminate the different contributions within the neighbor atoms which had different distances to the atom. In the cutoff radius r
c
u
t that specifies the domain around atom m, ω
n of the nearest neighbor atom is equal to 1.0 and gradually reduced to 0.0 with the increase of the distance between atom n and atom m. The weight function ω
n plays an important role on calculating the atomic strain. For each material, the weight function curve must be uniquely determined, but it is difficult to get the exact analytic expression of the weight function corresponding to the weight function curve, so many forms of ω
n can be assumed and the coefficients of the function should be optimized to fit the weight function curve. The S-curve model of biological population growth is proved to be one of the appropriate forms to fit the weight function curve [7], so the appropriate weight function ω
n can be assumed as(6)
ω
n
=
2
+
k
1
1
+
k
1
-
1
1
+
k
1
e
-
k
2
r
r
=
R
m
n
-
r
g
1
r
c
u
t
,where k
1 and k
2 are undetermined coefficients, R
m
n is the distance between atom n and atom m, and r
g
1 is the distance between the closest neighbor atom and atom m.The purpose of parameters identification is to find a set of solutions of discrete design variablesX to satisfy the objective function W
m
(
X
) as shown in (7)
o
b
j
:
m
i
n
W
m
X
=
∑
n
=
1
N
Δ
x
m
n
-
F
^
m
Δ
X
m
n
T
Δ
x
m
n
-
F
^
m
Δ
X
m
n
ω
n
X
=
F
^
11
m
,
F
^
12
m
,
F
^
13
m
,
F
^
21
m
,
F
^
22
m
,
F
^
23
m
,
F
^
31
m
,
F
^
32
m
,
F
^
33
m
,
k
1
,
k
2
T
s
.
t
.
:
X
~
i
m
i
n
≤
X
i
≤
X
~
i
m
a
x
.When the objective functionW
m
(
X
) reaches the minimum, the parameters F
^
i
j
m (i
,
j
=
1,2
,
3) are the optimal components of the optimal deformation gradient matrix F
^
m.WhenF
^
m is determined by the optimization model as shown in (7), the Lagrangian Green strain tensor E is defined with respect to reference coordinates in terms of this quantity as(8)
E
=
1
2
F
T
F
-
I
.The strain tensorE is a 3 × 3 matrix, and the effective strain of a point ε
e
f
f can be described by the second invariants of E, shown as(9)
ε
e
f
f
=
3
J
2
=
3
E
11
E
12
E
21
E
22
+
E
22
E
23
E
32
E
33
+
E
11
E
13
E
31
E
33
.


## 2.2. Parameter Sensitivity Analysis

The optimization model (see (7)) involves 11 parameters. To guarantee the reliability of subsequent parameters identification result and also to lighten the calculation workload of parameters identification, an overall analysis method of parameter sensitivity, based on Latin Hypercube Sampling method and Spearman rank correlation method, is proposed.Sensitivity analysis means picking out those sensitive parameters, which could significantly affect the reliability of result, among wide range of uncertain factors. A traditional analysis method is called single factor analysis method. And the basic procedures are as follows [8].(1)
Build systematic model, meaning function relation between system feature and factors,P
=
W
(
x
1
,
x
2
,
…
,
x
m
), and get the benchmark sets of parameters.
(2)
When analyzing the effect of one particular parameterx
i on feature P, keep the rest of parameters constant in basic value, and only allow parameter x
i to vary within the proper range. If minor change of x
i leads to drastic changes of P, that means P is very sensitive to x
i, and the sensitivity value of x
i is large; if x
i varies within a rather wide range, while the change in P is not obvious, that means P is not sensitive to x
i, and the sensitivity value of x
i is small.The optimization model (see (7)) obtained in this paper is nonlinear, and parameters are related to each other. All the parameters collectively affect the fitting accuracy of the objective function W
m
(
X
) as shown in (7). Hence, overall analysis of parameter sensitivity is required, which means allowing every parameter simultaneously to change within assigned ranges, to observe the effect of each parameter on objective function.To thoroughly analyze all the changing parameters, first, we need to pick up samples from the whole parameters space. Latin Hypercube Sampling (LHS) is a uniformed space sampling method, proposed earliest by McKay et al. in 1979 [9], which possesses advantages of fewer sampling times and effectively avoiding repeat sampling, and applies to complicated multidimensional parameters space sampling. The basic procedures of LHS are as follows.(1)
Set a parameter set consisting ofm parameters, and set the value range of each parameter as n nonoverlap zones with equal probability.
(2)
Randomly pick a sample of every parameterx
i within its range, meaning each parameter x
i consists of n sample values.
(3)
Randomly arrange then sample values of each x
i, to form m random arrangements.
(4)
Pick up one sample value for each parameter from them arrangements, which could form a sample parameter set X
k and pick n times in sequence, and then one will get n sample parameter sets {
X
1
,
X
2
,
…
,
X
n
}.The realization process schematic drawing of LHS method is shown in Figure2.Figure 2
The realization process schematic drawing of LHS method.Since the optimization model (see (7)) obtained in this paper is nonlinear, the output result of random sampling sample sets is uncertain. Therefore, it is necessary to adopt nonparameter statistical method to carry out the relativity analysis among random parameters. Spearman rank correlation method is a nonparameter statistical analysis method, proposed by Spearman in 1904 [10]. This method applies to relativity analysis among multiple parameters, to analyze the effect of input parameters on the output result. This is a very practical method [11]. The basic procedures of Spearman rank correlation method are as follows.(1)
Definition of rank: by arranging the sample sets{
A
1
,
A
2
,
…
,
A
n
} from large to small, an ordered sequence A is formed. Supposing A
i is at the order of R
i in the sequence A, then R
i is the rank of A
i in {
A
1
,
A
2
,
…
,
A
n
}.
(2)
By substitution ofn sets of sample parameters {
X
1
,
X
2
,
…
,
X
n
} into (7), the solution sets of objective function {
W
1
,
W
2
,
…
,
W
n
} can be obtained.
(3)
Supposing the rank of the parameterx
i
k in {
x
i
1
,
x
i
2
,
…
,
x
i
n
} is α
k and the rank of the parameter W
k in {
W
1
,
W
2
,
…
,
W
n
} is β
k, then the sensitivity of parameter x
i equals the absolute value of Spearman rank correlation coefficient, as shown in(10)
r
i
=
n
∑
k
=
1
n
a
k
b
k
-
∑
k
=
1
n
a
k
∑
k
=
1
n
b
k
n
∑
k
=
1
n
a
k
2
-
∑
k
=
1
n
a
k
2
n
∑
k
=
1
n
b
k
2
-
∑
k
=
1
n
b
k
2
.Arrange the sensitivity degreer
i of all the parameters from large to small, the larger sensitivity degree means greater impact on output result, which requires close attention in subsequent parameter identification analysis.The result of parameter sensitivity analysis obtained through this method forms the basic theory for the purpose of determining the density of optimization model parameter discrete zone. This method not only reduces calculation workload of subsequent parameter identification but also guarantees the reliability of the result of parameter identification.


## 2.3. Parameter Identification Method

Based on the result of parameter sensitivity analysis, in order to achieve identification of the undetermined parameters in the optimization model (see (7)), one possible method is genetic algorithm (GA). GA is a self-adapting overall optimization probability searching algorithm based on biological genetic and evolution process in the nature. Compared to traditional optimization algorithms, GA has better overall searching ability. The basic feature of GA is using overall evolution, which means to find the best solution through constant evolution of species.The controlling condition of basic GA includes individual coding method, fitness evaluation function, initial group, group size, the selection operator, the crossover operator, the mutation operator, and termination condition of genetic algorithm.However, there are many drawbacks in basic GA; for example, the accuracy, reliability, and calculating time cannot be satisfied simultaneously, and it is likely to encounter disadvantages such as earliness and poor local searching ability, so basic GA requires further update.This paper made a lot of effective update to basic GA, adopting the advanced niche genetic algorithm, suspicious peak value determination strategy, and local accurate searching techniques, so that the overall searching ability of genetic algorithm is improved, which can be shown as follows.(
1) Improvement to Traditional Niche Genetic Algorithm. First, build optimal individual preservation strategy. Each niche evolves individually and saves the current optimal individual after evolution is done for each generation. When evolution of this generation is done before the evolution of each generation starts, if the optimal individual is within the group, then leave it as it is; otherwise, use the copy of current optimal individual to replace the worst individual in current generation. This method guarantees that the optimal individual will not be eliminated during the process of evolution and will speed up the converge process.Otherwise, build data exchange mechanism among niches. Perform optimal individual exchange when each evolution is completed, meaning use the optimal individual of the first niche to replace the worst individual in the second niche. Under the premise of ensuring the diversity of the population, this method can improve the proportion of the good individuals and speed up the convergence.(
2) Building Suspicious Peak Value Determination Strategy. When independent evolution of each niche is done, it will converge to a peak value. When the amount of niche is large, then all the peak values must contain all the global optimization and local optimization points. Before the final global optimization point is determined, all the obtained peak values are considered as “suspicious peak value.” The optimization point of A amount of all the suspicious peak value point must be the global optimization point, which is referred to as objective function F
∗, and the corresponding variables are {
x
1
∗
,
x
2
∗
,
…
,
x
m
∗
}. Set the objective function of a suspicious peak value point B as F, and the corresponding variables are {
x
1
,
x
2
,
…
,
x
m
}. The condition that point B is another global optimization point different from A is(11a)
F
-
F
∗
F
∗
<
α
(11b)
d
i
=
∑
k
=
1
n
x
i
-
x
i
∗
2
>
β
i
i
=
1,2
,
…
,
m
,where α
(
0
<
α
<
1
) is a constant, meaning the relative searching range around the global optimization value; β
i is a constant, which is used to judge whether two points are neighboring; and d
i is the Hamming distance between two points. If d
i
≤
β
i, it means that points A and B are neighboring, and they are essentially the same global optimization point with different accuracy. Equation (11a) is used to determine that B is also a global optimization point, and (11b) is used to determine that B is another global optimization point different from A.(
3) Local Accurate Searching for All the Different Global Minimum. Supposing the variable sets of one global optimization point are {
x
1
∗
,
x
2
∗
,
…
,
x
m
∗
} and the range of each variable is x
i
l
≤
x
i
≤
x
i
u, the new searching range of each variable will be changed to(12)
x
i
∈
x
i
l
,
x
i
u
∩
x
i
∗
-
γ
i
,
x
i
∗
+
γ
i
i
=
1,2
,
…
,
m
,where γ
i is the searching range around the variable x
i.Finally, the genetic evolution operation will be performed to the variables within the new searching range. Until the termination condition is satisfied, the new solution got must have higher accuracy.The program flow diagram of advanced GA is shown in Figure3.Figure 3
The program flow diagram of advanced GA.


## 3. Results and Discussion

### 3.1. Example Strain Calculation

In order to compute the strain tensors of atoms and provide evidence to confirm the method proposed in this paper to be correct and feasible by comparing with the result got by another method from the existing reference [6], the atomic mechanics model of alpha titanium material is built first. The length, width, and thickness of the model are separately 250 Å, 100 Å, and 3 Å, and the tensile displacement load on both upper and lower surface of the model is 0.025 Å/ps, as shown in Figure 4.Figure 4
The atomic mechanics model.The coordinate vectors of atoms in the reference configuration and in the current configuration are separately got at the time oft
=
0 and t
=
10 ps. Taking the center atom numbered 175 with a cutoff radius of 2 nm as an example, the number of atoms within the cutoff radius around the center atom is 146.The weight function in [6] was assumed as(13)
ω
r
=
1
-
6
r
2
+
6
r
3
r
≤
0.5
2
-
6
r
+
6
r
2
-
2
r
3
0.5
<
r
<
1.0
0
r
≥
1.0
.Then the optimal deformation gradient matrixF
^ could be got by (14)
F
^
=
A
D
-
1
D
=
∑
n
Δ
X
m
n
Δ
X
m
n
T
ω
n
A
=
∑
n
Δ
x
m
n
Δ
x
m
n
T
ω
n
.According to (14) and (8), the optimal deformation gradient matrix and strain tensor of atom 146 at the time t
=
10 ps got by the method in [6] were(15)
F
^
175
=
0.9972
0.0107
-
0.1062
0.0042
1.0027
-
0.0255
-
0.0183
-
0.0096
1.0671
,
E
175
=
1.1102
e
-
016
-
5.4210
e
-
020
0
-
8.1315
e
-
020
1.1102
e
-
016
0
-
1.7347
e
-
018
-
8.6736
e
-
019
1.1102
e
-
016
.Then the effective strain of atom 175 wasε
e
f
f
175
=
3.3
×
1
0
-
16.The overall analysis method of parameter sensitivity based on the Latin Hypercube Sampling (LHS) method and Spearman rank correlation method proposed in this paper is useful to determine the range of parameters. According to the study in [8], the sensitivity of each parameter is related to its range and sampling times. Therefore, in order to balance the accuracy of each parameter, we should adjust the range of each parameter and the sampling times of all parameters to make all the parameter sensitivity values be close to each other and be less than 0.2. If a parameter sensitivity is much smaller than the sensitivities of other parameters, which means that the parameter accuracy is too high and not coordinated with the accuracy of the other parameters, then the computation quantity will be unnecessarily increased; if a parameter sensitivity is greater than 0.2, the accuracy of this parameter is insufficient.Considering the parameter sensitivity of each parameter in (7), the initial searching range of each parameter is set separately as F
i
j
∈
[
-
10,10
], k
1
∈
[
1.01,100
], and k
2
∈
[
0.01,20
]. By using the optimization method proposed in this paper with the weight function of (6), the atomic strain nephrograms at the time of t
=
10 ps, t
=
100 ps, and t
=
200 ps are shown in Figure 5. In Figure 5, the red dots and blue dots separately represent the largest and smallest strain of atoms.Figure 5
The atomic strain nephrograms.
(a)
t
=
10 ps
(b)
t
=
1
0
0 ps
(c)
t
=
20
0 psFigure5 can lead to the following conclusions.(1)
With the increase of time, the tensile displacement load grows linearly, and the maximum strain of atoms increases. These results are in accord with the theoretical situation.
(2)
When the load is very small, only a small part of the atomic strains are relatively large, as shown in a small amount of red dots in Figure5(a). With the increase of the load, plenty of metallic bonds break, and more and more atoms become disordered. The strains of these disordered atoms are relatively large, as shown in many red dots in Figures 5(b) and 5(c).At the timet
=
10 ps, the optimal deformation gradient matrix and strain tensor of atom 146 are(16)
F
^
175
=
-
0.9101
-
0.5269
-
0.1619
1.1913
1.1619
-
0.3205
0.3983
0.7830
1.9959
,
E
175
=
0
-
5.5511
e
-
017
-
5.5511
e
-
017
-
1.6653
e
-
016
-
5.5511
e
-
017
-
1.1102
e
-
016
-
6.9389
e
-
017
-
2.7756
e
-
017
0
.Then the effective strain of atom 175 isε
e
f
f
175
=
2.2
×
1
0
-
16, and the weight function is(17)
ω
n
=
2
+
44.5593
1
+
44.5593
-
1
1
+
44.5593
e
-
17.9041
r
.According to (17), the graph of weight function ω
n about the nondimensional parameter r can be got, as shown in Figure 6(a). Considering r
g
1
≈
0, then the change relation of the weight function ω
n with respect to the distance R
m
n between atom n and the center atom m can be got by (6), as shown in Figure 6(b).Figure 6
Theω
n-r and ω
n-R
m
n curves.According to Figure6(b), it can be found that the values of ω
n decrease from 1.0 to 0.0 with the increase of R
m
n. When R
m
n
→
r
c
u
t, ω
n
→
0; when R
m
n
>
r
c
u
t, ω
n
≈
0. Therefore, it is very convenient to calculate the objective function W
m
(
X
) shown in (5) by programming without judging whether the neighbor atom n is within the cutoff radius.


### 3.2. Weight Function Effects on Computed Strain

The comparison ofω
n-r curves (r
c
u
t
=
2 nm) of this paper and [6] is shown in Figure 7. The form and coefficients of the weight function ω
n play a key role on the atomic strain calculation. The coefficients of the weight function should be determined by (7) and cannot be arbitrarily assumed, so it is inappropriate for [6] to use the invariant weight function (see (13)) of the artificial assumption to calculate the discrete deformation gradient at the atom. According to Figure 7, it can be found that the optimized weight function curve is almost under the artificially assumed curve, so the effective atom strain got by the optimization method proposed in this paper is smaller and more reasonable than the results got in [6].Figure 7
Comparison ofω
n-r curves (r
c
u
t
=
2 nm) of this paper and [6].Then the effect of cutoff radiusr
c
u
t on the computed atom strain will be analyzed. When the cutoff radius is changed, the effective strain of atom 175 is shown in Table 1, and comparison of ω
n-r curves for different cutoff radii is shown in Figure 8.Table 1
The effect of cutoff radius on the effective strain of atom 175 in Figure4.
Cutoff radius
Atoms included
Effective strain of atom 175
k
1
k
2
1 nm
37
1.47 × 10−16
51.8432
16.5200
1.5 nm
84
3.33 × 10−16
59.9867
17.3910
2 nm
146
2.20 × 10−16
44.5593
17.9041
2.5 nm
223
3.40 × 10−16
87.0598
18.4370
3 nm
265
2.72 × 10−16
17.7613
19.2712Figure 8
Comparison ofω
n-r curves for different cutoff radii r
c
u
t by optimization method.According to Table1 and Figure 8, it can be found that when the cutoff radius varies from 1 nm to 3 nm, the optimized ω
n-r curves and effective strain of atom 175 are close to each other. When the nondimensional parameter r shown in (6) is more than 0.5, the effect of cutoff radius on the computed ω
n, or atom strain, can be almost neglected. The cutoff radius r
c
u
t was recommended to adopt the one of the potential functions.


## 3.1. Example Strain Calculation

In order to compute the strain tensors of atoms and provide evidence to confirm the method proposed in this paper to be correct and feasible by comparing with the result got by another method from the existing reference [6], the atomic mechanics model of alpha titanium material is built first. The length, width, and thickness of the model are separately 250 Å, 100 Å, and 3 Å, and the tensile displacement load on both upper and lower surface of the model is 0.025 Å/ps, as shown in Figure 4.Figure 4
The atomic mechanics model.The coordinate vectors of atoms in the reference configuration and in the current configuration are separately got at the time oft
=
0 and t
=
10 ps. Taking the center atom numbered 175 with a cutoff radius of 2 nm as an example, the number of atoms within the cutoff radius around the center atom is 146.The weight function in [6] was assumed as(13)
ω
r
=
1
-
6
r
2
+
6
r
3
r
≤
0.5
2
-
6
r
+
6
r
2
-
2
r
3
0.5
<
r
<
1.0
0
r
≥
1.0
.Then the optimal deformation gradient matrixF
^ could be got by (14)
F
^
=
A
D
-
1
D
=
∑
n
Δ
X
m
n
Δ
X
m
n
T
ω
n
A
=
∑
n
Δ
x
m
n
Δ
x
m
n
T
ω
n
.According to (14) and (8), the optimal deformation gradient matrix and strain tensor of atom 146 at the time t
=
10 ps got by the method in [6] were(15)
F
^
175
=
0.9972
0.0107
-
0.1062
0.0042
1.0027
-
0.0255
-
0.0183
-
0.0096
1.0671
,
E
175
=
1.1102
e
-
016
-
5.4210
e
-
020
0
-
8.1315
e
-
020
1.1102
e
-
016
0
-
1.7347
e
-
018
-
8.6736
e
-
019
1.1102
e
-
016
.Then the effective strain of atom 175 wasε
e
f
f
175
=
3.3
×
1
0
-
16.The overall analysis method of parameter sensitivity based on the Latin Hypercube Sampling (LHS) method and Spearman rank correlation method proposed in this paper is useful to determine the range of parameters. According to the study in [8], the sensitivity of each parameter is related to its range and sampling times. Therefore, in order to balance the accuracy of each parameter, we should adjust the range of each parameter and the sampling times of all parameters to make all the parameter sensitivity values be close to each other and be less than 0.2. If a parameter sensitivity is much smaller than the sensitivities of other parameters, which means that the parameter accuracy is too high and not coordinated with the accuracy of the other parameters, then the computation quantity will be unnecessarily increased; if a parameter sensitivity is greater than 0.2, the accuracy of this parameter is insufficient.Considering the parameter sensitivity of each parameter in (7), the initial searching range of each parameter is set separately as F
i
j
∈
[
-
10,10
], k
1
∈
[
1.01,100
], and k
2
∈
[
0.01,20
]. By using the optimization method proposed in this paper with the weight function of (6), the atomic strain nephrograms at the time of t
=
10 ps, t
=
100 ps, and t
=
200 ps are shown in Figure 5. In Figure 5, the red dots and blue dots separately represent the largest and smallest strain of atoms.Figure 5
The atomic strain nephrograms.
(a)
t
=
10 ps
(b)
t
=
1
0
0 ps
(c)
t
=
20
0 psFigure5 can lead to the following conclusions.(1)
With the increase of time, the tensile displacement load grows linearly, and the maximum strain of atoms increases. These results are in accord with the theoretical situation.
(2)
When the load is very small, only a small part of the atomic strains are relatively large, as shown in a small amount of red dots in Figure5(a). With the increase of the load, plenty of metallic bonds break, and more and more atoms become disordered. The strains of these disordered atoms are relatively large, as shown in many red dots in Figures 5(b) and 5(c).At the timet
=
10 ps, the optimal deformation gradient matrix and strain tensor of atom 146 are(16)
F
^
175
=
-
0.9101
-
0.5269
-
0.1619
1.1913
1.1619
-
0.3205
0.3983
0.7830
1.9959
,
E
175
=
0
-
5.5511
e
-
017
-
5.5511
e
-
017
-
1.6653
e
-
016
-
5.5511
e
-
017
-
1.1102
e
-
016
-
6.9389
e
-
017
-
2.7756
e
-
017
0
.Then the effective strain of atom 175 isε
e
f
f
175
=
2.2
×
1
0
-
16, and the weight function is(17)
ω
n
=
2
+
44.5593
1
+
44.5593
-
1
1
+
44.5593
e
-
17.9041
r
.According to (17), the graph of weight function ω
n about the nondimensional parameter r can be got, as shown in Figure 6(a). Considering r
g
1
≈
0, then the change relation of the weight function ω
n with respect to the distance R
m
n between atom n and the center atom m can be got by (6), as shown in Figure 6(b).Figure 6
Theω
n-r and ω
n-R
m
n curves.According to Figure6(b), it can be found that the values of ω
n decrease from 1.0 to 0.0 with the increase of R
m
n. When R
m
n
→
r
c
u
t, ω
n
→
0; when R
m
n
>
r
c
u
t, ω
n
≈
0. Therefore, it is very convenient to calculate the objective function W
m
(
X
) shown in (5) by programming without judging whether the neighbor atom n is within the cutoff radius.


## 3.2. Weight Function Effects on Computed Strain

The comparison ofω
n-r curves (r
c
u
t
=
2 nm) of this paper and [6] is shown in Figure 7. The form and coefficients of the weight function ω
n play a key role on the atomic strain calculation. The coefficients of the weight function should be determined by (7) and cannot be arbitrarily assumed, so it is inappropriate for [6] to use the invariant weight function (see (13)) of the artificial assumption to calculate the discrete deformation gradient at the atom. According to Figure 7, it can be found that the optimized weight function curve is almost under the artificially assumed curve, so the effective atom strain got by the optimization method proposed in this paper is smaller and more reasonable than the results got in [6].Figure 7
Comparison ofω
n-r curves (r
c
u
t
=
2 nm) of this paper and [6].Then the effect of cutoff radiusr
c
u
t on the computed atom strain will be analyzed. When the cutoff radius is changed, the effective strain of atom 175 is shown in Table 1, and comparison of ω
n-r curves for different cutoff radii is shown in Figure 8.Table 1
The effect of cutoff radius on the effective strain of atom 175 in Figure4.
Cutoff radius
Atoms included
Effective strain of atom 175
k
1
k
2
1 nm
37
1.47 × 10−16
51.8432
16.5200
1.5 nm
84
3.33 × 10−16
59.9867
17.3910
2 nm
146
2.20 × 10−16
44.5593
17.9041
2.5 nm
223
3.40 × 10−16
87.0598
18.4370
3 nm
265
2.72 × 10−16
17.7613
19.2712Figure 8
Comparison ofω
n-r curves for different cutoff radii r
c
u
t by optimization method.According to Table1 and Figure 8, it can be found that when the cutoff radius varies from 1 nm to 3 nm, the optimized ω
n-r curves and effective strain of atom 175 are close to each other. When the nondimensional parameter r shown in (6) is more than 0.5, the effect of cutoff radius on the computed ω
n, or atom strain, can be almost neglected. The cutoff radius r
c
u
t was recommended to adopt the one of the potential functions.


## 4. Conclusions

(1)
A common measure of deformation between atomic scale simulations and the continuum framework is provided and the strain tensors for multiscale simulations are defined in this paper. In order to compute the deformation gradient of any atomm, the weight function (see (6)) is proposed to eliminate the different contributions within the neighbor atoms which have different distances to atom m, and the weighted least squares error optimization model (see (7)) is established to seek the optimal coefficients of the weight function and the optimal local deformation gradient of each atom.
(2)
This paper proposed an overall analysis method of parameter sensitivity based on the Latin Hypercube Sampling (LHS) method and Spearman rank correlation method. The sensitivity of each parameter is related to its range and sampling times. In order to balance the accuracy of each parameter, we should adjust the range of each parameter and the sampling times of all parameters to make all the parameter sensitivity values be close to each other and be less than 0.2.
(3)
The advanced genetic algorithm proposed by this paper, based on the advanced niche genetic algorithm, suspicious peak value determination strategy, and local accurate searching techniques, effectively overcomes the drawbacks of basic genetic algorithm such as earliness and poor local searching ability, significantly improves global searching ability of genetic algorithm, and makes the solving process of the optimization model (see (7)) fast, accurate, and reliable.
(4)
The form and coefficients of the weight functionω
n play a key role in the atomic strain calculation, so the coefficients of the weight function should be determined for different materials by (7) and cannot be arbitrarily assumed.
(5)
When the nondimensional parameterr shown in (6) is more than 0.5, the effect of cutoff radius on the computed ω
n, or atom strain, can be almost neglected. The cutoff radius r
c
u
t was recommended to adopt the one of the potential functions.
(6)
The definition and optimization method of the atomic strain tensors for multiscale simulations proposed in this paper can be used for other nuclear power engineering materials.



---

*Source: 4649870-2016-08-18.xml*