Estimation of Isotropic Hyperelasticity Constitutive Models to Approximate the Atomistic Simulation Data for Aluminium and Tungsten Monocrystals

1 Introduction

Hyperelastic models in finite deformation formulation are most often used for rubber materials and in biomechanics[Ogden(1984);Jemioło(2002);Pandol fiand Vasta(2012)],although also for porous elastomers [Danielsson,Parks,and Boyce(2004)]and textile composites[Kueh(2012)]or even for auxetic materials[Ciambella and Saccomandi(2014)].Recently due to the interest in graphene and its mechanical properties,hyperelastic models have been applied to carbon nanotubes[Ling and Atluri(2007)]and single-layer graphene[Xu,Paci, Oswald,and Belytschko(2012)].The hyperelastic multiscale modelling technique has been used to predict the elastic properties of polymers,but the atomistic force fields usedin[Valavala,Clancy,Odegard,and Gates(2007)]are not very reliable.Based on molecular dynamics simulations of a metallic glass,a modified isotropic hyperelastic Hencky model was proposed in[Henann and Anand(2011)].For crystals and especially at the nanoscale,the linear theory dominates[Tadmor and Miller(2011)]although elastic deformations can be significant [Clayton(2011)],e.g.there is experimental evidence that the single-crystal quartz behaves elastically up to a pressure load of 10 GPa and 16%volumetric deformation[Wang,Mao,Jiang,and Duffy(2014)].Most crystals generally behave anisotropically.The proposal of anisotropic hyperelasticity on the basis of crystallographic motivated structural tensors can be found in[Schroeder,Neff,and Ebbing(2008)],but this proposal is not consistent with Hooke’s relation for small deformations,making it difficult to justify physically.A model for a subclass of cubic anisotropy material with limited anisotropy ratios was proposed by[Kambouchev,Fernandez,and Radovitzky(2007)],but there was no fitting of the parameters.An interesting model for describing cubic symmetry material with a fitting to abinitio data was presented in[Salvetti(2010)].Unfortunately,this model needs to determine 35 parameters and is not polyconvex.The lack of mathematical rigour and the results obtained do not justify the complexity of the model.The concept of nonlinear elasticity based on high order elastic constants quite popular in physics is described,for example,in[Clayton(2011)]and utilised for III-N compounds in[łepkowski,Majewski,and Jurczak(2005)].

It is known that symmetry considerations play an important role in the study of physical phenomena.From symmetry principle:if a crystal is invariant with respect to certain symmetry elements,any of its physical properties must also be invariant with respect to the same symmetry elements and Curie laws,it arises that the symmetries in the physical properties of a material cannot be smaller than the symmetry of the crystal,but can be bigger,so a cubic crystal can be elastically isotropic or nearly isotropic,see the International Tables for Crystallography[Paufler(2003)].Tungsten and aluminium although they are cubic crystals exhibit an elastic behaviour close to isotropy.It must also be remembered that classification of linear-elastic materials due to the symmetries leads only to eight classes of symmetry,while in crystallography are distinguished 32 crystallographic point groups and 7 Curie groups[Zheng(1994);Dimitrienko(2002);Bóna,Bucataru,and Slawinski(2007)].Distribution of linear-elastic materials by symmetry has nothing to do with crystallography.It results from the properties of fourth order Euclidean symmetric tensors(from linearity of Hooke’s relation and properties of three-dimensional Euclidean space)[Ostrowska-Maciejewska and Kowalczyk-Gajewska(2013)].

Molecular and ab simulations,in contrast to experiments,naturally allow for the direct determination of the internal energy function,from which we can determine the elastic energy,but then there are many different ways to calculate atomic-level stress tensor,with probably the most commonly used virial stress[Zhou(2003)]. Unfortunately,none of the atomic-level stresses are not the Cauchy or any other mechanical stresses[Shen and Atluri(2004)].Using this convenience of possessing the elastic energy and additionally taking into account a further critique of the atomistic definitions of the stress tensor[Murdoch(2007)],the sought for stored energy function is fitted directly to the atomistic elastic energy from simulations.

The present article first briefly discusses some physical,mathematical and numerical requirements which should,in our opinion,be satisfied by the constitutive models.Then we give the results of molecular numerical tests for aluminium and tungsten,being data for the parametrisation and isotropic average estimation of the obtained elasticity constants.Then,the main part of the article discusses the fitting of a hyperelastic model to the molecular data.Some conclusions are drawn in the final section.

2 Hyperelasticity constitutive model

2.1 Fundamental relations and constitutive model formulation

From the laws of the conservation of energy,linear momentum,and rotational momentum,together with the mass conservation rule,the following relations for the time derivative of the stored energy function(SEF)W can be obtained[Ogden(1984);Jemioło(2002)]:

where ρ0and ρ are density fields in,respectively,the reference and the current con figurations.Since W is purely elastic in origin,it is sometimes referred to as the elastic potential energy function[Ogden(1984)].The symmetrical stress tensorTappearing in Eq.(1)is the so called second Piola-Kirchhoff stress tensor,andSis the first Piola-Kirchhoff stress tensor(the nominal stress tensor),see[Bonet and Wood(2000)].However,using classical notion of work-conjugacy[Hill(1968)]many different work-conjugate stress-strain pairs can be set[Atluri(1984);Dłuzewski(2000)].It is worth recalling that the following identities between the stress tensors(σ is Cauchy’s stress tensor)hold:

The capital letter T as an upper index stands for tensor transposition.The tensorEin Eq.(1)is the Lagrange strain tensor,and can be expressed via the right Cauchy-Green stretch tensor

whereC=FTFandIis the identity tensor.From the local formulation of the energy conservation law,for sufficiently regular SEF W(F,X),there results the constitutive equation

In Eq.(4),Fis the deformation gradient,andXis the vector describing the location of the particle in the chosen coordinate system,see[Ogden(1984);Bonet and Wood(2000)].

The constitutive relation for the SEF in the form Eq.(5)can,in view of Eq.(1),in the Lagrange description be expressed in differential form

By further differentiation,we also get a fourth order tensor D,called the elasticity tensor in the material,in the Lagrange description:

ForE=0or,equivalently,the elasticity tensorcorresponds to Hooke’s tensor.The function Ŵ can depend onCthrough its invariants,for example:

where tr denotes the trace,cof the cofactor,and det the determinant of a tensor.In this paper,the following general family of SEF functions is considered[Ciarlet(1988)]:

in which Γ(J)is a convex function dependent only on volumetric changes and the ai,bjare material parameters.The function above,with appropriate restrictions on the parameters,automatically guarantees the fulfilment of the polyconvexity conditions[Ciarlet(1988)].There have been many proposals for Γ(J)functions for weak compressibility[Hartmann and Neff(2003)],but much less for full compressibility,here three of them are used[Jemioło(2002)],see Table 1.

Table 1:Γ(J)functions applied to the analysed problems[Jemioło(2002)].

2.2 Conditions for hyperelastic constitutive models

A detailed analysis of the physical,mathematical,and numerical requirements for the SEF function can be found,for example,in[Jemioło(2002);Ciarlet(1988);Silhavy(1997)],which discuss mathematical issues concerning the requirements for the function:convexity,polyconvexity,quasi-convexity,and rank-one convexity guaranteeing the existence of solutions to a hyperelasticity boundary value problem.Physically,these requirements translate into:the real value of the speed of the waves propagating in the elastic body,and SEF behaviour under extreme deformations.A standard requirement for the SEF is also the fulfillment of the condition of the energy-and stress-free natural state.A reasonable requirement from the standpoint of engineering is the consistency of the SEF with Hooke’s relation for small strains and small deformations,and according to Eq.(7)this is equivalent to consistency with the Saint Venant-Kirchhoff material model in the neighbourhood of the natural state[Jemioło(2002);Bonet and Wood(2000);Ciarlet(1988);Dłużewski(2000)].A very important aspect of the usefulness of a particular SEF is the ease of determining its parameters,e.g.,for Ogden material[Ogden(1984)]using the eigenvalues of the elongation tensor the parameter fitting procedure is a highly nonlinear process demanding very advanced optimization methods and with absolutely no guarantee of obtaining meaningful and unambiguous results[Ogden,Saccomandi,and Sgura(2004)].In the class of a general family of SEF function used herein,the fitting procedure for Eq.(9)is relatively straightforward and can be carried out with standard tools[Wolfram Research(2010)].Successively adding monomials improving the fitting is easy[Hartmann(2001);Steinmann,Hossain,and Possart(2012)].

For the present article,it has been decided to use fully compressible models.A common assumption appearing in the literature is that the corresponding strainenergy function can be decomposed additively into volumetric and isochoric parts.Even intuitively it easy to show that this decomposition is not physically realistic for anisotropic materials,e.g.,with this decomposition,any anisotropic cube under hydrostatic tension deforms into another cube instead of a hexahedron with non-parallel faces[Sansour(2008);NÃ Annaidh,Destrade,Gilchrist,and Murphy(2013);Vergori,Destrade,McGarry,and Ogden(2013)].The fallacy and the consequences of the assumption concerning the independence of the volumetric and isochoric(or deviatoric)behaviour of the material even for isotropy is unfortunately not well known[Henann and Anand(2011);Peng and Landel(1975)].

2.3 The simplest hyperelasticity constitutive models

Among the simplest models of isotropic hyperelasticity two deserve special attention.The first of these models is the Saint Venant-Kirchhoff(SVK)model.The stored energy function in the case of the SVK material model is defined by the function

where the constants λ0and µ0are identical as in the small deformation theory(Lamé’s constants)and depend on technical constants through the relations

The fact that the only parameters in the model are the same as for the linear elasticity model makes it unique.This model for small deformations is consistent with Hooke’s linear relation.It is possible then to formulate a condition that all the further considered hyperelasticity constitutive models of isotropic materials for small deformations should be consistent with the SVK model.This means that in each of the considered models,some combinations of the material parameters can be expressed by the parameters of a linear-elastic material.Determining the relation between the material parameters as in Eq.(9)and λ0,µ0can be done in different ways.In this paper,in order to determine the relation between the parameters of the selected family models,Eq.(9),the SEF functions were determined for the following tests:uniaxial strain,simple shear,volume deformation,and biaxial strain.The deformation gradient tensor representations in Cartesian coordinates for the atomistic tests are defined in Section 3.2.1-3.2.4.Then,the obtained SEF forms for each test have been expanded into a power series(in the case of the shear test,around zero;in other cases,around one).A comparison of these expansions with the expansions of the SEF for the SVK model allows interpreting the parameters λ0,µ0in the considered cases.The obtained relations are not written in an explicit form,due to its relatively complex form.However,for all analysed SEFs,the following algorithm for determining the material parameters was applied:

1.Provide a transition of the given model to the SVK model by evaluating in the above described manner the interpretation of the parameters λ0and µ0,which have the same meaning as the parameters in Hooke’s relation for small deformations;

2.Assumption of the λ0and µ0parameters on the basis of evaluations made as for the small deformation linear theory using atomistic simulations(compare the considerations presented in Section 4);

3.Determination of the rest of the parameters based on a nonlinear regression procedure with restrictions resulting among others from Section 2.2.

It should also be noted that the statistical evaluation in the case of all analysed constitutive models is performed under the assumption that the parameters λ0andµ0are constant and taken on the basis of the calculations made in Section 4.This means that to obtain a better fit of the model to the experimental data would probably be possible even in the case of the SVK model(the use of the optimization procedure to determine the parameters λ0and µ0),but achieved at the cost of compatibility of the model with the expectations of the classical Hooke’s relation.

There is another hyperelasticity model worth noting,in which there are only two constants,with interpretations identical to the case of the small deformations theory.This model usually is called a logarithmic(LN)model,for which the SEF is a function of the logarithms of the principal stretches,i.e.,the eigenvalues of the deformation tensorF:

where J= λ1λ2λ3,see(Jemioło,2002).SEF can also be written in a different version,i.e.as a function of the invariants of the lnUtensor:

There is a clear resemblance of the SEF in Eq.(13)with the SEF function in the case of a linear isotropic elastic material for the theory of small displacements(it is sufficient to replace the small strain tensor with the lnUtensor).The logarithmic model fulfils the Legendre-Hadamard condition in the moderately large deformation range for all positive material parameters,Hill’s inequalities with respect to Hencky strain,and the Baker-Ericksen inequalities in the whole deformation range[Bruhns(2001)],which makes it superior to other physically linear hyperelasticity models(e.g.,the Saint Venant-Kirchhoff),see[Jemioło(2002);Böhlke and Bertram(2002)].Unfortunately,the use of logarithmic models has a serious drawback:it is much more complicated in numerical implementation and requires a special approach[Itskov(2009)].

3 Atomistic modelling of materials

All molecular simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator(LAMMPS)[Plimpton(1995)].Due to interest in the static behaviour of the material,the Molecular Statics(MS)approach[Tadmor and Miller(2011);Maździarz,Young,and Jurczak(2011);Maździarz,Young,Dłużewski,Wejrzanowski,and Kurzydłowski(2010)]was used,which is more appropriate in such cases from molecular dynamics.The simulation region for all tests was chosen as 4x4x4 unit cells to properly enforce periodic boundary conditions(PBC)in all directions,fulfilling the requirement of the minimum-image convention and avoiding “self seeing”of the atoms[Tadmor and Miller(2011)].

3.1 Embedded atom model for Aluminium and Tungsten

Aluminium and Tungsten are metals arranged in a face-centered cubic(FCC)and body-centered cubic(BCC)structure,respectively.The potentials based on the Embedded Atom Model(EAM)have been successful in the description of metallic systems[Tadmor and Miller(2011)].

For FCC Aluminium,the EAM potential parametrised by[Winey,Kubota,and Gupta(2009)]was utilised and taken from[Becker(2011)].This potential reproduces the Aluminium FCC lattice constantthe cohesive energy Efcc=-2.646 eV,and the elastic constants d11=113.76 GPa&d12=61.71 GPa&

For BCC Tungsten,the EAM potential[Ackland and Thetford(1987);Becker(2011)]was used.This potential reproduces the Tungsten BCC lattice constantthe cohesive energy Ebcc=-8.91eV,and the elastic constants d11=532.61 GPa&d12=205.02 GPa&d44=163.20 GPa.

3.2 Uniform deformation tests from atomistic simulations

Using LAMMPS[Plimpton(1995)]in the Molecular Statics approach[Tadmor and Miller(2011);Maździarz,Young,and Jurczak(2011);Maździarz,Young,Dłużewski,Wejrzanowski,and Kurzydłowski(2010)]independently for FCC aluminium and BCC Tungsten,four numerical molecular homogenous deformation tests were performed,namely Uniaxial Strain(US),Simple Shear(SS),Volumetric Deformation(VOL)and Biaxial Strain(BS).The applied deformations correspond to the range of rank-one convexity of the energy function being equivalent to the validity of the Cauchy-Born hypothesis[Aghaei,Qomi,Kazemi,and Khoei(2009);Khoei,Qomi,Kazemi,and Aghaei(2009);Hudson(2011)].All tests were performed path forward and backward.Overlapping paths forward and backward indicate the elastic range of deformation.Each path of deformation was carried out in 50 steps.In each step,the iterations were terminated when the stopping criteria were satisfied,when the tolerance for energy and force reached etol=10-13&ftol=10-13.

3.2.1 Uniaxial strain

The deformation gradientFin crystallographic axes coinciding with Cartesian coordinate system axes for uniaxial strain without perpendicular deformations is defined by

3.2.2 Simple shear

The deformation gradientFin crystallographic axes for simple shear is classically defined by

where γ=tan(Φ),Φ is the angular change.The simulation box was sheared for aluminium by γ=25%and returned along the same path to the initial con figuration,see Table 8-9.

