Particle Modeling Based on Interatomic Potential and Crystal Structure: A Multi-Scale Simulation of Elastic-Plastic Deformation of Metallic Material

We formulate a macroscopic particle modeling analysis of metallic materials (aluminum and copper, etc.) based on theoretical energy and atomic geometries derivable from their interatomic potential. In fact, particles in this framework are presenting a large mass composed of huge collection of atoms and are interacting with each other. We can start from cohesive energy of metallic atoms and basic crystalline unit (e.g. face-centered cubic). Then, we can reach to interparticle (macroscopic) potential function which is presented by the analytical equation with terms of exponent of inter-particle distance, like a Lennard-Jones potential usually used in molecular dynamics simulation. Equation of motion for these macroscopic particles has dissipative term and fluctuation term, as well as the conservative term above, in order to express finite temperature condition. First, we determine the parameters needed in macroscopic potential function and check the reproduction of mechanical behavior in elastic regime. By using the present framework, we are able to carry out uniaxial loading simulation of aluminum rod. The method can also reproduce Young’s modulus and Poisson’s ratio as elastic behavior, though the result shows the dependency on division number of particles. Then, we proceed to try to include plasticity in this multi-scale framework. As a result, a realistic curve of stress-strain relation can be obtained for tensile and compressive loading and this new and simple framework of materials modeling has been confirmed to have certain effectiveness to be used in materials simulations. We also assess the effect of the order of loadings in opposite directions including yield and plastic states and find that an irreversible behavior depends on different response of the particle system between tensile and compressive loadings. How to cite this paper: Saitoh, K.-I. and Hanashiro, N. (2021) Particle Modeling Based on Interatomic Potential and Crystal Structure: A Multi-Scale Simulation of Elastic-Plastic Deformation of Metallic Material. World Journal of Nano Science and Engineering, 11, 45-68. https://doi.org/10.4236/wjnse.2021.113003 Received: July 6, 2021 Accepted: August 14, 2021 Published: August 17, 2021 Copyright © 2021 by author(s) and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access K.-I. Saitoh, N. Hanashiro DOI: 10.4236/wjnse.2021.113003 46 World Journal of Nano Science and Engineering


