Unified nonlocal rational continuum models developed from discrete atomistic equations
In this paper, a unified nonlocal rational continuum enrichment technique is presented for improving the dispersive characteristics of some well known classical continuum equations on the basis of atomistic dispersion relations. This type of enrichment can be useful in a wide range of mechanical problems such as localization of strain and damage in many quasibrittle structures, size effects in microscale elastoplasticity, and multiscale modeling of materials. A novel technique of transforming a discrete differential expression into an exact equivalent rational continuum derivative form is developed considering the Taylor’s series transformation of the continuous field variables and traveling wave type of solutions for both the discrete and continuum field variables. An exact equivalent continuum rod representation of the 1D harmonic lattice with the non-nearest neighbor interactions is developed considering the lattice details. Using similar enrichment technique in the variational framework, other useful higher-order equations, namely nonlocal rational Mindlin-Herrmann rod and nonlocal rational Timoshenko beam equations, are developed to explore their nonlocal properties in general. Some analytical and numerical studies on the high frequency dynamic behavior of these novel nonlocal rational continuum models are presented with their comparison with the atomistic solutions for the respective physical systems. These enriched rational continuum equations have crucial use in studying high-frequency dynamics of many nano-electro-mechanical sensors and devices, dynamics of phononic metamaterials, and wave propagation in composite structures.
keywords:Unified nonlocal rational continuum models, Molecular dynamics, Continuum enrichment, Dispersion, Wave propagation
Experimental observations at micro/nanoscales have shown the inability of classical continuum mechanics in providing crucial insight into the basic phenomena of deformation, crack growth, and phase transformation in numerous material systems. Therefore, there is a need for convenient and generalized continuum models which can elucidate the detailed mechanisms of various micro/nano-electromechanical sensors and devices widely used in nano-technology, material science, and biology. Moreover, an efficient continuum equation can very accurately approximate the mean part of atomistic solutions and thus, can substantially improve the accuracy and economy of a multiscale analysis. However, a potent continuum equation should be able to describe certain mechanical phenomena which are due to the intrinsic characteristic of microstructural details of a material system. Unfortunately, as a coarse scale model, the conventional continuum theory fails to do so. In reality, however, no material is an ideal continuum but is discrete in nature with complicated spatially varying internal microstructure. Therefore, the conventional continuum description is no longer adequate to describe mechanical details of any material at these small scales. Thus, enrichment of conventional continuum theories is required to improve the agreement between the continuum predictions and the experimental observations.
Any existent microstructural heterogeneity of crystalline materials affects the dispersion of elastic waves, which is a real physical phenomenon that can be observed and studied experimentally. The shape of a general propagating wave changes over time due to the different propagating speeds of different harmonic components in it. This phenomenon is commonly known as dispersion. In the dynamics context, however, the most important motivation is to improve the dispersive characteristics of continuum equations as compared to the actual atomistic dispersion phenomena. In general, the dispersive characteristic of a material system is mainly dominated by the discrete characteristics of their crystal lattices and the range of interatomic interaction forces prevailing in the system (Brillouin, 2003; Kittel, 2008; Jirasek, 2004). However, these vital dispersion properties of short elastic waves in actual discrete atomistic systems are lost due to the standard homogenization procedure involved in the development of the continuum models (Jirasek, 2004).
Over the last century, many researchers have expended their valuable efforts in constructing several enrichments of the standard continuum theory to capture the important features of the microstructure of any heterogeneous material system. All these variety of enrichment models are classified under generalized nonlocal continuum models. Most of these enrichment techniques are inspired by the seminal idea of the Cosserat brothers (Cosserat and Cosserat, 1909). All these generalized Cosserat theories consider additional fields which are independent of displacement field and provide supplementary information on the small scale kinematics. One such major class of enrichments is based on the idea that the total elastic energy density is a function of higher order gradients (, , , etc.) as well as the first gradient of the displacement fields () (Toupin, 1962, 1964; Mindlin and Tiersten, 1962; Mindlin, 1964; Mindlin and Eshel, 1968; Koiter, 1964, 1969). This class of enrichment is known as strain-gradient theory of elasticity in the field of nonlocal continuum mechanics. The simplest enrichment of this class assumes the strain energy density as the function of and . This results in a wave equation, which has an extra fourth order partial derivative term as compared to the classical wave equation. The dispersion relation of this enriched wave equation, for the positive values of higher-order elastic modulus, becomes erroneously opposite to that of the discrete atomistic model and the experimental observations at higher frequency (Jirasek, 2004). Therefore, this strain-gradient elasticity model can give a reasonable approximation of the actual dispersion characteristics only for some negative values of higher-order elastic modulus. However, the assumption of negativeness of the higher-order elastic modulus gives physically absurd results like the loss of convexity of elastic potential and unstable behavior of short elastic waves. Therefore, the enriched wave equation obtained by this model is termed as the “bad Boussinesq problem”.
Another class of nonlocal enrichments is known as models with mixed spatial-temporal derivatives which ultimately tries to overcome the snags of the “bad Boussinesq problem”. Fish et al. (2002b) proposed a method in which they replace the fourth order partial derivative term by a term with mixed derivative, . This model gives reasonable accuracy in dispersion prediction at moderate frequency and overcomes the problem of instability for very short waves. This is known as the “good Boussinesq problem” and its extension to multiple dimensions is complicated (Fish et al., 2002a). Metrikine and Askes (Metrikine and Askes, 2002; Askes and Metrikine, 2002, 2005) proposed a different continualization method of this class which includes an additional heuristic model parameter. With a proper choice of parameter, this model can reasonably approximate dispersion and can avoid instability for short elastic waves. However, this model has a fourth order derivative term and higher-order boundary conditions are required to solve the problem.
The popular integral type nonlocal elasticity is based on weighted spatial averaging. This is also an elegant generalization of Cosserat theories (Cosserat and Cosserat, 1909). Eringen and Edelen (Edelen et al., 1971; Eringen, 1965, 1972, 1983) have, mainly, advanced the theories of integral type nonlocal elasticity. In this type of enrichments, the constitutive relation at a point of continuum involves weighted averages of the state variables over a certain neighborhood of that point (Bazant and Jirásek, 2002). In their constitutive relation, they use a dimensionless even function called nonlocal weight function which determines the accuracy of the model. Eringen (1983) derived the approximate expression of this nonlocal weight function by studying the relationship between 1D atomic lattice with nearest-neighbor interaction and the nonlocal integral models. This has inspired a plethora of research articles on the extension of this model for higher dimensions by several other researchers (Peddieson et al., 2003; Wang, 2005; Wang et al., 2006; Reddy, 2007; Reddy and Pang, 2008; Aydogdu, 2009; Reddy, 2010; Narendar and Gopalakrishnan, 2009, 2010a, 2010b, 2010c; Narendar et al., 2011). This model has good reasonable accuracy in predicting the dispersion relation in 1D. However, this too is a reasonably approximate method, but not an exact method. Eringen’s integral model would reproduce the dispersion relation for D atomistic model only for wavelengths larger than , where is the lattice spacing (Jirasek, 2004). The expression of weight function, proposed by Eringen, is very complicated and not very practical to use (Jirasek, 2004).
Observing the difficulties in determining microstructure-dependent length scale parameters (Lam et al., 2003; Maranganti and Sharma, 2007), Park and Gao (2006, 2008) used the modified couple stress theory proposed by Yang et al. (2002) to develop nonlocal Euler-Bernoulli beam which has only one material length scale parameter. Ma et al. (2008) developed a microstructure-dependent Timoshenko beam model based on the same modified couple stress theory. Afterwards, Reddy (2011) developed microstructure-dependent nonlinear Euler-Bernoulli and Timoshenko beam theories for functionally graded beams using the principle of virtual displacements, which has motivated several applications.
In the current paper, we propose a novel and elegant technique of transforming a discrete differential expression into its exactly equivalent rational continuum derivative form. Using the concept of Taylor’s series expansion of sufficiently smooth displacement field we show that, for a traveling wave type of solution, a discrete mechanical model can exactly be realized as a partial derivative term premultiplied by a linear differential operator with constant coefficients which incorporates the discrete microstructural information. Using this technique, we have derived the exact equivalent wave equation of a 1D harmonic lattice with nearest neighbor interaction as well as with non-nearest neighbor interactions. Most interestingly, the gradient type nonlocal elasticity models (Toupin, 1962, 1964; Mindlin and Tiersten, 1962; Mindlin, 1964; Mindlin and Eshel, 1968) and the integral type nonlocal elasticity models (Edelen et al., 1971; Eringen, 1965, 1972, 1983) are shown to be the special cases of this generalized nonlocal rational continuum formulation. Development of some other higher order continuum equations namely nonlocal rational Mindlin-Herrmann rod and nonlocal rational Timoshenko beam equations are shown to be straightforward using this novel idea. These nonlocal rational continuum models have only one microstructure-dependent length scale parameter. However, these equations show substantially improved agreement in the dispersion characteristics with the linearized atomistic equations. Most importantly, these equations are simple and can be solved analytically and/or semi-analytically in the Laplace transform based spectral finite element method (NLSFEM) (Murthy et al., 2011; Patra et al., 2014) framework with ease. Moreover, simple approximations in the nonlocal rational continuum formulations generate some well known gradient type and/or integral type nonlocal elasticity equations which can be solved analytically and/or using the conventional finite element method.
The rest of this paper is organized as follows. The novel nonlocal enrichment of continuum equations on the basis of atomistic dispersion relation is presented in Section . In Section , an exact equivalent nonlocal rational rod realization of 1D harmonic lattice with nearest neighbor interaction as well as with non-nearest neighbor interactions and the derivation of several discrete dispersion differential operators are described in detail in Section . The variational formalism for the general derivation of nonlocal rational continuum equations with appropriate boundary conditions is presented in Section . Section describes the details of nonlocal rational Mindlin-Herrmann rod approximation of a 2D harmonic lattice. The details of nonlocal rational Timoshenko beam approximation of a 2D harmonic lattice system is presented in Section . Section presents the analytical and numerical studies. Analytical solutions for dispersive characteristics of a 1D harmonic lattice with nearest neighbor interaction as well as with next-nearest neighbor interactions are presented in Section . The comparison of numerical solutions by both the molecular dynamics (MD) and NLSFEM is presented in this section. The dispersive characteristics of nonlocal rational Mindlin-Herrmann rod and nonlocal rational Timoshenko beam are compared with their respective atomistic predictions in Section and Section , respectively. The small scale effects on the natural frequencies of a simply supported nonlocal rational Timoshenko beam are studied in Section . The paper concludes with a summary in Section .
2 Unified nonlocal rational continuum enrichment on the basis of atomistic dispersion relation
The dispersion prediction of elastic waves deteriorates when an actual discrete system is modeled as a homogeneous continuum. Here we describe a novel technique to incorporate the discrete dispersion characteristics of an actual atomistic system in its equivalent continuum equations. In this section, the procedure of finding proper discrete dispersion differential operators and the variational derivation technique of several unified nonlocal rational continuum equations are described.
2.1 An equivalent nonlocal rational nod for 1D harmonic lattice: Exact discrete dispersion relation
The enriched 1D wave equation is developed directly from the atomistic equation of 1D lattice. For the continuum enrichment, an equivalent 1D harmonic lattice model (Fig. 1) of any 1D lattice with nonlinear potential is considered. The lattice consists of identical atoms of equal mass connected by identical, massless, lossless, perfectly linear springs of stiffness . The initial position and deformed position of the atom are given by
where, is the lattice constant and is the out-of-equilibrium displacement. Force exerted on the atom by the atom in the harmonic lattice is proportional to the difference . First, the formalism is developed considering nearest-neighbor interactions only. Afterwards, the formalism is extended for non-nearest-neighbor interactions. The homogeneous equation of motion of the atom () in the harmonic lattice obeys Newton’s second law as
Considering the traveling wave solution of the type , the dispersion relation and the expression of group velocity , and phase velocity are given by (Kittel, 2008)
where and are the longitudinal wavenumber and the circular frequency, respectively and denotes absolute values. Now assuming the 1D chain of atoms as a continuum rod of uniform cross-section with Young’s modulus and density , letting and , traveling wave solution will be of the type . Here, the discrete solution is a subset of the continuum solution . Considering this, one can rewrite equation (2) as
Or, the homogenized Eq. (5) can be expressed as
where denotes the right-hand side external forcing and L denotes the left-hand side homogeneous differential operator in Eq. (5). The dispersion relation and the expression of group velocity , and phase velocity associated with Eq. (5) are given as
The deterioration in the above expressions of dispersion relation and group velocity (Eq. (2.1)) from its original atomistic one (Eq. (2.1)) is mainly due to the limiting consideration of . However, the atomic separation , which is of the order of m in reality, is many order bigger unit of length as compared to the smallest possible Plank’s length m considered in physics. Therefore, the consideration of is wrong for a realistic atomistic systems. The limiting equation (5) gives erroneous results for actual 1D harmonic lattices other than the long wave condition.
To overcome this lacuna, the exact equivalent continuum equation is derived in a different manner. Considering the Taylor’s series expansion of the sufficiently smooth function and for , one can write the discrete differential portion of the right-hand side of Eq. (2) as
Here, is an infinite sum and cannot be neglected when . Therefore, the participation of is taken into account by premultiplying a discrete dispersion operator with the partial derivative and the Eq. (8) is rewritten as
Here the subscript of denotes the order of the partial derivatives of the field variables with which it is pre-multiplied. Considering the continuum traveling wave solution , Eq. (9) can be reduced to the form
Now manipulating the Eq. (10) as in the following manner, the modified equation is rewritten as
The Eq. (2.1) can be reduced as
From which, it follows that
The linear differential operator contains infinite terms with constant coefficients. It is seen from Eqs. (8) - (12) that the discrete differential expression in Eq. (8) can be transformed into an exact derivative term with a multiplier linear differential operator with constant coefficients without any approximation. This novel transformation is true for any non-zero value of .
The solution below the resolution of atomistic solution is physically meaningless. Therefore, considering as an equivalent representation of and using the development of Eqs. (8) - (13), Eq. (2) is rewritten as
From which, it follows that
Or, the nonlocal rational continuum rod Eq. (15) can be expressed as
where denotes the corresponding nonlocal rational differential operator in Eq. (16) and for no external forcing. This new enriched wave equation (15) has exactly identical dispersion relation as the atomistic Eq. (2.1) and can exactly represent its atomistic counterpart for all wavenumbers/frequencies. The novel nonlocal rational rod model contains one additional microstructural length scale parameter . The enriched continuum equation (14) reduces to it’s classical form at the long wave limit i.e., when . This equation is very simple to solve in the frequency domain and can be solved semi-analytically using NLSFEM to give very accurate predictions about discrete atomistic systems.
For small frequency/long wave length condition (i.e., when is small), neglecting higher degree terms above the term containing , Eq. (12) can be approximated as
This approximate Eq. (18) resembles the wave equation of strain-gradient elasticity also known as “bad Boussinesq problem” (Jirasek, 2004) but with the corrected sign. This shows that the corrected wave equation of strain-gradient elasticity is a small wavenumber approximation of the novel nonlocal rational continuum model.
Most of the brittle materials have short range interactions. Therefore, an atomistic model with nearest-neighbor or next-nearest-neighbor interactions can very satisfactorily describe the system dynamics of brittle materials. In case of metals, the range of interactions observed are of quite long range and have been found between as many as nearest neighbors (Kittel, 2008). For non-nearest-neighbor interactions, the generalized atomistic dispersion relation can be written as (Kittel, 2008; Jirasek, 2004)
where is the equivalent harmonic spring constant of nearest neighbor interaction. For any nonlinear interatomic potential as considered in Rafii-Tabar and Mansoori (2004), magnitude of decreases with the increasing values of . For example, in case of a Lennard-Jones type potential, the magnitude of is approximately order less with respect to the magnitude of . Therefore, by straightforward extension, the generalized enriched continuum wave equation for a harmonic lattice with non-nearest-neighbor interactions can be written as
Or, in simplified form, the enriched wave equation is obtained as
where is the stiffness ratio and is the second order discrete dispersion operator for the nearest neighbor interaction. Eq. (21) is the exact equivalent continuum representation of the 1D harmonic lattice with non-nearest neighbor interactions. For the sake of clarity, we rewrite Eq. (20) for the next-nearest neighbor interaction (i.e., for ) in detail as
Expanding the term, one can simplify the Eq. (22) as
After similar manipulation as in Eq. (2.1), the exact continuum form of the discrete 1D harmonic lattice with next-nearest neighbor interaction is obtained as
From which, it follows that
Eq. (25) clearly shows that the interactions between the non-nearest neighbor atoms are taken into account exactly by the inclusion of higher order derivative terms in the equation. It can be seen that for the limiting condition of and for nearest neighbor interaction, Eq. (25) reduces to the classical rod equation (5). However, the long wave consideration (i.e., ) reduces the Eq. (24) to a correct form of strain-gradient elasticity equation which resolves the sign paradox of the “bad Boussinesq problem”. Following similar approach as above, exact equivalent continuum wave equations for distant neighbor interactions can be obtained from Eq. (21).
Eq. (25), in fact, includes a dependance on the immediate as well as non-nearest neighborhood of the point under consideration. In broad sense, continuum models involving a characteristic length scale parameter or a weighted summation/intregal are generally considered as nonlocal continuum models (Bazant and Jirásek, 2002). These enriched wave equations contain the linear differential operator of infinite terms with constant coefficients, which are functions of the lattice parameter . Therefore, these type of enriched continuum formulations (i.e., Eqs. (14), (25)) can be considered as generalized nonlocal continuum models as they contain a characteristic length parameter (lattice constant) and the weighted nonlocal summations. The nonlocal expression proposed by Eringen (1983) can be shown as an approximation of Eq. (15).
2.1.1 Discrete dispersion differential operators
In the similar fashion, the expressions of some important discrete dispersion differential operators s associated with partial derivatives of different orders are derived. The field variable is a very smooth function. Therefore, one can consider the Taylor’s series expansion of about any atom locations in the lattice and write
For the same traveling wave solution of the type , Eq. (26) reduces to the form
Manipulating the Eq. (27) in the following manner, one can obtain
Therefore, the first order discrete dispersion differential parameter is obtained as
Following the similar approach as above, one can show that
From the above details, it is reasonable to conclude that
These s can be used to derive some other useful unified nonlocal rational continuum equations very efficiently.
2.1.2 Variational formalism: energy functional for general formulation
The well-posedness of the problems, which involves the requirement of appropriate boundary conditions, can be obtained adopting variational formalism. Therefore, the governing equations with appropriate boundary conditions are to be derived from an appropriate elastic potential functional using Hamilton’s principle. The Hamilton’s principle for deformable conservative systems can be stated as
where is the kinetic energy, is the elastic potential energy, and is the work done by external loads on the system. Here , are the initial and final times, respectively. Observing the characteristic of the Eqs. (15) and (25), the elastic potential energy functional for a nonlocal linear elastic system with small-strain condition can be assumed as
where and are the generalized nonlocal stiffness parameters which are functions of discrete dispersion differential operators and the local strains of the order , respectively. The function is a function of , which takes account of the microstructural interaction of the continuum. Here, the expression of depends upon the order of gradients of the independent displacement field variables in the strain field and the degree of that order partial derivative term in the quantity . The nonlocal rational stress tensor , associated with the order strains, can be defined as
The local strains are related to the local Cauchy stress tensor by the relation
Considering the criteria , the exact expression in the Eq. (36) can be approximated as
where is another nonlocal parameter other than the length scale parameter . From which, it can be shown that the nonlocal constitutive equation for higher dimensions will be as
The above derivations clearly show that the popular strain gradient elasticity models Toupin (1962); Mindlin (1964) and integral nonlocal elasticity models (Eringen, 1983) are special approximations of the novel unified nonlocal rational continuum model.
The next two sections describe the development procedures of some higher-order nonlocal rational continuum equations with improved dispersive characteristics and/or nonlocal properties in general.
2.2 A nonlocal rational higher-order rod approximation of 2D harmonic lattices
Any D atomistic system can actually be idealized as D higher-order rod and/or higher-order beam continuum with a very high accuracy for static and low frequency dynamic problems. However, these homogenized classical continuum models in the form of partial differential equations fail to represent the discrete atomistic systems satisfactorily in case of high frequency dynamic conditions. Therefore, with proper enrichment of these local continuum equations, the discrepancies between the continuum predictions and experimental observations can be resolved satisfactorily.
A simple D triangular lattice system with hexagonal interaction (Fig. 2) is considered for modeling wave propagation in D atomistic system. The triangular lattice system is chosen because it represents as either the basal plane of many hexagonal close-packed (hcp) lattice systems or the () plane of many face-centered cubic (fcc) lattice systems (Kittel, 2008). There are an enormous number of other materials (not only Graphene, Graphyne, Germanene, palladium, Molybdenum disulfide, etc.) that are 2D materials and are very interesting from the basic science and applications point of view. For many other 2D lattice systems, the formalism can be extended easily. Therefore, at first, the atomistic equations and dispersion relations for the D system is described. Then, higher enriched order rod and beam approximations for 2D triangular lattice system are obtained simultaneously using atomistic information. The atom in column and row location in the D lattice system, where if is odd and if is even, is indexed by . The position of the atom in the - plane at time is given by (Fig. 2)
Here , (Fig. 2), where is the equilibrium atomic separation and are out-of-equilibrium atomic displacements along the respective directions in the - plane. Considering nearest-neighbor interactions (Fig. 2) at first, the Hamiltonian for the system can be written as (Friesecke and Matthies, 2003)
Here denotes the Euclidian norm on , denotes the momentum vector of the atom, and are the lattice basis vectors and , respectively. denotes the interatomic potential function of any arbitrary type. However, in the present work, we consider potential incorporating stretching mode only. Although, the nearest-neighbor interaction is assumed initially, however, this formalism can be easily extended for non-nearest neighbor interactions. The potential may be any arbitrarily differentiable function, but for continuum approximations, it is reasonable to consider equivalent harmonic potential of the form
where is the equivalent harmonic spring constant and denotes the distance of separation between any two nearest neighbor atoms. For any other nonlinear potentials, the constants of equivalent harmonic springs are obtained from the curvature of the nonlinear potential. The equations of motion of the atom can be written as
Expanding Eq. (43), the equations of motion can be written as
Where the forcing terms are defined as . Here, the Hamiltonian equations of motion are still nonlinear, due to the frame-indifference of the interatomic forces (Friesecke and Matthies, 2003). This is analogous to the geometric nonlinearity in continuum elasticity. Therefore, Eq. (2.2) must be linearized for studying the phonon dispersion characteristic of the 2D lattice. The linearized equations can be obtained from Eq. (2.2) as
Where , , and are matrices as given below
Consider the D atomistic system as a deep axial waveguide under loading in X-direction only (Fig. 2). Let the wave propagate in the X-direction as considered in Friesecke and Matthies (2003). Motivated by the atomistic solutions of Eq. (2.2), the approximate continuum deformation field variables are assumed to be of the form
where and are displacement field variables corresponding the neutral axis along X-direction and takes account of the Poisson’s effect. The local strain fields corresponding to the field variables in Eq. (47) are computed as
Here, , , and are the second moment of area of cross-section, the mass density, and the area of cross-section of the 2D atomistic lamina, respectively. In Eq. (2.2), and are adjustable parameters (Doyle, 1997; Gopalakrishnan, 2000) and , , and are stiffness parameters, respectively. For generality, the continuum parameters , can be obtained using the following relation (Kittel, 2008)
where is potential energy per atom in a unit cell under uniform tensile deformation, is the effective surface area occupied by the atom in the X-Y plane, is the thickness of the D atomistic lamina and , are strains, respectively. From Eq. (50), the parameters for the triangular lattice are obtained as , . The parameter is obtained from and . The associated boundary conditions for Eq. (2.2) are specified as
The equations of motion in (2.2) are the so called local Mindlin-Herrmann rod (LMHR) equations (Mindlin and Herrmann, 1951; Doyle, 1997). This classical Mindlin-Herrmann rod (LMHR) model (Eq. (2.2)) shows poor dispersive characteristics as compared to its atomistic counterpart (Eq. (2.2)), because of the homogenization procedure involved in the formulation of continuum equations. To overcome this snag, the homogenized L operator in Eq. (2.2) is to be replaced with the nonlocal operator containing appropriate s. However, the nonlocal elastic potential functional form corresponding to the unified nonlocal rational Mindlin-Herrmann rod (NRMHR) equations can be written according to Eq. (33) as