3.2.3 Volumetric deformation

The deformation gradientFin crystallographic axes for volumetric deformation is defined by

where λ=L/L0is the principal stretch/compression ratio.The simulation box was stretched volumetrically by 19%for aluminium,returned along the same path to its initial con figuration,then compressed volumetrically by 19%and again returned along the same path to the initial con figuration.For tungsten,the same test was performed,but the stretch/compression ratio achieved was only 11%,see Table 8-9.

3.2.4 Biaxial strain

In this case,the deformation gradientFin crystallographic axes for biaxial strain without normal deformation(Plane Strain)is defined by

4 The elastic constants from atomistic simulations

Using an orthonormal basis(Itskov,2009;Man and Huang,2011),Hooke’s relation for cubic symmetry can be written in the form

A popular measure of anisotropy is the so-called anisotropy ratio(Kambouchev,Fernandez,and Radovitzky,2007),de fined as

and for the analysed materials,it has valueswhich is quite close to isotropy.According to(Man and Huang,2011;Kowalczyk-Gajewska,2009)the isotropic average estimates for the elastic constants for FCC aluminium and BCC tungsten of Voigt,Reuss and Voigt-Reuss-Hill(VRH)were used,see Eq.(18)&Tables 2-3.It is worth mentioning that the elastic constants d12,d44correspond exactly to the Lamé constants λ0,µ0,respectively,while d11=λ0+2µ0.

Table 2:Elastic constants of FCC aluminium[GPa]and its isotropic average estimates[GPa].

Table 3:Elastic constants of BCC tungsten[GPa]and its isotropic average estimates[GPa].

5 Fitting parameters for hyperelasticity to data from atomistic simulations

On the basis of the atomistic calculations carried out in Section 3,the results for four homogeneous deformations:volume deformation,shear,uniaxial and biaxial deformation in the case of aluminium and tungsten for quite significant deformations were obtained.These results in this section will serve as data to determine the free(it should be emphasized that two of the parameters are determined and adopted so as to ensure that the model predictions for small deformations are the same as for Hooke’s classical relation)material parameters in the analysed constitutive models.

The parameters λ0and µ0are taken on the basis of the Voigt isotropic average estimates from Section 4.For the statistical evaluation of the obtained approximation on the basis of the given model and a set of material parameters,the cumulative determination coefficient CR(CR∈(0,1))was introduced("RSquared" coefficient of determination R2[Wolfram Research(2010)]),which has a standard interpretation,with the only difference that is calculated taking into account three out of the four tests(with equal weights),i.e.,the volumetric deformation test,shearing test,and uniaxial deformation test.

In the case of the SVK and LN models,for the material data as in Section 4,e.g.µ0=29160 MPa,λ0=59620 MPa for aluminium,see Table 2 and µ0=163440 MPa and λ0=205260 MPa for tungsten,see Table 3,the cumulative values of determination coefficients are given in Table 4.

Analysis of the results presented in Table 4 shows that the obtained approximations are very good in the range of deformation considered as elastic for aluminium and tungsten except for the SVK model approximation for aluminium.At this stage,however,it is worth indicating the disadvantages of these two basic models and suggesting others that are free of these defects.The main disadvantage of the SVK model mentioned earlier is the fact that during volumetric compression to volumes close to zero,the energy does not goes to infinity but takes a constant value,whichfrom the physics standpoint is wrong,see Fig.2b.On the other hand,the LN model is free from that defect,but its numerical implementation is difficult because of the need of calculating the logarithm of the tensor and its derivatives(Itskov,2009).Neither of the models meets the requirements of polyconvexity.

Table 4:Cumulative values of determination coefficient in case of SVK and LN models.

Figure 2:Contour graphs of SEF for plane strain case as a function of principal stretches for tungsten:a)for relatively small deformations,b)for large deformations(solid line-SVK model,dashed line-LN model)[MJ/m3].

Choosing models that are free from the defects of the LN and SVK models and the fulfilment of the polyconvexity requirements altogether with the requirement that the material parameters included in the model should be involved linearly,optionally with the exception of the volumetric function,were addressed.The requirement that the material parameters should be involved linearly in the SEF,is not fulfilled in the case of Agarwal’s function proposed in Table 1.The functions Γ(J)describing the volume changes,however,as will be seen later in this article,have a crucial impact on improving the quality of the approximation.

One of the first proposed hyperelasticity models of compressible isotropic materials meeting the polyconvexity requirement and situated in a class of Eq.(9)models with the functions of volume changes as in Table 1 is Ciarlet’s model with the SEF which can be brought to the form[Jemioło(2002)]

In this model,in comparison to the LN and SVK models,only one additional parameter f needs to be determined.In this paper,it is assumed that the free material parameters will be determined using the linear/non-linear regression procedure for three tests simultaneously(i.e.,volume deformation,shearing and uniaxial strain).The fourth test is used for further verification of the correctness of the determined approximation.However,it is not always possible or advisable to use the results of all tests,as for example in the case of Ciarlet’s model,the f parameter is only present in the SEF for volume and shearing deformation tests.Consideration of the other tests does not change anything.Of course,the determination of the parameter f on the basis of only the shearing test will provide an excellent approximation to that test while having a negative impact on the volume deformation test approximation,and vice versa(see the results presented in Table 5,where the cumulative determination coefficients are presented together with the determination coefficients for the single tests:RVOL,RSS,RUS,for both aluminium and tungsten).

