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

Ken-Ichi Saitoh, Naoya Hanashiro

Department of Mechanical Engineering, Faculty of Engineering Science, Kansai University, Suita, Japan.

DOI: 10.4236/wjnse.2021.113003                

Abstract

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.

Keywords

Molecular Dynamics, Particle Method, Elasticity, Plasticity, Nanostructure, Multi-Scale Modeling

1. 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⁰ - 10⁻¹⁰ meter (in space) or 10⁰ - 10⁻¹⁵ 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 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.

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.

2. Theory and Method

2.1. 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\times L\times L=V$ as shown in Figure 1, it is replaced by N particles. These particles have adequate mass,

$M_i=\rho \frac{V}{N}$, (1)

so that the original density of the material ρ is retained. Motion of each particle is described by Newtonian equation of motion,

$M_i\frac{{\text{d}}^2{r}_i}{\text{d}{t}^2}=F_i$, (2)

where $F_i$ includes interparticle force as well as additional forces induced by viscous friction and thermal fluctuation. Indeed, this equation is the same as that of molecular dynamics (MD) simulation. Interparticle interaction can be configured from pairwise potential function ${\phi }_E\left(r\right)$, and then the conservative force is derived by

$F_i={\displaystyle {\sum }_{j=1}^{Neighbor}{F}_{ij}}={\displaystyle {\sum }_{j=1}^{Neighbor}\frac{\partial {\phi }_{E,ij}}{\partial r_{ij}}\frac{r_{ij}}{r_{ij}}}$, (3)

where ${\phi }_{E,ij}={\phi }_E\left(r=r_{ij}\right)$ means potential energy acting between two particles i and j. Separation between particles $r_{ij}$ is the absolute value of interparticle difference vector $r_{ij}=r_j-r_i$.

Figure 1. Concept of macroscopic particle modeling (MPM).

From the former study of the original “particle modeling” method [19] [20], a power function expressed by

$F_{ij}=\frac{G}{r_{ij}^{p+1}}-\frac{H}{r_{ij}^{q+1}}$ (4)

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 $F_{ij}$ in equation (4) should have the unit of force, i.e. Newton (N). From energy-conservative mechanics, a potential function integrated from force function Equation (4) is given by

${\phi }_{E,ij}=-\frac{G}{p{\left(r_{ij}\right)}^p}+\frac{H}{q{\left(r_{ij}\right)}^q}$ . (5)

2.2. 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,

$M_i\frac{{\text{d}}^2{r}_i}{\text{d}{t}^2}=F_{i,C}+F_{i,D}+F_{i,R}=-\frac{\partial \Phi }{\partial r_i}+F_{i,D}+F_{i,R}$, (6)

where $F_{i,D}$ and $F_{i,R}$ are force vectors called dissipative and random forces, respectively. The conservative force on one particle is derived from

$\Phi ={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{i>i+1}^N{\phi }_{E,ij}}}$, (7)

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 $v_i$, $v_j$, by

$F_{i,D}={\displaystyle {\sum }_{j=1}^{Neighbor}\mu \left(v_j-v_i\right)}$. (8)

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

$F_{i,D}=\gamma v_i$. (9)

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.

2.3. 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 abinitio (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³), a lot of atoms are assigned to one macroscopic particle. For example, in solid state of fcc metal, each atom has cohesive energy e₀ 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,

$P_{\text{atom}}=N_{\text{atom}}{e}_0$. (10)

when all the macroscopic particle have a spherical shape and they gather to make fcc structure of lattice constant a₀ (its unit is meter), N_atom is given by,

$N_{\text{atom}}=\frac{\rho \pi }{6}\frac{r_0^3}{m}\frac{1}{0.74}$, (11)

where r₀ 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₀ and a₀ as well as structural type of crystal are easily available from experimental fact or abinitio evaluation. Therefore, these e₀, a₀ 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 $r=r_0$, their interparticle force should all vanish, and Equation (4) gives the condition,

$F_{ij}\left(r=r_0\right)=\frac{G}{r_0^{p+1}}-\frac{H}{r_0^{q+1}}=0$. (12)

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.

$G=f\left(p,q\right)r_0^p,H=f\left(p,q\right)r_0^q,f\left(p,q\right)=\frac{P_{\text{atom}}}{6}\frac{pq}{p-q}$. (13)

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₀; 2) equilibrium distance (energy minimum) at r₀; 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₀ as follows.

Young’s modulus E can be estimated as the curvature of potential curve at equilibrium distance r₀, then,