Introduction
Recently, multi-scale modeling of materials behavior with hierarchical approach has attracted much interest in research and development of materials science. In surveying various computational methodologies, it is recognized that molecular dynamics (MD) has been well established based on microscopic view [1], while traditional continuum-based approaches such as finite element (FE) analysis [2] [3] or particle methods [4] [5] (or sometimes called "mesh-less" methods) are formulated basically based on macroscopic viewpoint. In addition to these two separate approaches, in these days, a variety of combined methods between them has been and is being proposed extensively. For example, in the field of solid-state metallic material, it is guessed to be a very difficult task to model mechanical behavior with multi-scale/trans-scale (we would like to unify this conceptual word into "multi-scale" from now on) viewpoint, because such modeling has to contain both macroscopic and microscopic (atomistic) viewpoints simultaneously. As is often pointed out in the study of structural materials, discrepancy in space and time, when it comes from the atomistic level to the continuum level, is at least 10 0 -10 −10 meter (in space) or 10 0 -10 −15 second (in time) respectively. These differences are too wide to make concurrently combined computation model easily. Only in the case when the mechanical behavior of a solid material can be limited to just small strain or small deformation regime, it may become quite possible by constructing the model with elastic and linear constitutive law and parameters. Such approaches include quasi-continuum method [6], or perfectly macroscopic particle methods such as smoothed particle hydrodynamics (SPH) method [7] and Moving Particle Simulation (MPS) method [8]. In recent years, the new approach called "Peridynamics" (PD) [9] becomes the most popular particle method and is being focused on in the research field of solid mechanics. In the theory of PD, atomistic and non-local conceptual view is brought into and embedded on the formulation of macroscopic framework and a good affinity with MD is expected. But, due to its innovative but developing nature, it seems to be still on the way of waiting the time when it will be affirmatively applied to multi-scale problems in materials science and engineering.
In conventional approach so far to construct a multi-scale view including both the atomistic and the continuum, there is a try that potential energy function between atoms in solid (metal and crystal) is directly connected to macroscopic elastic constants. In that case, a harmonic approximation that is valid only in small displacement of elastic regime should be assumed. However, for plastic World Journal of Nano Science and Engineering deformation appearing in much larger loading level, the simple approximation rather becomes quite hard task to link directly between microscopic behavior (i.e. slip motion between atoms or atomic diffusion) and macroscopic constitutive law. Plasticity, super-elasticity as well, which is usually observed in many solid materials, is one of difficult issues, and so adequate multi-scale modeling in theory and computational set-up should be required.
Nowadays, it is recognized that not only macroscopic modeling but also microscopic (atomistic) one is inevitable for further development of engineering materials [10]. We have known that, in fact, total mechanical behavior of materials is controllable by introducing state-of-art nanoscale structures. Thus, what is required at this point is a feasible methodology capable of analyzing solid and crystal materials, in which microscopic and macroscopic methods are coupled. In recent years, particle methods (or sometimes called "mesh-less" methods), as an auxiliary method for on-mesh methods such like FE method, have been developed and used extensively, by virtue of the enhancement of computation power. In those new methods, "particles" are likely to be treated in various senses, which correspond to their size scales. On the smallest scale, particles should be real substances as atoms or molecular groups (which is as in MD method). On the other hand, in larger scale, particles are only representative points on continuous body of the material, which are needed for the approximation of discrete-type computation. In every case of particle methods, an essential feature, on the whole, is that we are solving equations of motion and will chase the motion of particle step by step in the time development. Therefore, we can expect that the combination between different scales in particle methods would be nicely accomplished by unified theory along with adequate modeling using particles.
The authors recognize that there are two approaches having opposite direction in modeling with particles as follows. One approach is that macroscopic variables are immediately connected to the evaluation of energies and deformations obtained in atomistic simulation. For example, stress and strain are absolutely macroscopic evaluations that are only defined under continuum hypothesis and relation, but they can be transferred into atomic simulations with the name of "atomic stresses" [11] or "atomic strains" [12] [13]. Another approach is that, to the contrary, an atomistic (microscopic) variable is utilized in macroscopic modeling for materials behavior. In fact, this idea has already been implemented as some new computational theories, such as quasi-continuum (QC) method [14] [15], virtual atomic cluster (VAC) method [16], and so on [17]. In our opinion, the latter approach seems challenging and interesting for materials modeling. This idea also leads to a new trial, where an interatomic potential function that is directly derived from microscopic relation can be conversely applied to macroscopic particle modeling (MPM). The equations of macroscopic particles to be solved are of the same type as those in atomic system, say, in MD method, while the size of particles, i.e. the interparticle distance, is not equivalent to atomic size but is much bigger one. World Journal of Nano Science and Engineering In this paper, we propose a new framework of MPM method stated above. This MPM will be constructed on microscopic (atomic) potential energetics and dynamics which are obtained from atomistic information. This method is capable of reproducing macroscopic materials behavior by adopting Langevin-type equation of motion, where dissipative force, random force indicating thermal fluctuation, and conservative potential force derived from atomic system are well integrated. In the context of our study, the conservative potential force is being expressed by a simple power-law relevant in atomistic simulation of condensed matter, such as the Lennard-Jones potential [18]. We will see that a simple expression is sufficiently useable to obtain macroscopic materials behaviors. By using this method, we can prospect that mechanical properties and thermal properties of solid (metallic) materials will be evaluated very simply. Moreover, it can lead to more sophisticated multiscale approach in the future. In this paper, as the first stage of the study, we would like to exhibit the detail in derivation of equation and parameters settings. Then, we will verify it through computational results of simple uniaxial loading in elastic-plastic regime.