Table 5:The values of the parameter f,the cumulative values of the determination coefficient,and the determination coefficient for each individual tests for Ciarlet’s model

Additionally it should be noted that even if the model is linear with respect to its parameters,but restrictions on their values are active,then from the numerical point of view it is necessary to use a non-linear regression procedure.In determining the parameter f for aluminium and tungsten,it was assumed that f∈ (0,1)and λ0＞µ0(1-f),due to the polyconvexity requirement for SEF,Eq.(20).In determining the material parameter f for aluminium and tungsten in two cases,the resulting parameter f was close to zero(the first constraint on the parameter had activated),which means that there is no influence of the first invariant,see Table 5.

It can be noticed that in the case of aluminium,comparable results were obtained determining the parameter f on the basis of the volumetric deformation test or on two tests simultaneously(i.e.the volumetric deformation test and the shear test),whereas in the case of tungsten,the best approximation of the atomic calculation results was obtained using data from two tests simultaneously.The cumulative determination coefficients in these cases are higher than in the case of the SVK model and somewhat lower than in the case of the LN model.Ciarlet’s model is free from the defects of the SVK and LN models,as shown in Fig.3a.In addition,in Fig.3b a contour graph of the SEF for data from the atomic calculations for aluminium(and the assumption of plane strain)for Ciarlet’s and LN models is presented.In this figure,fairly considerable differences for both models for the stretches up to 1.3 are visible.

Figure 3:SEF for plane strain case as a function of principal stretches for aluminium a)the cross-sections of SEF with λ1=λ2= plane,b)SEF’s contour graphs[MJ/m3].

Another group of constitutive models,which are a special case of the function(9),are the Neo-Hookean(NH)type models with SEF in the form

The application of a function of volume changes as in the second row of Table 1,i.e.,the function proposed by Agarwal,meeting the requirements of the natural state and consistency with SVK model leads to the following relations on the parameters(the same for all models of the NH type),c1= µ0andIn the case of the function of volume changes proposed by Simo and Pister,the parameters are equal to,respectively,and c2=µ0.In the third case,i.e.,Ogden’s function of volume changes,it is possible to obtain the relationsand

Since,for each of the NH type models,obtaining its consistency with the SVK model and the fulfilment of its natural state assumption leads to the conclusion that the parameters,as for the linear Hooke’s relation,entirely determine the material behaviour,so also for the material parameters of aluminium and tungsten obtained in Section 4,the evaluation of the approximation to the experimental results as described previously was performed.The resulting cumulative values of determination coefficient are given in Table 6.

Table 6:Cumulative value of determination coefficient for the NH model with different functions of volume changes.

From the results presented in Table 6 it can be concluded that in all cases the best approximations were obtained with Agarwal’s function describing the material volume changes.Additionally,it can be stated that the test results obtained on the basis of the atomic calculations in the case of tungsten are easier to approximate.It is worth noting that the NH-type models do not have the defects of the LN and SVK models,and the cumulative determination coefficients in all cases have higher values than in the case of the hyperelasticity models previously analysed.An additional advantage of these models is that they can be used,for example,in the analysis of the metals in boundary-value problems without having to determine any additional material parameters than in the linear elasticity theory with the classical Hooke’s relation.Fig.4a-c present the graphs of the first Piola-Kirchhoff tensor S11component in:volumetric,uniaxial and biaxial tests for different NH model variants and SVK model for aluminium,and hence different values of cumulative determination coefficients.The graphs of the first Piola-Kirchhoff stress tensor S11component of the SVK model con firm its defects,particularly evident in Fig.4c,where for λ≈0.6,the value of the S11component stops growing.It is worth noting that the differences in stress prediction using the NH type model with different Γ(J)functions for the largest predicted stretches are significant.In addition,Fig.4d compares the NH model with the predictions using Agarwal’s volume changes function,for three individual tests(all data taken as for aluminium).

NH type models do not include the second invariant of the stretch tensor.The simplest model in which there is a dependence on the second invariant is the Mooney-Rivlin model(MR),which can be obtained by assuming M=N=1 in Eq.(9).Regarding the function of volume changes in Eq.(9),on the basis of the results obtained for the NH model it was assumed that in the further considerations always an Agarwal’s function is chosen.After applying the assumption of the consistency with the SVK model,only one parameter in MR remains to be determined in the optimization procedure.It was assumed that this parameter is k.

Figure 4:The component S11of the first Piola-Kirchhoff stress tensor in tests with:a)volume,b)the uniaxial c)biaxial deformations for different NH model variants and d)comparison of the response of NH model with Agarwal’s volume changes function for the three individual tests(data taken as for aluminium)[MPa].

This parameter was determined on the basis of multi-criteria optimization(i.e.,matching the three tests at the same time with the limitations on parameters that a1＞0,b1＞0,and k＞0.5,which provides SEF polyconvexity).In the case of aluminium,the parameter k=1.022 was obtained,which in turn leads to b1→0,and thus to eliminating the influence of the second invariant.It is different in the case of tungsten,where k=0.5552,and the influence of the second invariant is visible.It should be noted,however,that the cumulative determination coefficient in this case has only increased by 1.3%compared to the NH model.