$E=\frac{1}{2r_0}\frac{{\partial }^2{\phi }_{E,ij}}{\partial r^2}\left(r=r_0\right)$. (14)

Substituting Equation (5) into this formula, it gives,

$E=\frac{\pi \rho \frac{1}{0.74}pq}{72}\frac{e_0}{m}$. (15)

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₀, Young’s modulus is estimated. They are tabulated on dependency on p, q as shown in Table 1.

Actual Young’s modulus obtained by uniaxial tensile testing are approximately 70 - 73 GPa for aluminum and 110 - 130 GPa for copper. Referring these values,

Table 1. Relation between potential parameters (p, q) and calculated Young’s modulus.

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.

2.4. 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. ${\phi }_U$ is the function representing “unloading behavior” after plastic deformation proceeds, which is expressed by

${\phi }_U\left(r_{ij}\right)=-\frac{G}{p{\left(r_{ij}-\Delta r_0\right)}^p}+\frac{H}{q{\left(r_{ij}-\Delta r_0\right)}^q}+C_{\phi }$, (16)

where $\Delta r_0$ and $C_{\phi }$, which have units of length (m) and energy (J) respectively, are undetermined values until unloading starts. So, many different ${\phi }_U$ functions exist per unloading events. The functions ${\phi }_U$ connect each other as well as elastic potential ${\phi }_E$ by another function ${\phi }_P$, which prescribes a plastic route for between pairs of particles. The function ${\phi }_P$ used here is just a polynomial of interparticle separation $r_{ij}$. The function ${\phi }_P$ should be determined microscopically (from MD simulation) or empirical values obtained by experimental elastic-plastic testing.

Actually, only the elemental framework of simulating in plastic regime has been shown here. It is true that it requires further verification and discussion, but let us to omit those detailed description in the present paper at this time, to avoid complication of the context.

Figure 2. Interatomic potential functions representing elasticity and plasticity.

2.5. 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,

$F_D=-\alpha r_i$, (17)

and random force is given by

$F_R=F_R\left(\sigma \right)$, (18)

where α is a viscosity for atom and σ is the square root of variance. In MD scale, α and σ have the relation,

${\sigma }^{micro}={\left(\frac{2{\alpha }^{micro}{k}_{B}T}{\Delta t}\right)}^{\frac{1}{2}}$, (19)

and

${\alpha }^{micro}=m\frac{\pi }{6}{\omega }_D$, (20)