Particle Modeling Method: Basic Concept
Generally, in particle modeling, the material which occupies a certain space is replaced by assembly of discrete particles. This concept is universal and common to that in the former literature [19] [20] [21] [22] and is also used in usual SPH (smoothed particle hydrodynamics) method [7] [23] [24] [25] [26] [27]. For example, when total volume is confined in a cube with dimension L L L V × × = as shown in Figure 1, it is replaced by N particles. These particles have adequate mass, so that the original density of the material ρ is retained. Motion of each particle is described by Newtonian equation of motion, where ( )  From the former study of the original "particle modeling" method [19] [20], a power function expressed by has been used for interparticle force. In atomic simulation, the famous Lennard-Jones potential has the same function form as Equation (4) (corresponding to p = 12, q = 6). These macroscopic potential parameters, G, H, p, q, have to be determined. The units of parameters, G and H, will be determined after two non-dimensional parameters, p and q, are determined so as that ij F in equation

Introducing of Langevin Equation for Particle Modeling
In atomic system, Equation (2) is supposed to have only conservative force so that the total energy of the system should be conserved in principle. However, macroscopic particles' system should have additional treatment due to invisible effect from many disappeared degrees of freedom. Therefore, Langevin-type equation of motion can be applied to this particles' system, where energy dissipation and production of thermal energy by fluctuation are integrated. This type of equations have already used to dissipative particle dynamics (DPD) [28] and Brownian dynamics (BD) [29] simulations, which can analyze material behavior in slightly larger scale than atomic size. The Langevin equation is given by, World Journal of Nano Science and Engineering where , i D F and , i R F are force vectors called dissipative and random forces, respectively. The conservative force on one particle is derived from which is total conservative energy of the system. The dissipative force term can be formulated by considering inherent resistance of a moving particle. In the meso-or micro-scale (relatively in smaller scale), that term is supposed to be caused by interparticle friction. When a particular viscosity μ is provided, the dissipative force should be estimated, from particles' velocity On the other hand, in meso-and macro-scale (relatively in larger scale), the motion of particle should be through some other media, it is supposed that the dissipative force depends only on absolute velocity. In such a case, the friction force may be estimated by Random force induces fluctuation of particles and is formulated based on its randomness.
Thus, in computation, numerically produced random numbers (pseudo-random numbers) with Gaussian distribution are utilized. The dissipative and random forces are working for a particle motion, so to speak, as brake and accelerator, respectively. All these three types of force in Langevin equation (Equation (6)) balance each other in equilibrium, or they evolve to make structural change (deformation and catastrophic fracture) in finite temperature condition.