Figure 5:Graphs of the SEF functions for different NH model variants for tungsten(line markings as shown in Fig.4)[MJ/m3].

Table 7:Cumulative value of determination coefficient for different MRL models with Agarwal’s function of volume changes.

Further generalization of the MR model is done by increasing the number of elements in the series(9).In order to distinguish between these models,it was assumed that,for example for M=N=4,the Mooney-Rivlin like(MRL)model will be denoted by MR4.For MRL models as before,Agarwal’s function of the volume change is assumed,because it was considered as most appropriate.It should be noted that in the MR2 model(and higher)it is possible to approximate the shear test better,because in all models analysed so far in the shearing test,onlyµ0was present.However,caution must be exercised,because an improvement in the shear test predictions can result in significant deterioration in the predictions of the other tests,see Fig.6a.This figure shows,in the case of aluminium,the predictions of the MRL model with two sum elements of the shearing test,for which the material parameters are obtained in two cases:using the data for the shear test,or using all three tests.This allows a better approximation of the shear test,but in the case of the other tests the predictions are completely wrong,because the parameter k→0 creates a singularity in the volume part of the SEF.In turn,Fig.6b shows how the number of sum terms taken in the MRL model affect the quality fo the approxi-mation of the SEF,which is also measured by the CR coefficient and presented in the case of the volumetric deformation test for tungsten.In this case,the material parameters were determined analogously as before on the basis of three tests at the same time with all the necessary constraints.It can be seen that the best approximations of the results from the atomistic simulations can be obtained by using the MR4 model,although the CR increases only 1.3%with respect to the MR model,see Table 7.The contour graph presented in Fig.7 shows the SEF for the plane strain case and MRL models with two and four sum elements,for moderate stretches(i.e.,in the range in which the atomistic calculations were made,providing the data to determine the material parameters),and in the case of very large deformations.In the range of stretches similar to those achieved in the atomistic calculation,the compatibility of the MR2 and MR4 models is very good,but outside of this range,significant differences between the potentials are visible.

The optimized parameters for the MRL models for aluminium and tungsten and the correspondence between the parameters obtained on the basis of the natural state restriction with Hooke’s relation can be found in Table 10 and Eq.(22).

Figure 6:Graphs of the SEF functions for different MRL model variants:a)for aluminium shearing test in case of MR2,b)for tungsten volumetric deformation test for subsequent MRL models[MJ/m3].

6 Conclusions

Figure 7:Contour graphs of SEF for plane strain case as a function of principal stretches for tungsten:a)for relatively small deformations,b)for large deformations(solid line-MR2 model,dashed line-MR4 model)[MJ/m3].

Table 8:Calculated strain energy density data curves:Volumetric Deformation,Simple Shear,Uniaxial Strain and Biaxial Strain for Aluminium.

Table 9:Calculated strain energy density data curves:Volumetric Deformation,Simple Shear,Uniaxial Strain and Biaxial Strain for Tungsten.

Table 10:Parameters for MRL models for Aluminium and Tungsten.

Analysing the results of atomistic calculations, some differences between monocrystalline aluminium and tungsten can be observed.The first material exhibits a higher degree of anisotropy,but still it can be regarded as close to being isotropic.Therefore,to approximate the results of atomistic calculations,the hyperelasticity constitutive relations for isotropic materials were applied.These relations were selected to meet the fundamental physical and mathematical requirements and also to be relatively simple and compatible with the linear Hooke’s relation for isotropic materials for the strain tensor,e.g.,E→0.Hyperelastic isotropic polyconvex SEF functions of general form Eq.(9)were fitted to the strain energy function coming from molecular simulations of aluminium and tungsten.The resulting approximations based on the simultaneous matching of the three independent tests have been evaluated,by determining the cumulative correlation coefficient.It was shown in examples that it is always possible to obtain a better fit of a hyperelasticity constitutive model for a single test,but at the expense of a significant deterioration of the approximations to the remaining tests.In summary,it can be deduced that:

1.Among the simplest models, the logarithmic one(LN)is superior to the SVK,but is unfortunately difficult to numerically implement.There is evident,however,a clear difference in the cumulative coefficients of determination between the aluminium and tungsten.In both cases,the SVK model gives a less accurate approximation to the results of atomistic simulations than logarithmic one.

2.Polyconvex SEF functions do not give worse results than the simplest SEF functions.Furthermore,the simplest polyconvex ones,the Neo-Hookean(NH)type,do not even need fitting any extra parameters beyond the regular Lamé constants from the small deformation theory.

3.The quality of the SEF is decided mainly by the proper choice of volumetric function,not the choice of the number of elements in the series(9)to be taken.Therefore,it is sometimes worth abandoning the assumption that the material parameters should enter the model linearly,as was the case in the Agarwal function considered here.Keep in mind that the process of optimizing the parameters should stay controllable.

4.For fully compressible SEF functions,there is no volumetric-isochoric split,furthermore the explicit dependence on J imposes no restrictions on the scope of their application as opposed to the model proposed in[Henann and Anand(2011)].

5.For the polynomial type of polynonvex SEF functions,it is easy to fit their parameters and implement them in Finite Element Method codes,even commercial ones,e.g.ABAQUS through the user subroutine UMAT or UHYPER[Hibbit,Karlsson,and Sorensen(2014)].