where ${\omega }_D=k_B\theta \bar{h}$ is Debye frequency ( $\theta$, $\bar{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 ${\eta }_{CG}$ for α and ${\zeta }_{CG}$ for σ, and resulted expressions are

${\sigma }^{macro}={\eta }_{CG}{\left(\frac{2{\alpha }^{macro}{k}_{B}T}{\Delta t}\right)}^{\frac{1}{2}}$, (21)

and

${\alpha }^{macro}={\zeta }_{CG}m\frac{\pi }{6}{\omega }_D$. (22)

Verifying these coarsening parameters, ${\eta }_{CG}$ and ${\zeta }_{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.

2.6. 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(a). The parameters required to arrangement are interparticle distance r₀ and lattice constant a₀. They have a certain relation, $r_0=a_0/\sqrt{2}$. The longitudinal length of the specimen $l_x$ is configured by $l_x=r_0+n_x\cdot a_0$, where $n_x$ is the number of division in x-direction and it also determine the length of the fcc unit cell. First, as shown Figure 3(b), a rectangular parallelepiped with

Figure 3. Configuration of macroscopic particles and modeling the specimen. (a) Configuration of macroscopic particles in the fcc lattice; (b)Sculpting from rectangular parallelepiped to resulted cylinder.

the dimensions $l_x$ and $l_y$ is arranged by particles, and then a cylinder with diameter of $l_x=D$ is sculpted from it. Once the diameter of cylinder $l_y$ is determined, division number by particles in diameter (width) direction is also determined by the same interparticle distance as the longitudinal direction. In loading simulations, the strain increase is prescribed by the constant velocity of moving particles which are located at two ends along longitudinal (z) direction. 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.

2.7. 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_{\phi }$ and $n_h$, 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.

$n_{\varphi }=2^i\left(i=1~4\right)$ (23)

$n_h=2^j\cdot n_{\varphi }\left(j=-3~3\right)$ (24)

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_{\phi }$ and $n_h$ are shown in Table 3. The arrangement of particles (view of external form and at cross-section) when the ratio of division parameters ( $n_h/n_{\varphi }$ ) is $2^2=4$ are shown in Figure 4.

Table 2. Calculation conditions for each models by the choice of the number of division in diameter (width) direction $n_{\varphi }$.

(a) (b)

Table 3. The number of macroscopic particles and actual height of the cylindrical specimens. (a) The number of macroscopic particles; (b) Actual height of the cylinder [mm].

Figure 4. The arrangement of macroscopic particles when the ratio of division parameters ( $n_h/n_{\varphi }$ ) is 2² = 4.

3. Results and Discussion

3.1. Elastic Analysis

3.1.1. Results of Tensile Simulations

Stress-strain curves obtained for $\left(n_h/n_{\varphi }\right)=2^0=1$ are shown in Figure 5. In every cases, they exhibits almost the same linear response because of elastic region. Young’s moduli E are estimated at strain 0.1% for each $n_{\varphi }$ s: E = 61.3, 58.6, 59.9, and 60.0 GPa for $n_{\varphi }=2,4,8,16$, respectively. These values are slightly lower than the theoretical value 68.7 GPa, which is derived directly from potential function used here. However, realistic elastic responses for the material are well reproduced by macroscopic particles. The relation between division number $n_{\varphi }$ and Young’s modulus (estimated up to 0.1%) are shown in Figure 6. When

Figure 5. Stress-strain diagrams (for the cases $\left(n_h/n_{\varphi }\right)=2^0=1$ ).

Figure 6. Relationship between Young’s modulus/interparticle strain and the number of grids.

$n_{\varphi }$ is 2¹ (=2), the largest modulus is obtained. As shown in Figure 6 at the same time, the dependency on division number $n_{\varphi }$ 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 $n_{\varphi }=2$, 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, $n_h$, also shows the same tendency of varying stiffness as explained for $n_{\varphi }$.

The dependence of Poisson’s ratio ν on the division ratio $n_h/n_{\varphi }$ is shown in Figure 7. When $n_{\varphi }$ 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 $n_h/n_{\varphi }=4$. The pictures are drawn at cross-sections at $x=0$, 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.

3.1.2. 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 $n_{\varphi },n_h$ and their ratio $n_h/n_{\varphi }$ 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.

Figure 7. Relationship between Poisson’s ratio and $n_h/n_{\varphi }$.

Figure 8. Stress distribution inside the cylindrical specimen (normal stress components inz direction; for $n_h/n_{\varphi }=2^2$ ; viewing on cross-section at $x=0$ ; in tensile testing).

Figure 9. Relationship between Young’s modulus and $n_h/n_{\varphi }$.

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, ${\epsilon }_{ij}$, to the global strain as for the specimen, ${\epsilon }_{\text{system}}$, is estimated approximately to be ${\epsilon }_{ij}/{\epsilon }_{\text{system}}=\left(1/2\right)\cdot \left(1/\sqrt{2}\right)=0.353$, where the ratio of lattice constant $a_0$ and interparticle distance $r_0$ is theoretically $r_0/a_0=1/\sqrt{2}$. As the number of division increases, the strain values between neighbor particles are approaching to 0.35, which agrees with theoretical consideration explained here.

(a)(b)

Figure 10. Relationship between longitudinal strain for neighbor particles and $n_h/n_{\varphi }$ (in compressive testing). (a) The dependency of strain on division parameters; (b) The fcc lattice and geometrical relationship between length parameters.

The stress distribution inside the specimen for the ratio between the numbers of division, $n_h/n_{\varphi }=4$, 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).

3.2. Plastic Analysis

3.2.1. 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 $n_{\varphi }=2^4=16$ and $n_h/n_{\varphi }=4$. The time increment $dT=4.58\times {10}^{-8}\text{s}$ is used for all these cases. The calculation conditions are summarized in 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

Table 4. Calculation conditions including plastic regime (tensile or compressive test).

Figure 11. Stress distribution inside the cylindrical specimen (normal stress components in z direction; for $n_h/n_{\varphi }=2^2$ ; viewing on cross-section at x = 0; in compression testing).

Figure 12. Stress-strain diagram including plastic regime.

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(b) visualizes the increase of stress value by strain increment of 0.5%. Tensile stress

Figure 13. Stress distribution in tensile loading (viewing the cross-section at x = 0). (a) Stress distribution at each tensile strain; (b) Changes of stress between strains.

inside the specimen increases almost uniformly, but it is always lower than those in surface particles.

3.2.2. 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. Yielding event is observed in the compression case in the same way as in tensile case. After yielding (approximately −1% of nominal strain), the magnitude of compressive stress drastically decreases when compressive strain reaches −3.5%, when sliding movement of neighboring particles is clearly observed. It is just like atomistic slip of fcc crystal structure, and those slip planes would correspond to (111) indicated by Miller’s index. Those slips move in oblique direction in the specimen and eventually reach both free surfaces. After slip events, the position of particles is being rearranged, during which some oscillation of stress is observed.

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 ${\epsilon }_z=-3.0\text{\%}$ to ${\epsilon }_z=-3.5\text{\%}$, there occurs slippage between neighbor particles in 45 degrees from compressive loading axis, which is not observed in tensile simulations. After slip behavior takes place, any fracture does not occur but stress relaxation behavior along multiple of slip planes is observed. The three-dimensional configuration of slip planes exactly corresponds to (111) plane of fcc structure.

Conceptual explanation of force acting between particles inside of the specimen is displayed in Figure 16, for (a) tensile and (b) compressive cases, respectively. In those pictures, a red line segment means tensile forces and a blue line

Figure 14. Stress-strain diagram.

Figure 15. Stress distribution in compressive loading (viewing the cross-section at x = 0). (a) Stress distribution at each compressive strain; (b) Changes of stress between strains.

Figure 16. Distribution of bonding forces (in the case of (a) tensile, or (b) compressive deformation). (a) Tensile loading case; (b) Compressive loading case.

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.

3.2.3. 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.

Figure 17. Stress-strain diagram (first tension, then compression). (a) Overall processes; (b) Magnified view around the second yield.

Figure 18. Stress-strain diagram (the case of first compression, then tension). (a) Overall processes; (b) Magnified view at the second yield.

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.

4. 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 framework 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.

Acknowledgements

The authors greatly acknowledge Nippon Steel Corporation (Kimitsu) for their financial support and for the support of research facilities.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Liu, W.K., Karpov, E.G. and Park, H.S. (2006) Nano Mechanics and Materials: Theory, Multiscale Methods and Applications. Wiley, Hoboken. s://doi.org/10.1002/0470034106
[2]	Ibrahimbegovic, A. (2010) Nonlinear Solid Mechanics: Theoretical Formulations and Finite Element Solution Methods (Solid Mechanics and Its Applications). Springer, Berlin.
[3]	Fish, J. and Belytschko, T. (2007) A First Course in Finite Elements. Wiley, Hoboken. s://doi.org/10.1002/9780470510858
[4]	Li, S. and Liu, W.K. (2007) Meshfree Particle Methods. Springer, Berlin.
[5]	Liu, G.R. (2003) Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific, Singapore. s://doi.org/10.1142/9789812564405
[6]	Tadmor, E.B., Ortiz, M. and Phillips, R. (1996) Quasicontinuum Analysis of Defects in Solids. Philosophical Magazine A, 73, 1529-1563. s://doi.org/10.1080/01418619608243000
[7]	Monaghan, J.J. (1992) Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics, 30, 543-574. s://doi.org/10.1146/annurev.aa.30.090192.002551
[8]	Koshizuka, S. (2018) Moving Particle Semi-Implicit Method. Academic Press, Cambridge.
[9]	Silling, S. (2000) Reformulation of Elasticity Theory for Discontinuities and Long-Range Forces. Journal of Mechanics and Physics of Solids, 48, 175-209. s://doi.org/10.1016/S0022-5096(99)00029-0
[10]	Liu, W.K. and McVeigh, C. (2007) Predictive Multiscale Theory for Design of Heterogeneous Materials. Computational Mechanics, 42, 147-170. s://doi.org/10.1007/s00466-007-0176-8
[11]	Saitoh, K. and Liu, W.K. (2009) Molecular Dynamics Study of Surface Effect on Martensitic Cubic-to-Tetragonal Transformation in Ni-Al Alloy. Computational Materials Science, 46, 531-544. s://doi.org/10.1016/j.commatsci.2009.04.025
[12]	Ushida, H., Ogata, S. and Kimizuka, H. (2011) Stability Analysis of Fcc Crystal under Local Shear Deformation by Molecular Dynamics. Journal of the Society of Materials Science, Japan, 60, 71-78. s://doi.org/10.2472/jsms.60.943
[13]	Saitoh, K. and Dan, K. (2012) The Method of Microscopic Strain Analysis Based on Evolution of Atomic Configuration for the Simulation of Nanostructured Materials. Journal of Society of Materials Science, Japan, 61, 162-168. s://doi.org/10.2472/jsms.61.162
[14]	Curtin, W. and Miller, R.E. (2003) Atomistic/Continuum Coupling in Computational Materials Science. Modelling and Simulation in Materials Science and Engineering, 11, R33-R68. s://doi.org/10.1088/0965-0393/11/3/201
[15]	Knap, J. and Ortiz, M. (2001) An Analysis of the Quasicontinuum Method. Journal of Mechanics and Physics of Solids, 49, 1899-1923. s://doi.org/10.1016/S0022-5096(01)00034-5
[16]	Qian, D., Wagner, G.J. and Liu, W.K. (2004) A Multiscale Projection Method for the Analysis of Carbon Nanotubes. Computational Methods in Applied Mechanical Engineering, 193, 1603-1632. s://doi.org/10.1016/j.cma.2003.12.016
[17]	Lu, H., Daphalapurkar, N.P., Wang, B., Roy, S. and Komanduri, R. (2006) Multiscale Simulation from Atomistic to Continuum. Coupling Molecular Dynamics (MD) with the Material Point Method (MPM). Philosophical Magazine, 86, 2971-2994. s://doi.org/10.1080/14786430600625578
[18]	Allen, M.P. and Tildesley, D.J. (1987) Computer Simulation of Liquids. Oxford University Press, Oxford.
[19]	Greenspan, D. (1991) Supercomputer Simulation of Liquid Drop Formation on a Solid Surface. International Journal of Numerical Methods in Fluids, 13, 895-906. s://doi.org/10.1002/fld.1650130707
[20]	Greenspan, D. (1997) Particle Modeling. Birkhauser, Basel. s://doi.org/10.1007/978-1-4612-1992-7
[21]	Wang, G., Al-Ostaz, A., Cheng, A.H.-D. and Mantena, P.R. (2009) Hybrid Lattice Particle Modeling: Theoretical Consideration for a 2D Elastic Spring Network for Dynamic Fracture Simulations. Computational Materials Science, 44, 1126-1134. s://doi.org/10.1016/j.commatsci.2008.07.032
[22]	Wang, G., Al-Ostaz, A., Cheng, A.H.-D. and Mantena, P.R. (2009) Hybrid Lattice Particle Modeling of Wave Propagation Induced Fracture of Solids. Computational Methods in Applied Mechanical Engineering, 199, 197-209. s://doi.org/10.1016/j.cma.2009.10.002
[23]	Lucy, L.B. (1977) A Numerical Approach to the Testing of the Fission Hypothesis. Astronomical Journal, 82, 1013-1024. s://doi.org/10.1086/112164
[24]	Libersky, L.D., Petschek, A.G., Carney, T.C. and Allahdadi, F.A. (1993) High Strain Lagrangian Hydrodynamics. Journal of Computational Physics, 109, 67-75. s://doi.org/10.1006/jcph.1993.1199
[25]	Hoover, W.G. and Posch, H.A. (1996) Numerical Heat Conductivity in Smooth Particle Applied Mechanics. Physical Review E, 54, 5142-5145. s://doi.org/10.1103/PhysRevE.54.5142
[26]	Randles, P.W. and Libersky, L.D. (1996) Smoothed Particle Hydrodynamics: Some Recent Improvements and Applications. Computational Methods in Applied Mechanical Engineering, 139, 375-408. s://doi.org/10.1016/S0045-7825(96)01090-0
[27]	Li, S. and Liu, W.K. (2002) Meshfree and Particle Methods and Applications. Applied Mechanics Review, 55, 1-34. s://doi.org/10.1115/1.1431547
[28]	Hoogerbrugge, P.J. and Koelman, J.M.V.A. (1992) Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhysics Letters, 19, 155-160. s://doi.org/10.1209/0295-5075/19/3/001
[29]	Satoh, A. (2011) Introduction to Practice of Molecular Simulation. Elsevier, Amsterdam. s://doi.org/10.1016/B978-0-12-385148-2.00003-3
[30]	Maruyama, S. and Choi, S.H. (2001) Molecular Dynamics of Heat Conduction through Carbon Nanotube. Thermal Science and Engineering, 9, 17-24.
[31]	Kinjo, T. and Hyodo, S.-A. (2007) Equation of Motion for Coarse-Grained Simulation Based on Microscopic Description. Physical Review E, 75, Article ID: 051109. s://doi.org/10.1103/PhysRevE.75.051109
[32]	Lei, H., Caswell, B. and Karniadakis, G.E. (2010) Direct Construction of Mesoscopic Models from Microscopic Simulations. Physical Review E, 81, Article ID: 026704. s://doi.org/10.1103/PhysRevE.81.026704