Determination of Potential Parameters and Mechanical Properties
In the present MPM method, conservative force almost determines elastic properties. Therefore, based on existing elastic constants available from experiment or ab initio (first principle quantum chemistry) calculation, macroscopic potential parameters are determined. At first, total absolute amount of potential energy contained in one macroscopic particle is estimated. When atoms are condensed as solid material with density ρ (its unit is kg/m 3 ), a lot of atoms are assigned to one macroscopic particle. For example, in solid state of fcc metal, each atom has cohesive energy e 0 which is calculated as the summation of twelve pairwise potential energies (its unit is J). If the number of atoms in one macroscopic particle is N atom , collective energy in this macroscopic particle should be, when all the macroscopic particle have a spherical shape and they gather to make fcc structure of lattice constant a 0 (its unit is meter), N atom is given by, where r 0 is equilibrium distance between macroscopic particles (its unit is m), m is atomic mass (its unit is kg), the factor 0.74 means occupancy rate in fcc lattice, which becomes non-dimensional. In most cases, e 0 and a 0 as well as structural type of crystal are easily available from experimental fact or ab initio evaluation. Therefore, these e 0 , a 0 and ρ are recognized as microscopic parameters, but N atom and P atom come to macroscopic variables to be configured next. When macroscopic particles are separated by equilibrium distance 0 r r = , their interparticle force should all vanish, and Equation (4) gives the condition, Equivalence between microscopic and macroscopic energies (Equation (10)) and equilibrium condition (Equation (12)) are both used to solve unknown potential parameters G, H by following equations.
Two undetermined exponents (powers) p, q which determine net shape of potential curve are needed and will be determined next. Thus, energetic link between microscopic and macroscopic systems has been carried out. However, for mechanical response, the curve shape of potential function is crucial. Although there may be many choices for function type, we think that Equation (5) is reasonable because it furnishes three basic characteristics necessary for macroscopic particles. They are: 1) divergent feature of energy in closer separation than r 0 ; 2) equilibrium distance (energy minimum) at r 0 ; and 3) convergence to zero-energy in far larger separation. The undetermined parameters p, q can be adjusted from elastic constant (i.e. Young's modulus, in the unit of Pa) in the vicinity of equilibrium distance r 0 as follows.
Young's modulus E can be estimated as the curvature of potential curve at equilibrium distance r 0 , then, Substituting Equation (5) into this formula, it gives, This means that Young's modulus depends on the product between p and q. For the fcc metals (aluminum, copper, etc.), using each cohesive energy e 0 , Young's modulus is estimated. They are tabulated on dependency on p, q as shown in p = 4 and q = 9 may be the most suitable. Of course, we have another choice, but by this procedure we are able to set p, q to configure elastic response of macroscopic particles system.

Introducing Plasticity in the Present Particle Modeling Framework
It is generally difficult to express both microscopic and macroscopic plasticity mechanisms at the same time. In macroscopic view, constitutive relation (sometimes it exhibits very complicated formula) can be responsible for plasticity. But, it contains lots of macroscopic parameters which are not immediately connected to microscopic dynamics or parameters.
Energy dissipation or heat transport in plasticity, for example, has been one of difficult issues. We are challenging to introduce plasticity mechanism into macroscopic particles' system based on microscopic parameters, extending elastic method described above. Figure 2 shows potential function including plastic regime. U ϕ is the function representing "unloading behavior" after plastic deformation proceeds, which is expressed by

Coarsening Dissipative Force and Random Force Terms
In the field of MD study, the Langevin equation is sometimes used and is including Debye model for thermal conduction [30]. Dissipative term is formulated by, and random force is given by where α is a viscosity for atom and σ is the square root of variance. In MD scale, α and σ have the relation, and where D B k h ω θ = is Debye frequency ( θ , h are Debye temperature and Planck's constant, respectively), k B is Boltzmann's constant, Δt is time increment of MD and T is an equilibrium temperature. This MD system has to be coarsened and transformed to macroscopic particles' system, by virtue of scale factors. There is a trend of study to link microscopic parameters to macroscopic relation [31] [32]. However, here, to start with, we assume very simple relations as follows.
The scale factors are CG η for α and CG ζ for σ, and resulted expressions are Verifying these coarsening parameters, CG η and CG ζ , should be required, but it will be done in further studies.
In this paper, as a first stage of this study, we think that the potential force (conservative force) is the most important factor for evaluation of material properties in solid. Therefore, we can omit those dissipative and random forces temporarily. It means that MPM simulation results discussed in the following section are carried out in zero Kelvin condition.
However, in principle, the MPM simulation in finite temperature T is also feasible by the basic framework shown in this section.

Arrangement of Macroscopic Particles for the Cylindrical Specimen
In this study, it is assumed that a metallic material is subjected to the external loading. The specimen is composed of pure aluminum and cylindrical shape. It is applied uniaxial tensile or compressive loading along the longitudinal direction. The macroscopic particles are arranged in face-center-cubic (fcc) lattice as shown in Figure 3 Only two centered particles at each end part are fixed in their initial position also in x and y directions as well, so as to realize uniaxial loading to the specimen. The longitudinal strain ε is estimated from updated distances between particles at both ends. Stress tensor components are estimated from virial formula, which is conceptually an average of the product of interparticle forces and difference of position vectors of all interacting particles inside the system. This evaluating method of stress is usually used in molecular dynamics calculation, which deals with atomistic particles' system.

The Choice of Height and Diameter of Cylindrical Models.
The cylindrical specimen with the same diameter of 14 mm made of aluminum is being divided neatly into macroscopic particles. In dividing, we designate two integers, n ϕ and h n , which are numbers of division in the direction of diameter (width) and height, respectively. These numbers of division are chosen as the power of 2 as expressed in Equations (23) and (24) using other integers i and j, and they are excluding cases of the division number less than 2.
( ) The numbers of division and calculation conditions are summarized in Table   2. The number of particles and resulted actual height of the specimen depending on the choice of n ϕ and h n are shown in Table 3. The arrangement of particles (view of external form and at cross-section) when the ratio of division parameters ( h n n φ ) is 2 2 4 = are shown in Figure 4.

Results of Tensile Simulations
Stress-strain curves obtained for ( ) 0 2 1 h n n φ = = are shown in Figure 5.  n φ is 2 1 (=2), the largest modulus is obtained. As shown in Figure 6 at the same time, the dependency on division number n φ of the strain of the specimen which is evaluated from the change of interparticle distances exhibits quite similar trend to Young' moduli. Estimation of the strain tends to be larger for smaller division like 2 n φ = , and it means that in such cases a rough division by particles results in larger stiffness than with finer division. It is confirmed that another division parameter in the direction of the specimen's height, h n , also shows the same tendency of varying stiffness as explained for n φ .
The dependence of Poisson's ratio ν on the division ratio h n n φ is shown in Figure 7. When n φ increases, Poisson's ratio saturates and reaches to 0.33, which is almost the same as the actual value of pure aluminum.
Next, stress distributions comparing between strain value 0% and 0.1% are shown in Figure 8, for the case of a division ratio 4 h n n φ = . The pictures are drawn at cross-sections at 0 x = , so as to see inside of the specimen. They are normal stress component in tensile direction (z), where color ranges from green to red according to tensile stress value. As shown in Figure 8, inside of the specimen almost uniform distribution of stress is obtained, but particles at the surface exhibit smaller value which is caused by reduction of the coordination number of the particles there.

Results of Compression Simulation
The set-up of compression simulations is the same as that of tensile testing, except for the moving direction of end regions. The dependence of Young's modulus on the numbers of division , h n n φ and their ratio h n n φ in these compression testing is shown in Figure 9. The tendency is almost the same as the tensile case, but it is slightly lower for the compressive condition. It is because in the tensile case particles at the surface tend to be constrained in longitudinal direction, but in the compression case the motion toward free surface are less constrained and quite allowable and thus the stress at the surface shows lower values.   The strain values when compressive strain is 0.1% are compared as in Figure   10(a). These strain values are slightly smaller than those in tensile testing (not shown here). The averaged strain ranges from 0.30 to 0.35, which is affected by initial configuration of particles, that is, fcc lattice structure as shown Figure   10(b). The ratio of the local strain between neighbor particles, ij ε , to the global strain as for the specimen, system ε , is estimated approximately to be The stress distribution inside the specimen for the ratio between the numbers of division, 4 h n n φ = , is shown in Figure 11. The component of the normal stress in compressive direction is shown and its color ranges from green to blue according to the compressive value. Inside the specimen, almost uniform distribution is observed, but particles at the surface show smaller compressive stress (as absolute value).

Results for Simple Tensile Elastic-Plastic Loading
The tensile simulation including not only elastic regime but also plastic one is conducted. The specimen is the same as that used in elastic simulation, which is built with particle division parameters  Table 4.
The obtained stress-strain curves including plastic regime is shown in Figure   12. After yielding (approximately 1% of nominal strain), stress still increases as   strain is increasing, which means strain hardening occurs. The stress distribution when the specimen is committing strain hardening is shown in Figure   13(a). The color of particles ranges continuously from green to red in relative manner, where red-colored particles have larger tensile stress. Figure 13  inside the specimen increases almost uniformly, but it is always lower than those in surface particles.

Results for Simple Compressive Elastic-Plastic Loading
The compression testing including not only elastic regime but also plastic one is also conducted. The obtained stress-strain curves for the selected division case are shown in Figure 14. The stress distribution obtained at each strain is shown Figure 15(a). The color ranges from green to blue, where blue means larger stress value. The increment of stress during change per −0.5% strain increase is also shown in Figure 15(b).
During the period of strains from 3.0%    segment does compressive forces. In tensile simulation, component of bonding force between neighbor particles which is parallel to the loading direction is responsible for plastic strain, but that perpendicular to tensile direction becomes compressive. To the contrary, in compressive simulation, the direction of bonding force is reversed. Due to the formulation of plasticity we have established, only tensile strain accompanies stress relaxation in yielding event. In compressive deformation as shown in Figure 16(b), component of bonding force perpendicular to loading axis exhibits yielding and shearing force also occurs and it causes slip between particles.

Consideration on Bachinger's Effect (the Case of Reversing
Loading Directions) It is experimentally known that metals will usually show a kind of Bauchinger's effect, where reduction of yield stress occurs when the direction of applied plastic stress or strain is reversed. But, simple constitutive modeling in numerical simulation is formulated with a fixed yield condition, and it is agreed that very complicated theoretical formulation must be required for reproducing Bauchinger's effect. At this point, we will see whether Bauchinger's effect is reproduced by using the present particle modeling method.
First, uniaxial tensile strain up to 3% is applied to the cylindrical specimen, and then subsequently uniaxial compressive strain up to 3% is applied to the same specimen. The sequence of stress-strain diagrams is shown in Figure 17. Figure 17(a) is showing overall process, in that solid line means continuous tensile and compression processes and dashed line is for only compression simulation from unloaded state, for comparison. Figure 17(b) is the view magnifying around the yielding (that strain is approximately 1.2%). When comparing two lines, the case of compression from unloaded state exhibits yield stress smaller than that obtained in reversed deformation. Specifically, at compressive strain −0.3%, they are −0.575 GPa (solid line) and −0.552 GPa (dashed line), which means the contrary behavior to Bauchinger's effect. This is because yielded state provided by tensile process is just realized by stress shift due to mathematical function, and is not accompanying any topological change of particles' configuration such as slip motion. The arrangement of particles has not disturbed during tensile plastic deformation. But, in turn, the next compressive yield event and subsequent plastic deformation must accompany a topological change (that is, slip motion) of particles' configuration, and they should require larger force or energy than simple compressive plasticity.
To the contrary, the reversed case, in that compressive load is first applied then tensile one is applied, is resulted as shown in Figure 18. The magnified view around the yield event is shown in Figure 18(b). The secondary tensile loading (solid line) shows slightly smaller value of yield stress than the case from undeformed state (dashed line). The specified values at tensile strain 0.3% are 0.574 GPa (solid line), and 0.588 GPa (dashed line). It means Bauchinger's effect. The difference of yield stress in this case is explained in the same way as the case of first tensile and secondary compression which was stated above.  These simulations with reversed deformation display that any symmetrical yielding behavior is not included in the present formulation which uses shift of potential energy between neighbor particles. It is concluded that an accurate repetitive yielding behaviors for cyclic loading and Bauchinger's effect must take non-symmetrical interaction between particles into account and the formulation of potential function should be reorganized. These considerations will be studied further as our future work.

Conclusion
In this study, we propose a computational framework by using macroscopic particles method (MPM). Although the proposed MPM method is simply constructed from microscopic parameters such as cohesive energy and density of metallic materials, it provides high feasibility on implementing a simple frame-World Journal of Nano Science and Engineering work for multi-scale simulation and in discussing the deformation of metallic materials. As seen in the examples of rod-shaped models we presented here, the evaluation by the present MPM method depends on the division number of particles. In the elastic regime, Young's modulus and the Poisson's ratio can be sufficiently reproduced. On the other hand, it is understood that the plastic modeling in MPM method still remains the difficult issue. The plasticity mechanism is absolutely important for material modeling, and as the first challenge, we have discussed an irreversible change of the configuration of particles. We found that there is different stability in between tensile and compressive loadings for the total elastic-plastic regime. We can conclude that further sophistication and improvement of the MPM method especially for plastic mechanism are quite expectable and they are worthwhile studying continuously in the future.