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

A common measure of deformation between atomic scale simulations and the continuum framework is provided and the strain tensors formultiscale 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 atomm, 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.


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.

Deformation Gradient and Strain Tensors.
In order to describe the positions of atoms at the initial time  0 and at the current time , we assume a fixed Cartesian coordinate system as shown in Figure 1.
In Figure 1, X and x are, respectively, the coordinate vectors of atoms in the reference configuration Ω 0 at the initial time  0 and in the current configuration Ω 1 at the current time , and  maps the atoms from X to x.Consider Assuming sufficient continuity, the local deformation at X is characterized as the gradient of the motion and can be defined as The deformation of an infinitesimal segment x at one point in the reference configuration can be expanded in a Taylor series.If the higher-order terms are omitted, x can be written as Then the deformation gradient F 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  is characterized by the changes in the relative position of its neighbors.Atom  is located at the position X  in the reference configuration Ω 0 and position x  in the current configuration Ω 1 .Then the relative position of neighboring atom  and the deformation gradient F  at atom  are given as The deformation gradient of atom  is related to its neighboring atoms, and the deformation gradient of each neighboring atom with respect to atom  is different, so the deformation gradient F  of atom  cannot generally be got by a single atom .Therefore, we should seek an optimal local deformation gradient F , which can make the weighted squares error   shown in the following equation minimized: where  is the number of neighboring atoms and   is a weight factor.The optimal deformation gradient F of the atom  is determined by the deformation gradients of all the neighboring atoms with respect to atom , but the contribution of each atom  to F is different and related to the distance between atom  and atom .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  cut that specifies the domain around atom ,   of the nearest neighbor atom is equal to 1.0 and gradually reduced to 0.0 with the increase of the distance between atom  and atom .The weight function   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   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   can be assumed as where  1 and  2 are undetermined coefficients,   is the distance between atom  and atom , and  1 is the distance between the closest neighbor atom and atom .
The purpose of parameters identification is to find a set of solutions of discrete design variables X to satisfy the objective function   (X) as shown in When the objective function   (X) reaches the minimum, the parameters F  (,  = 1, 2, 3) are the optimal components of the optimal deformation gradient matrix F .
When F 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 The strain tensor E is a 3 × 3 matrix, and the effective strain of a point  eff can be described by the second invariants of E, shown as 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].
( 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   (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 of  parameters, and set the value range of each parameter as  nonoverlap zones with equal probability.
(2) Randomly pick a sample of every parameter   within its range, meaning each parameter   consists of  sample values.
(3) Randomly arrange the  sample values of each   , to form  random arrangements.
(4) Pick up one sample value for each parameter from the  arrangements, which could form a sample parameter set X  and pick  times in sequence, and then one will get  sample parameter sets {X 1 , X 2 , . . ., X  }.
The realization process schematic drawing of LHS method is shown in Figure 2.

Science and Technology of Nuclear Installations
Sets of sample parameters: Sampling and random permutation of parameter x 1 : Sampling and random permutation of parameter x 2 : Sampling and random permutation of parameter x i : Sampling and random permutation of parameter x m : 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 { 1 ,  2 , . . .,   } from large to small, an ordered sequence  is formed.Supposing   is at the order of   in the sequence , then   is the rank of   in { 1 ,  2 , . . .,   }.
(3) Supposing the rank of the parameter    in { 1   ,  2  , . . .,    } is   and the rank of the parameter   in { 1 ,  2 , . . .,   } is   , then the sensitivity of parameter   equals the absolute value of Spearman rank correlation coefficient, as shown in Arrange the sensitivity degree   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.

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  amount of all the suspicious peak value point must be the global optimization point, which is referred to as objective function where   is the searching range around the variable   .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 Figure 3.

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.
The coordinate vectors of atoms in the reference configuration and in the current configuration are separately got at the time of  = 0 and  = 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 Then the optimal deformation gradient matrix F could be got by According to ( 14) and ( 8), the optimal deformation gradient matrix and strain tensor of atom 146 at the time  = 10 ps got by the method in [6] were Then the effective strain of atom 175 was  175 eff = 3.3×10 −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  ∈ [−10, 10],  1 ∈ [1.01, 100], and  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  = 10ps,  = 100 ps, and  = 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 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 Figure 5(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 time  = 10 ps, the optimal deformation gradient matrix and strain tensor of atom 146 are Then the effective strain of atom 175 is  175 eff = 2.2 × 10 −16 , and the weight function is According to (17), the graph of weight function   about the nondimensional parameter  can be got, as shown in Figure 6(a).Considering  1 ≈ 0, then the change relation of the weight function   with respect to the distance   between atom  and the center atom  can be got by (6), as shown in Figure 6(b).
According to Figure 6(b), it can be found that the values of   decrease from 1.0 to 0.0 with the increase of   .When   →  cut ,   → 0; when   >  cut ,   ≈ 0. Therefore, it is very convenient to calculate the objective function   (X) shown in (5) by programming without judging whether the neighbor atom  is within the cutoff radius.

Weight Function Effects on Computed Strain.
The comparison of   - curves ( cut = 2 nm) of this paper and [6] is shown in Figure 7.The form and coefficients of the weight function   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].
Then the effect of cutoff radius  cut 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   - curves for different cutoff radii is shown in Figure 8.
According to Table 1 and Figure 8, it can be found that when the cutoff radius varies from 1 nm to 3 nm, the optimized   - curves and effective strain of atom 175 are close to each other.When the nondimensional parameter  shown in (6) is more than 0.5, the effect of cutoff radius on the computed   , or atom strain, can be almost neglected.

Figure 1 :
Figure 1: General motion in the neighborhood of a discrete atomic particle.

Figure 3 :Figure 4 :
Figure 3: The program flow diagram of advanced GA.

Figure 8 :
Figure 8: Comparison of   - curves for different cutoff radii  cut by optimization method.
When analyzing the effect of one particular parameter   on feature , keep the rest of parameters constant in basic value, and only allow parameter   to vary within the proper range.If minor change of   leads to drastic changes of , that means  is very sensitive to   , and the sensitivity value of is large; if   varies within a rather wide range, while the change in  is not obvious, that means  is not sensitive to   , and the sensitivity value of   is small.
, and the corresponding variables are { * 1 ,  * 2 , . . .,  *  }.Set the objective function of a suspicious peak value point  as , and the corresponding variables are { 1 ,  2 , . . .,   }.The condition that point  is another global optimization point different from  is Supposing the variable sets of one global optimization point are { * 1 ,  * 2 , . . .,  *  } and the range of each variable is   ≤   ≤   , the new searching range of each variable will be changed to *  ∑ =1 (  −  *  ) 2 >   ( = 1, 2, . . ., ) , (11b) where  (0 <  < 1) is a constant, meaning the relative searching range around the global optimization value;   is a constant, which is used to judge whether two points are neighboring; and   is the Hamming distance between two points.If   ≤   , it means that points  and  are neighboring, and they are essentially the same global optimization point with different accuracy.Equation (11a) is used to determine that  is also a global optimization point, and (11b) is used to determine that  is another global optimization point different from . (3) Local Accurate Searching for All the Different Global Minimum.

Table 1 :
The effect of cutoff radius on the effective strain of atom 175 in Figure4.