Ackland,G.J.;Thetford,R.(1987):An improved N-body semi-empirical model for body-centred cubic transition metals.Philosophical Magazine,Part A,vol.56,pp.15-30.

Aghaei,A.;Qomi,M.A.;Kazemi,M.;Khoei,A.(2009): Stability and sizedependency of Cauchy-Born hypothesis in three-dimensional applications.International Journal of Solids and Structures,vol.46,no.9,pp.1925-1936.

Atluri,S.N.(1984): Alternate stress and conjugate strain measures,and mixed variational formulations involving rigid rotations,for computational analyses of finitely deformed solids,with application to plates and shells-I:Theory.Computers and Structures,vol.18,no.1,pp.93-116.

Becker,C.(2011): Atomistic simulations for engineering:Potentials and challenges.Tools,Models,Databases and Simulation Tools Developed and Needed to Realize the Vision of ICME,ASM.

(2002): On the Ellipticity of Finite Isotropic Linear Elastic Laws.Preprint.Univ.,Fak.für Maschinenbau.

(2007):Coordinate-free Characterization of the Symmetry Classes of Elasticity Tensors.Journal of Elasticity,vol.87,no.2-3,pp.109-132.

Bonet,J.;Wood,R.D.(2000): Nonlinear continuum mechanics for finite element analysis.Cambridge University Press.

Bruhns,O.T.(2001): Constitutive inequalities for an isotropic elastic strainenergy function based on Hencky’s logarithmic strain tensor.Proceedings of the Royal Society A:Mathematical,Physical and Engineering Sciences,vol.457,no.2013,pp.2207-2226.

Ciambella,J.;Saccomandi,G.(2014): A continuum hyperelastic model for auxetic materials.Proceedings of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,vol.470,no.2163.

Ciarlet,P.(1988): Mathematical Elasticity Volume I:Three-Dimensional Elasticity.North-Holland,Amsterdam.

Clayton,J.D.(2011): Nonlinear Mechanics of Crystals.Springer Dordrecht Heidelberg London New York.

Danielsson,M.;Parks,D.;Boyce,M.(2004):Constitutive modeling of porous hyperelastic materials.Mechanics of Materials,vol.36,no.4,pp.347-358.

Dimitrienko,Y.I.(2002): Tensor Analysis and Nonlinear Tensor Functions.Springer Netherlands.

Dłużewski,P.(2000):Anisotropic Hyperelasticity Based Upon General Strain measures.Journal of elasticity and the physical science of solids,vol.60,no.2,pp.119-129.

Hartmann,S.(2001):Parameter estimation of hyperelasticity relations of generalized polynomial-type with constraint conditions.International Journal of Solids and Structures,vol.38,no.44-45,pp.7999-8018.

Hartmann,S.;Neff,P.(2003):Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility.International Journal of Solids and Structures,vol.40,no.11,pp.2767-2791.

Henann,D.L.;Anand,L.(2011): A Large Strain Isotropic Elasticity Model Based on Molecular Dynamics Simulations of a Metallic Glass.Journal of Elasticity,vol.104,no.1-2,pp.281-302.

Hibbit;Karlsson;Sorensen(2014): ABAQUS/Standard Analysis User’s Manual.Hibbit,Karlsson,Sorensen Inc.,USA.

Hill,R.(1968):On constitutive inequalities for simple materials-I.Journal of the Mechanics and Physics of Solids,vol.16,no.4,pp.229-242.

Hudson,Thomas,O.C.(2011): On the stability of Bravais lattices and their Cauchy-Born approximations*.ESAIM:Mathematical Modelling and Numerical Analysis,vol.46,no.1,pp.81-110.

Itskov,M.(2009): Tensor Algebra and Tensor Analysis for Engineers:With Applications to Continuum Mechanics.Springer Publishing Company,Incorporated,2nd edition.

Jemioło,S.(2002):Studium hiperspreżystych własności materiałów izotropowych.Prace Naukowe Politechniki Warszawskiej.Budownictwo,vol.z.140,pp.3-308.

Kambouchev,N.;Fernandez,J.;Radovitzky,R.(2007):A polyconvex model for materials with cubic symmetry.Modelling and Simulation in Materials Science and Engineering,vol.15,no.5,pp.451.

Khoei,A.;Qomi,M.A.;Kazemi,M.;Aghaei,A.(2009):An investigation on the validity of Cauchy-Born hypothesis using Sutton-Chen many-body potential.Computational Materials Science,vol.44,no.3,pp.999-1006.

Kowalczyk-Gajewska,K.(2009):Bounds and self-consistent estimates of overall properties for random polycrystals described by linear constitutive laws.Archives of Mechanics,vol.61,no.6,pp.475-503.

Kueh,A.B.H.(2012): Fitting-free hyperelastic strain energy formulation for triaxial weave fabric composites.Mechanics of Materials,vol.47,pp.11-23.

łepkowski,S.P.;Majewski,J.A.;Jurczak,G.(2005):Nonlinear elasticity in III-N compounds:Ab initio calculations.Phys.Rev.B,vol.72,pp.245201.

Ling,X.;Atluri,S.(2007): A hyperelastic description of single wall carbon nanotubes at moderate strains and temperatures.CMES-Computer Modeling in Engineering and Sciences,vol.21,no.1,pp.81-91.

