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 atom , the weight function is proposed to eliminate the different contributions within the neighbor atoms which have different distances to atom , 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.
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 time and at the current time , we assume a fixed Cartesian coordinate system as shown in Figure 1.
In Figure 1, and are, respectively, the coordinate vectors of atoms in the reference configuration at the initial time and in the current configuration at the current time , and maps the atoms from to . Consider
Assuming sufficient continuity, the local deformation at is characterized as the gradient of the motion and can be defined as
The deformation of an infinitesimal segment at one point in the reference configuration can be expanded in a Taylor series. If the higher-order terms are omitted, can be written as
Then the deformation gradient 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 in the reference configuration and position in the current configuration . Then the relative position of neighboring atom and the deformation gradient 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 of atom cannot generally be got by a single atom . Therefore, we should seek an optimal local deformation gradient , 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 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 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 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 aswhere and are undetermined coefficients, is the distance between atom and atom , and 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 to satisfy the objective function as shown in
When the objective function reaches the minimum, the parameters () are the optimal components of the optimal deformation gradient matrix .
When is determined by the optimization model as shown in (7), the Lagrangian Green strain tensor is defined with respect to reference coordinates in terms of this quantity as
The strain tensor is a 3 × 3 matrix, and the effective strain of a point can be described by the second invariants of , 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].(1)Build systematic model, meaning function relation between system feature and factors, , and get the benchmark sets of parameters.(2)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.
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 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 and pick times in sequence, and then one will get sample parameter sets .
The realization process schematic drawing of LHS method is shown in Figure 2.
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 from large to small, an ordered sequence is formed. Supposing is at the order of in the sequence , then is the rank of in .(2)By substitution of sets of sample parameters into (7), the solution sets of objective function can be obtained.(3)Supposing the rank of the parameter in is and the rank of the parameter in 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.
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 amount of all the suspicious peak value point must be the global optimization point, which is referred to as objective function , and the corresponding variables are . Set the objective function of a suspicious peak value point as , and the corresponding variables are . The condition that point is another global optimization point different from iswhere 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. Supposing the variable sets of one global optimization point are and the range of each variable is , the new searching range of each variable will be changed towhere 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.
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.
The coordinate vectors of atoms in the reference configuration and in the current configuration are separately got at the time of and 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 could be got by
According to (14) and (8), the optimal deformation gradient matrix and strain tensor of atom 146 at the time ps got by the method in [6] were
Then the effective strain of atom 175 was .
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 , , and . By using the optimization method proposed in this paper with the weight function of (6), the atomic strain nephrograms at the time of ps, ps, and ps are shown in Figure 5. In Figure 5, the red dots and blue dots separately represent the largest and smallest strain of atoms.
(a) ps
(b) ps
(c) ps
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 ps, the optimal deformation gradient matrix and strain tensor of atom 146 are
Then the effective strain of atom 175 is , 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 , 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 , ; when , . Therefore, it is very convenient to calculate the objective function shown in (5) by programming without judging whether the neighbor atom is within the cutoff radius.
3.2. Weight Function Effects on Computed Strain
The comparison of - curves ( 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 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. The cutoff radius 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 atom , the weight function (see (6)) is proposed to eliminate the different contributions within the neighbor atoms which have different distances to atom , 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 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 parameter shown in (6) is more than 0.5, the effect of cutoff radius on the computed , or atom strain, can be almost neglected. The cutoff radius 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by the National Defense Basic Scientific Research Program of China through the Contract of B1520132013-1.