Man,C.-S.;Huang,M.(2011):A Simple Explicit Formula for the Voigt-Reuss-Hill Average of Elastic Polycrystals with Arbitrary Crystal and Texture Symmetries.Journal of Elasticity,vol.105,no.1-2,pp.29-48.

Maździarz,M.;Young,T.D.;Dłużewski,P.;Wejrzanowski,T.;Kurzydłowski,

K.J.(2010):Computer Modeling of Nanoindentation in the Limits of a Coupled Molecular-Statics and Elastic Scheme.Journal of Computational and Theoretical Nanoscience,vol.7,no.6,pp.1172.

Maździarz,M.;Young,T.D.;Jurczak,G.(2011):A study of the affect of prerelaxation on the nanoindentation process of crystalline copper. Archives of Mechanics,vol.63,no.5-6,pp.533.

Murdoch,A.(2007): A Critique of Atomistic Definitions of the Stress Tensor.Journal of Elasticity,vol.88,no.2,pp.113-140.

(2013):Deficiencies in numerical models of anisotropic nonlinearly elastic materials. Biomechanics and Modeling in Mechanobiology,vol.12,no.4,pp.781-791.

Ogden,R.W.(1984): Non-Linear Elastic Deformations.Dover Publications.

Ogden,R.W.;Saccomandi,G.;Sgura,I.(2004):Fitting hyperelastic models to experimental data.Computational Mechanics,vol.34,no.6,pp.484-502.

Ostrowska-Maciejewska,J.;Kowalczyk-Gajewska,K.(2013):Rachunek tensorowy w mechanice ośrodków ciagłych.Biblioteka Mechaniki Stosowanej:Monogra fie.Wydawnictwo Instytutu Podstawowych Problemów Techniki PAN.

Pandolfi,A.;Vasta,M.(2012): Fiber distributed hyperelastic modeling of biological tissues.Mechanics of Materials,vol.44,no.0,pp.151-162.

Paufler,P.(2003): International Tables for Crystallography.Volume D:Physical Properties of Crystals.by André Authier(Editor).Published for the International Union of Crystallography by Kluwer Academic Publishers:Dordrecht/Boston/London,2003. Acta Crystallographica Section A Foundations of Crystallography.

Peng,S.T.J.;Landel,R.F.(1975):Stored energy function and compressibility of compressible rubberlike materials under large strain.Journal of Applied Physics,vol.46,no.6,pp.2599-2604.

Plimpton,S.(1995): Fast Parallel Algorithms for Short-Range Molecular Dynamics.Journal of Computational Physics,vol.117,no.1,pp.1-19.

Salvetti,M.(2010): Hyperelastic Continuum Modeling of Cubic Crystals Based on First-principles Calculations.Massachusetts Institute of Technology,Department of Mechanical Engineering.

Sansour,C.(2008): On the physical assumptions underlying the volumetricisochoric split and the case of anisotropy. European Journal of Mechanics-A/Solids,vol.27,no.1,pp.28-39.

Schroeder,J.;Neff,P.;Ebbing,V.(2008):Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors.Journal of the Mechanics and Physics of Solids,vol.56,no.12,pp.3486-3506.

Shen,S.;Atluri,S.N.(2004):Atomic-level Stress Calculation and Continuum-Molecular System Equivalence.CMES-Computer Modeling in Engineering and Sciences,vol.6,no.1,pp.91-104.

Silhavy,M.(1997): The mechanics and thermodynamics of continuous media.Springer;1997 edition(November 28,1996).

Steinmann,P.;Hossain,M.;Possart,G.(2012): Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for Treloar’s data.Archive of Applied Mechanics,vol.82,no.9,pp.1183-1217.

Tadmor,E.;Miller,R.(2011): Modeling Materials:Continuum,Atomistic and Multiscale Techniques. Modeling Materials: Continuum, Atomistic, and Multiscale Techniques.Cambridge University Press.

Valavala,P.;Clancy,T.;Odegard,G.;Gates,T.(2007):Nonlinear multiscale modeling of polymer materials. International Journal of Solids and Structures,vol.44,no.3-4,pp.1161-1179.

Vergori,L.;Destrade,M.;McGarry,P.;Ogden,R.(2013):On anisotropic elasticity and questions concerning its finite element implementation.Computational Mechanics,vol.52,no.5,pp.1185-1197.

Wang,J.;Mao,Z.;Jiang,F.;Duffy,T.(2014):Elasticity of single-crystal quartz to 10 GPa.Physics and Chemistry of Minerals,pp.1-10.

Winey,J.M.;Kubota,A.;Gupta,Y.M.(2009): A thermodynamic approach to determine accurate potentials for molecular dynamics simulations:thermoelastic response of aluminum. Modelling and Simulation in Materials Science and Engineering,vol.17,no.5,pp.055004.

Wolfram Research,I.(2010):Mathematica,2010.

Xu,M.;Paci,J.T.;Oswald,J.;Belytschko,T.(2012):A constitutive equation for graphene based on density functional theory.International Journal of Solids and Structures,vol.49,no.18,pp.2582-2589.

Zheng,Q.S.(1994):Theory of Representations for Tensor Functions-A Unified Invariant Approach to Constitutive Equations.Applied Mechanics Reviews,vol.47,no.11,pp.545-587.

Zhou,M.(2003): A new look at the atomic level virial stress:on continuummolecular system equivalence. Proceedings of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,vol.459,no.2037,pp.2347-2392.