Z.D.Han and S.N.Atluri

1　Balance laws for the Cauchy Stress σ,the first Piola-Kirchhoff Stress P,the second Piola-Kirchhoff Stress S,and the Eshelby Stress T

We consider the finite deformation of a solid,wherein a material particle initially atX,moves to a locationx.We use a fixed Cartesian coordinate system with base vectorsei,such that(X=XIeiandx=xiei).The displacement of the material particle is,andThe deformation gradient tensor isThere are in finitely many possible definitions of a stress-tensor in a finitely deformed solid[see,for instance Atluri(1984)].Among the more commonly used ones are:the Cauchy stress tensorσ;the first Piola-Kirchhoff stress tensorP;and the second Piola-Kirchhoff stressS,which are related to each other,thus[see Atluri(1984)]:

whereJis‖F(xiàn)‖,and(·)tdenotes a transpose.

Considering a generalanisotropic hyperelasticsolid,with the strain energy per unit initial volume being denoted asW,the constitutive relation forPmay be written as[see Atluri(1984)]:

IfWis a frame-indifferent function ofF,it should be a function only ofFt·F.Thus[see Atluri(1984)],

where

asCis the right Cauchy-Green deformation tensor,Ebeing the Green-Lagrange Strain tensor.

The equations of Linear Momentum Balance(LMB)and Angular Momentum Balance(AMB)can be written equivalently in terms ofσ,P,andS[see Atluri(1984)],as:

whereρ0is the mass density(per unit initial volume).For a homogeneous solid,ρ0is not a function ofX,but is a constant.

The equivalence of Eqs.(6)a and(6)b can be seen from the following[using Eq.(1)],

However,as noted by Shield(1967)and Ogden(1975),we have the purely geometric identity for any finite deformation,that:

Thus,from Eqs.(6)(d,e),we have:

Thus,the geometric identity(6)e guarantees the equivalence of Eqs.(6)a and(6)b,in any exact solution.However,in any computational solution,such an equivalence may not always be assured.

WhenFis derived from an objectiveW,as in Eq.(4),P≡S·Ftby definition[Sis symmetric],and hence the AMB forFin Eq.(6)b is inherently embedded in the structure ofW,viz,thatWis a function ofEonly.

Following the seminal work of Noether(1918),and the extraordinarily important contributions of Eshelby(1957,1975),the Eshelby stress tensor is de fined,for finite elasto-static deformations,as

or

whereIis an identity tensor,Wis the strain energy density(per unit initial volume),Pikis the first Piola-Kirchhoff stress tensor,anduk,J=?uk/?XJ.WhileTis often referred to as the “Energy-Momentum Tensor”,it clearly has the dimensions of“Stress”and was also independently derived by Atluri(1982)[see Eq.(29)in Atluri(1982)].It was also discussed in Atluri(1984)that any function ofPandFis a stress-measure(in finite deformation solid mechanics),such asTwhich is also a function ofPandF.We may also w riteT,equivalently,as[see Eq.(27)in Atluri(1982)]:

It can be seen from Eqs.(7)a and(8)thatTis in general an unsymmetric tensor,unlessSandCare co-axial,which is the case only for isotropic materials.We note thus,that the stress tensorTis symmetric for isotropic materials.We note thatPis a two-point tensor field such that

where(NdA)is a vector of an oriented-area in the undeformed con figuration anddfis a vector of force acting on an oriented area(nda)(the image of(NdA))in the deformed con figuration.Likewise,Fis a two-point tensor field such that

whereGLare covariant base vectors in the undeformed con figuration,andgkare covariant base vectors in the deformed con figuration.

Hence it follows thatP·Fisa tensorde fined in the undeformed con figuration.SinceWis the strain energy density per unit volume in the undeformed con figuration,it is seen that the tensorTin Eq(7)isa tensor defined entirely in the undeformed con figuration,similar to the second Piola-Kirchhoff stress tensorS,except thatSis symmetric,andTis unsymmetric.

If we consider in finitesimal deformations,then the difference betweenP,Sandσdisappears.Further,if we consider linear elastic behaviors,along with infinitesimal deformation,P(orSorσ)become linear inF.Even for linear elastic in finitesimal deformations,the Eshelby Stress TensorTis a quadratic function ofF.This poses some difficulties in the computational mechanics of in finitesimal deformations of even linear elastic solids,if the Eshelby stress tensor is used directly as a variable.The “traction vector”t?,corresponding to the stress-tensorT,at a boundary with a unit normalNin the initial con figuration,is written as:

or

Figure 1:a close volume V includes a crack tip w ithin Vε

such that,the integral of the traction vectort?over the small volumeVεenclosing a crack-tip in 2-D,or a crack-segment in 3-D[see Fig.1],may be denoted as the vectorT?,and represented thus:

or

Follow ing the celebrated works of Eshelby(1957,1975),we callT?the vector of“Force on the Defect”in finite-deformation elasto-statics of general(anisotropic)hyperelastic solids.More precisely,in the case of cracks,it is a vector that quantifies the singular nature of stress and strain fields inVε.We note thatis free from any defects and singularities,and we may w rite:

in which the follow ing equations are used:

Thus,for the problem of a hyperelastic solid,containing a crack-like defect1In the case of a sharp-crack-like defect,in the integrand in T?in?Vε,on the right hand side of Eq.(14)a can have singularities of the type r?1/2,such that,in a 2-D problem,ds=rdθ,and T?has a finite value even when r→ 0 in Vε.,the vector of“concentrated Force on the Defect”T?is given in Eq.(12).For an elastostatic problem of a defective homogenous solid,with no body forces or crack face tractions,we have a “path-independent”representation forT?:

This forceT?in Eq.(12)has been called the “Force on the Defect”by Eshelby(1951,1957,1975).

In as much as the Stress TensorT[Eq.(7)]is de fined as:

We see that in[Eq.(13)],even for finite deformations,

If the material is homogeneous,[i.e.,Wdoes not explicitly depend onX],and the body forces do not exist,and if the deformation is static,the stress tensorTis divergence free in a volumeV?Vε,which is free of singularities:

Otherwise,in general,we may w rite the generalized balance law forT:

Eq.(19)a implies that the balance law forTinherently involves the weak-form of the linear momentum balance law forP,multiplied by the“test function”F.

Eq.(19)a may be written more conveniently as,

where the vectorbis “the distributed force on the Defect”.Thus,each of the 3 balance lawsTis equal to a combination of the 3 balance laws forP.If the material is hypo-elastic or elastic-plastic,we may consider an objective rate of the stresstensorT,and an incremental“Strength of the Singularity”,ΔT?,as contemplated in Atluri(1982).

It is interesting to observe that,the “generalized weak-form”,

whereVis a volume free of singularities and defects,andis the virtual deformation gradient of a comparison state with displacement fieldhas been used by[Okada,Rajiyah and Atluri(1989a,b),Han and Atluri(2003),and Qian,Han and Atluri(2004)]in deriving very novel non-hyper-singular integral equations for stresses in solid-mechanics,and gradient fields in acoustics,etc.These novel non-hyper-singular integral equations lead to extremely convenient algorithms for traction boundary value problems,and 3-dimensional fracture and fatigue mechanics problems[Han and Atluri(2002)].

Tis in general anunsymmetric tensor defined in the undeformed con figurationfor anisotropic materials,and the balance law forT,even in finite deformation as stated in Eq.(19)is a set of linear PDEs in the undeformed coordinates,XI.It has been stated earlier that,for isotropic materials,Tis symmetric and,similar toS,is a tensor entirely defined in the undeformed con figuration.Thus,for in finitesimal deformations in an isotropic solid,Tstill satisfies the balance law:

Fora general anisotropic hyperelastic solid,one may find the displacement-like quantitiesvI.We may requirevKto satisfy the same boundary conditions as foruiat?(V?Vε).ThusvKmay be thought of as“the variation of the undeformed body”which is compatible with the prescribed boundary conditions forui,and correspond to the compatible strains derived from the stress tensorT.If we assume the displacementu=u(X)is one-to-one mapping,its inverse mapping can be taken as the displacements[Knowles and Sternberg(1972)],as,

Thus the deformed con figuration has been mapped back to the undeformed con figuration.In the other words,vKcan also be considered as“the displacements of the deformed body”.

For finite deformations,the inverse deformation gradients,F?1,is de fined as,

One may define the strain energy in the deformed con figuration,denoted as,as a frame-indifferent function ofF?1.The corresponding first Piola-Kirchhoff stress tensor of the inverse deformation,denoted ascan be de fined accordingly,as done in Eq.(3),as,

The Eshelby stress tensorTcan also be de fined alternately[Eshelby(1975),Chadw ick(1975)],as

It shows the duality of the Eshelby Stress and the Cauchy Stress which is discussed in Section 3.

The corresponding left Cauchy-Green deformation tensor of the inverse deformation gradients can be w ritten as,

which is also the inverse of the right Cauchy-Green deformation tensor.It has been addressed as the “Finger Deformation Tensor”in the chemistry community for handling various physical fields in the current con figuration,which is not common in applied mechanics.The corresponding strain tensor can be de fined as

One may observe that for a linear elastic material,

whereEbeing the Green-Lagrange Strain tensor in Eq.(5).

2　A variety of weak-form s for the balance laws for T,in finite-deformation,anisotropic hyperelastic ity

For finite deformations of a non-isotropic and non-homogeneous hyperelastic solid as shown in Fig.1,the balance laws for the unsymmetric Eshelby stress tensorTare:

which are a set of strong-form linear partial differential equations in the undeformed coordinatesXIin the initial con figuration in which the unsymmetic tensorTis entirely defined.

If we choose arbitrary but differentiable test functionswJ(X),we may w rite the weak-form of Eq.(28),as:

SinceV?Vεis free of any singularities,an application of the divergence theorem to Eq.(29)leads to:

On the other hand,if we choose the gradientswJ,K(X)as the test functions,the weak-form of Eq.(28)may be written in a vector-from,as:

Or,equivalently,as:

Now we consider a class of test functions in Eq.(30),such that2We already considered the general test functioncorresponding to a virtual displacement fieldin deriving non-hyper-singular integral equations for stresses in solid mechanics,in Eq.(20).,

whereBJKandCJare constants.

We now consider several simple cases of the constantsBJKandCJ.

Case(A):BJK=0and

In this case,

Use of Eq.(34)in Eq.(30)results in:

For arbitraryCJ,Eq.(35)leads to:

In Knowles and Sternberg(1972),and in Eshelby(1975),the notion of a“conservation law”is used only when the integral over the volumeV?Vεis zero;thus in the view of Knowles and Sternberg(1972)and Eshelby(1975),“path-independentintegrals”are identities expressed as integrals over the surface?(V?Vε)only.However,in a computational sense,we allow here that a path-independent-integral may involve both surface and volume integrals;such that an integral over?Vεmay be,for computational purposes,expressed as an integral over?Vplus another integral overV?Vε[See Nikishkhov and Atluri(1987)].In the present sense of a path-independent-integral,for finite deformations of nonhomogeneous nonisotropic solids,the “path-independent”integrals have the representation:

or

For elasto-static problems and for homogeneous anisotropic materials with zero body forces,Eq.(37)reduces to:

Case(B):andCJ=0.

We consider the test functionwJsuch thatBJKis skew-symmetric,and can be expressed as:

whereeJLKis the permutation tensor,such that

The weak form of Eq.(30)can now be written as:

Eq.(40)a may be written for arbitraryωLas:

Thus,for finite deformations of an anisotropic hyperelastic material,with body forces,the generalized L-Integral has the representation:

If one considers only an isotropic material,TIJ=TJIand thusTMJeJLM=0.Further,if one restricts to infinitesimal elasto-statics,and zero body forces,the integrands in the volume integral overV?Vεin Eq.(41)vanishes,the generalized L-Integral can be reduced as,

The L-Integral in Eq.(41)is equivalent to its alternative form in term of the displacements,as given by Knowles and Sternberg in 1972.By definition,Eq.(40)a can be written as,

However,the angular momentum balance for the stressP[NdS·P=dfin the deformed con figuration]states that:

or

Thus,Eq.(42)leads to:

Using Eq.(43)into Eq.(40)c,we obtain:

Eq.(44)is the conservation law in terms of the displacements for finite deformations of a hyperelastic,anisotropic solid.Eq.(44)can be reduced to the original form identified by Knowles and Sternberg in 1972,for small deformation linear elasticity without body forces,

Thus,with Eq.(44),the generalized L-Integral in Eq.(41)has the alternative representation in terms of the displacements,as

Case(C):BJK=CδJKandCJ=0.

such that

whereCis a constant.Use of Eq.(47)in Eq.(30)results in a scalar identity:

Thus,for finite deformations of an anisotropic hyperelastic material,with body forces,the generalized M-Integral has the representation:

The conservation law in Eq.(48)and the generalized M-Integral in Eq.(49)may also have the representation in term of the displacements.Eq.(48)may be written as:

However,

Thus,one is lead to the general conservation law:

where

For semi-linear anisotropic hyperelastic materials,one may postulateWsuch that:

Thus,for semi-linear anisotropic hyperelastic materials,one may write:

For semi-linear anisotropic materials,one may write the conservation law for finite deformations:

However,

Thus,Eq.(55)a may be written,for finite deformations,of sem i-linear anisotropic hyperelastic materials,as:

If we consider only in finitesimal deformations of anisotropic linear-elastic materials,and when deformations are independent of time and body forces are absent,Eq.(56)reduces to:

where,for in finitesimal deformations,PLkbecomes synonymous with the Cauchy stress tensor,i.e.,σik.Eq.(57)has been identified as a conservation law,leadingto the now so-called Mintegral given by Know les and Sternberg(1972)for infinitesimal deformations of linear elastic anisotropic materials.

Thus,the generalized M-Integral in Eq.(49),for finite deformations of an anisotropic hyperelastic material,with body forces,has the representation in term of the displacements.

Case(D):andCJ=0.

Wenow consider test functionswJ=BJKXK,whereBJKis asymmetric deformation matrix,with the constraint that

Thus,

With the polar decomposition,the Eshelby stress tensorTcan be written as

whereHIJ,I=TIJ,I=ρ0bJis curl-free andGIJ,I=0 is divergence free.With zero body forces and using Eq.(59)in Eq.(30),we may write:

Thus we obtain the generalized conservation law for finite deformation in ananisotropic hyperelastic solid,as:

By applying Stoke’s Theorem to Eq.(61),the generalized conservation law can be written for the potential function of the Eshelby stress tensor,as

Thus Eq.(63)may be considered to lead to the G-Integral:

Note that by using other arbitrary test functions,which may be arbitrary polynomials inXK,we may obtain an arbitrary number of generalized conservation laws.However,each of these four special cases(A-D)discussed above has its own physical meaning and is corresponding to its own conservation law,as

?Case A,the linear momentum conservation law;

?Case B,:the angular momentum conservation law;

?Case C,:the divergence theorem;

?Case D,:Stokes’theorem;

Note that the G-Integral is the fourth path-independent integral of Noether’s type[Noether(1918)],besides the three integrals reported by[Know les and Sternberg(1972)and Eshelby(1975)].Here these four conservation laws provide 12 independent equations for 3-D problems,and 6 independent equations for 2-D problems.has been w idely applied to incompressible problems with the pressure as an independent variable.can be applied for the independent shear stresses for problems with materials such as liquid crystals or meta-materials,or anti-plane problems assuming no shear strains.

In addition,the Eshelby stress tensor can also be extended to the gradient theory of solids if the strain energy function in Eq.(3)is also dependent on the second derivatives of the displacements,ui,JK[Eshelby(1975)].The corresponding conservation laws can also be obtained simply by writing the corresponding weak forms of both macroscopic and microscopic momentum balance laws,as well as their forms weighted by the “test function”of the first and second derivatives of the displacements,following Eq.(17).The concept of“force of defects”can also be determined through conservation laws[Gurtin(2002)].Another extension of the Eshelby tensor is to micropolar materials,in which the strain energy function is dependent on the deformation curvature.Similar conservation laws can also be derived by involving the curvature terms[Lubarda and Markenscoff(2003)].

3　The use of the Eshelby Stress Tensor in com putational finite deformation solid mechanics

Since the concept of“the force on defect”was introduced by Eshelby in 1951,the Eshelby stress tensor(or the energy-momenturm tensor)has been extended for continuum mechanics of solids by Eshelby(1975),and independently in Atluri(1982)in which the expression of the Eshelby tensor has been given in term of the strain energy in the deformed con figuration.The duality of the Eshelby Stress and the Cauchy Stress was also discovered by Eshelby(1975)and Chadwick(1975),and was extended to finite strain.With the use of the states of inverse deformation[Shield(1967)]and of the dual reciprocal states in finite elasticity[Ogden(1975)],the Eshelby stress tensor has been widely explored.The concept of“Eshelbian Mechanics”was introduced by Maugin(1995)based on the con figuration invariance of the energy conservation law of Noether’s Theorem,as the Lagrangian-Ham iltonian-Noetherian formulation.In contrast,Newtonian mechanics is based on the conservation laws of linear and angular momenta,and leads to the energy conservation law in the rate form in the undeformed con figuration.As shown in Eqs.(17)and(19),the balance law forTinherently involves the“weighted”weak-form of the momentum balance law forP.

The Eshelby stress tensor and its alternate forms have been widely used in developing numerical methods,especially for problems with singularities or inhomogeneities.On the other hand,by its definition in Eq.(7),the Eshelby stress tensor is a quadratic function of the deformation gradient tensor even for linear elastic or seme-linear(involving a linear relation betweenPandF)materials undergoing in finitesimal deformations.It becomes very difficult to develop numerical methods explicitly based on using the Eshelby stress tensor as a direct variable in solving problems,sinceTis a nonlinear function ofFeven for small-strain linear elastic behavior.A few exceptions include the exact use of the Eshelby tensor for inhomogeneous inclusions(Eshelby 1957,1959),and in the boundary integral forms for micromechanics[Mura(1991)].The Eshelby stress tensor has so far been widely used only ina post processing computationto evaluate the forces on defects(or the con figurational forces),especially for mesh-based numerical methods,once the stress and deformation are already computed.

We now review the difference between the equilibrium equations of the Cauchy stress,and the balance laws for Eshelby stress.The strong form of the momentum balance equations for the Cauchy stress tensorσare in the deformed con figuration as in Eq.(6)a.The Cauchy stress tenorσin the deformed con figuration is analogous to the Eshelby stress tensor in the undeformed con figuration,and the Cauchy stress tensor is also a quadratic function of the inverse of the deformation gradient tensor,i.e.even for semi-linear elastic solids.It implies that both the equilibrium equations based on the Cauchy stress tensor and the balance law of the Eshelby stress tensor are not suitable for linearization for solving the corresponding strong forms even for linear elastic materials undergoing in finitesimal deformations.The well-known equilibrium equation based on the first Piola-Kirchhoff stress tensorPin Eq.(6)b can be obtained through the coordinate transformation as in Eq.(10).It can be easily linearilized in the undeformed con figuration.Its weak forms over the solution domain have been w idely used for developing numerical methods including the finite element methods.Without losing generality,the weak form of Eq.(6)b can be w ritten with a test functionw jas,

and for a continuous functionw j,

In order to make the strong form of the Cauchy stress tensor satisfied through the weak forms in Eq.(65)a,the trial functions of the first Piola-Kirchhoff stress tensorPneed to satisfy its strong form in Eq.(6)b.It also requires that the deformation gradient tensor satisfy the geometric identity[Shield(1967),Ogden(1975)],

as already discussed in Eqs.(6)(d-f).

If Eq.(66)is not satisfied in a numerical computation,various attendant numerical issues need special treatment,such as the well-known finite-element mesh compatibility issues of strong displacement continuity and strong traction reciprocity.For example,many “l(fā)ocking”issues in mesh-based numerical methods are related to the weak-form conditions,

Because of which,the extra con figurational forces are introduced.The crack problem is another example of a “strong singularity”at the crack-tip,as the value ofis in finite withinVε.

Hence,the Eshelby stress tensor as well as its alternate forms have been widely used in computing the con figurational forces in the scalar or vector forms.The attendant path-independent integrals have been computed using the contour integral method,the domain integral method[Nikishkov and Atluri(1987)],or the interaction integral method,as well as theT?integral for dynamic nonlinear problems[Atluri(1982),Nishoka and Atluri(1983)].Its vector form in Eq.(12)provides the directional strength of the singularities for crack initialization and propagation[Gurtin&Podio-Guidugli(1996),Kienler&Herrmann(2002)].It has also been used to evaluate and correct the incompatibility of the mesh-based trial functions,in order to avoid con figuration forces,such as in the selective integration scheme,in the assumed strain field methods,and in the use of high-order elements with proper terms etc.One of the recent applications has been to develop the locking-free mesh-based methods for Hamiltonian systems bychoosing only con figurational force-free termsunder high speed rotation[Garcia-Vallejo,Mikkola and Escalona(2007),Sugiyama,Gerstmaya,Shabana(2006),Zhao and Ren(2012)].

It is impossible to make the deformation gradient tensor satisfy Eq.(66)through mesh-based trial functions.The trial functions through the use of the meshless interpolations,such as the moving least squares approximations and the radial basis aroximations,have been widely studied.However,the high-order continuity of such meshless trial functions,through the global solution domain,does not imply that Eq.(66)is satisfied any better.It becomes even more computationally costly because a high-order numerical quadrature scheme is required if the global Galerkin approach is adopted.In contrast,the local meshless trial functions as in the Meshless Local Petrov Galerkin(MLPG)methods of Atluri et al(1998,2004),can satisfy Eq.(66)better,especially with the use of low order polynomial basis.Various test functions can be also chosen for computational efficiency,through the Meshless Local Petrov Galerkin(MLPG)approach[Atluri,et al(1998,2004)].It should be pointed out that the mixed MLPG method[Atluri,Han and Rajendran(2004)]becomes even more promising since the strain or stress can be interpolated independently along with the displacements as“the generalized degrees of freedom”.Hence,Eq.(66)can be satisfied in a better way.By mapping the deformation gradient or stress variables back to the nodal displacements,the“l(fā)ocking-free MLPG method”has been developed by[Atluri,Han and Rajendran(2004)].Since the balance laws for the Eshelby stress tensorTare essentially “weighted”forms of the momentum balance laws for the first Piola-Kirchhoff stressP,we may also make use of the Eshelby stress tensor to remove the restriction in Eq.(66),and the MLPG method can be extended to allow discontinuity in deformation.The inhomogeneity can also be included as “the distributed force on the defect”in Eq.(19).We call the resulting computational approach to solve for the displacements and the stress in a finitely deformed solid as the MLPG-Eshelby Method for computational solid mechanics.This represents a radical departure from the current state of computational solid mechanics.While the MLPG-Eshelby method for general computational finite deformation solid mechanics will be fully described in our forthcoming papers,a simple one-dimensional example is provided in the next section.

4　A simple example of the application of the Eshelby Stress Tensor in Computational Solid M echanics

4.1　Formulation

For illustration purposes,we use the identity for the Eshelby stress tensor,derived independently in Atluri(1982)[Eqs.(18)&(19)in Atluri(1982)].While computational methods for finite deformations of hyperelastic solids will be considered in our for thcoming papers,we consider now only a homogeneous linear elastic bar in a one-dimensional domain ?,with a boundary,subjected to the continuous body forceb(X)and the surface tractionand undergoing infinitesimal deformations.The bar is descretized intonsegments,each of an arbitrary length,as shown in Fig.2.

One may write the Eq.(18)in[Atluri(1982)]for one-dimensional elasto-static problems,as,

Figure 2:a bar divided into non-uniform segments

and the traction boundary conditions as

It is clearly seen that the combination of the first two terms in the left-hand side in Eq.(68)a leads to the definition of the Eshelby stress tensor in the present on edimensional problem[i.e.T=W?P·F].The corresponding strong form balance law for the Eshelby stress tensor can be w ritten as,

Letu(X)be the trial function for the displacements,over the solution domain,and let it be piece-wise continuousexcept at several points X=Zi.Eq.(69)is simply the weighted weak form of the equilibrium equation for the Piola-Kirchhoff stressPfor the one-dimensional problems,as

in which the constitutive relation forPin Eq.(3)is applied.One may choose the trial functionu(X)to be piece-wise continuous over each local sub-domain(exceptX1andXnhave only one side),but not necessarily continuous atZi=Xi?lifori=2,n.Hence,there aren?1 gaps between the local sub-domains,at the pointsZi.

One may chooseδX[a “variation”in the values of the coordinates of a material particle of the solid,in its initial con figuration]to be a test function,in order to write the weak-form of Eq.(70).We takeδXto be piece-wise continuous,except at several pointsX=Zi.The equilibrium equations for the Eshelby stress,as in Eq.(70),can be written in a weak form,over each local domain ?i,as

or by applying the divergence theorem,

forand

in which by definition,

From Eq.(72)it is clear thatδx(X)is not in the same class of functions as the trial function,u(X).Thus Eq.(71)for the Eshelby stress tensor necessarily implies a Petrov-Galerkin approach.

The summation of Eq.(71)b over the solution domain should be zero only if there are no gaps of the trial function[i.e.u(X)is continuous at every point].However,the trial function is not continuous but has gaps atX=Zi.The alternate test functionδx(X)also becomes non-continuous over the gapsX=Zidue to the discontinuity of the displacement gradients,F(X).If the actual solution contains no “real gaps”atX=Zi[i.e.δXis continuous],the corresponding gaps are introduced simply because of the discretization error,as well as the alternate test functionδx(X).However the weak form in Eq.(71)a is still valid over the domainZi?ε≤X≤Zi+ε.For illustration purposes,we may use this weak form over the gaps in the undeformed con figuration to preserve the balance of energy,rather than applying the balance law of the Eshelby stress tensor,to derive the corresponding weak form over the gaps.By ignoring the body force over the gaps and taking a small segmentcentered at the gapsX=Zi,one may take a linear continuous test function over the gaps as,

Thus one may have the weak form for the gaps as

By adding Eqs.(74)and(71)b,the weak form of the equilibrium equations can be obtained for the numerical discretization problem as:

On the other hand,if defects exist at the gapsX=Zi,the corresponding“forces on the defects”can be computed through Eq.(68).The test functionδXhas to be discontinuous atX=Zi.The average movement of the defects can be defined as

With the use of Eq.(27)or the weak forms in Section 3 forVε,the Eshelby tractiont?can be computed in terms of the strain energy changes ofVVεwhich also be used to drive the defects or break elements[LSTC(2013)].The applications of the Eshelby stress for discontinuous mechanics will be discussed in our following papers.

Eq.(75)may,at first glance,appear to be quite similar to the discretized equilibrium equations for the hybrid finite element methods[Atluri(1975)],which,however,are based on the direct“weak-form”of a momentum balance law forP,rather than on the balance law forP,weighted a priori withF.However,the concepts behind Eq.(75),and the hybrid finite element methods,are quite different.For a given set of nodes in the undeformed initial con figuration,{Xi},and a trial function for the displacementu(X),the corresponding stressP(X)can be computed following the standard finite element procedures.The“weak-form”of the momentum balance law forPstates thatu(X)is the solution,if the deformed domainx=X+u(X)renders the “total potential energy”stationary.By taking any admissible test functionδu(X),the corresponding stresscomputed from the neighboring solutionmust also satisfy the same“weak-form”of the momentum balance law forP.It is also clear,thatin the usual Galerkin finite element methods,the test function δu(X)[or the“variation”of u(X)]belongs tothe same class of functions as u(X).In contrast,in present case of writing a“weak form”of the balance law for the Eshelby stress tensor,namely for Eq.(70),one may keep the trial functionu(X)unchanged but change the coordinates in the initial con figuration,namelyXi,to admissible neighboring pointsYi=Xi+δXi,thus resulting in a new node set{Yi}.Keeping the same trial functionu(X),the new trial function for the nodeYican be computed through Eq.(72),as

or within a piece-wise continuously defined sub-domain,

It is clear that the test functionδx(X)used in Eq.(71)is dependent on both{Xi}andu(X),instead of being an admissible function of{Xi}only,as in the usual finite element methods.It is the “change of the trial function”caused by “mesh changes”.Within the general Galerkin approach,δx(X)is replaced with the displacement variation and thus the weak form in Eq.(75)is reduced to the momentum balance law forP.The conservation law for the energy of the system can not be preserved in the usual finite element method.In other words,the Galerkin approach enforces that the system’s energy conservation law is preserved only in the undeformed con figuration,instead of in any other con figurations.Noether’s theorem[Noether(1918)]states that the energy conservation law must be con figuration invariant.Hence,the test function in term ofδx(X)[as an alternate form ofδX]should be chosen differently from the trial function ofu(X)which essentially leads to the Petrov-Galerkin approach.

With a continuous test functionδX,the corresponding stressP(Y)can also be computed in a same way.The“weak-form”of the balance law forTstates thatu(X)is the solution if all“computed”solutionsu(Y)andP(Y)also satisfy the same weak form.If not,the new ly computed unbalanced nodal forces are the so-called con figurational forces.

If the mapping operation between the mesh changesδXand the changes of the trial functionδx(X)=FδXis invertible,through using a continuous trial function,any admissible test functionδx(X)can be used in Eq.(75),instead ofδXexplicitly.However,such compatible invertible relations can be not de fined for 2-D or 3-D problems for a global mesh-based or even a meshless interpolation.One may chooseδx(X)to be same asδu(X),which has been done in most global meshless Galerkin methods by ignoring the actual inverted functionδXthroughF.It may cause some numerical issues,such as i)the continuity requirement not allow ing any discontinuity or defects within the whole solution domain;ii)higher order quadrature schemes for high-order nonlinear integrands in the domain integrals requiring more material integration points which is computationally time costly for nonlinear materials;and iii)difficulty in enforcing the essential boundary conditions asδXis mapped to different class of funtions.On the other hand,a compatible and invertible relations betweenδx(X)andδXcan be easily constructed in a closed form within “a local spatial patch”[or a local sub domain].Such techniques have been widely used in the error estimation and mesh adaptivity with higher order accuracy.Hence the weak-forms within a local sub domain become more convenient through the Meshless Local Petrov-Galerkin(MLPG)approach[Atluri(1998,2004)].

4.2　Numerical implementation

There are many ways to choose the trial functionu(X),over the solution domain with various orders of continuity,in the initial coordinates.First the moving least squares(MLS)approximation is used to construct the trial function based on the fictitious nodal value[Atluri(2004)],as

The continuity of the trial functionuMLSis dependent on the weight functionsw(i)(X)in the MLS interpolation[Atluri(2004)].In the present study,we choose the fourth-order spline function as the weight function,which leads to a continuous trial function.

Secondly,the mixed interpolation for bothF(X)andu(X)can also be used to construct the trial functions based on the fictitious nodal valueandF(i)[Atluri,Han and Rajendran(2004)]over each segment,as

It is clear that the trial functionuMIX(X)is piece-wise linear,but discontinuous at pointsZi(i.e.with gaps).

The test functionδXdoes not explicitly appear in Eq.(75),instead onlyδxappears.Hence,one may de fineδxover the solution domain,and the correspondingδXmay be computed through Eq.(72)which is not involved in the numerical computation.It needs to be pointed out that Eq.(72)can not be defined globally(except in their global boundary/domain weighted integral forms)rather than within a local sub-domain.It is another reason which precludes the Eshelby stress tensor from being implemented through the global mesh-based methods[Such as the usual finite element methods,or the global Galerkin methods].Eq.(75)becomes nonlinear inδxif the test functionδXis assumed first,andδxthen is computed through Eq.(72).

The test functionδxcan be chosen to be piece-wise linear for a simple mapping relations betweenδXandδx.First the element-based simple polynomial shape functions are chosen to interpolate the test function based on the nodal valuesv(i)over each segment,as,

in which the shape functionsandare linear.The test functionδxFEMis piece-wise linear and possessesC0continuity(i.e.no gaps),and a linear relation betweenδXandδxcan be obtained within elements if the trial function is continuous,such asuMLS.

On the other hand,the mixed MLS interpolation is also chosen to construct the discontinuous trial functionuMIX.A simple continuous test functionδXcan be chosen as discussed in Section 2,as:

whereδX(i)andδλ(i)are two independent nodal variables.

Thus a linear relation betweenδXandδxcan be computed within each local subdomain,by definition,as

in which the nodal valuesδx(i)andδF(i)are independent variables and different from those used in Eq.(80).

The MLPG approach based on the Eshelby stress tensor[hereafter labeled as the MLPG-Eshelby Method]is presented here first by choosinguMLSas the trial function,andδxFEMas the test function.We call this the “Primal MLPG-Eshelby Method”.The domain integrals in Eq.(75)are performed over the solution domain without any singular gaps.

The second MLPG-Eshelby method is formulated by choosinguMIXas the trial function andδxMIXas the test functions,and labeled as the “Mixed MLPG Eshelby method”.The domain integrals in Eq.(75)can be simplified by evaluating the constant terms,as

It is interesting that no interpolation is involved in Eq.(84)other than evaluating the nodal values,including the displacement gradients and stresses.Essentially Eq.(84)becomes a “Particle”method which is computationally efficient.In addition,the second order term inδF(i)can be omitted and the system can be written as,

4.3　Numerical results

The bar is fixed at the left hand end and subjected to four loading conditions

i)uniform tension with zero body force(constant stress),as

b(X)=0 andP(Xn)=1

ii)constant body force as gravity load(linear stress),as

b(X)=1 andP(Xn)=0

iii)linear body force as centrifugal force(second order nonlinear stress),asb(X)=XandP(Xn)=0

iv)second order body force(third order nonlinear stress),as

b(X)=X2andP(Xn)=0

The bar is discretized regularly into 10 sub-domains,or irregularly with maximum 30%random variation from the regular sub-domains.All nodal coordinates are listed in Table 1.

The normalized relative displacement errors are shown in Figs.3-6 for the regular sub-domains,and in Figs.7-10 for the irregular sub-domains.The numericalresults show that the “Primal MLPG Eshelby Method”and the“M ixed MLPG Eshelby Method”pass the patch test and are quite stable.

Table 1:Node coordinates of a bar

Acknow ledgement:The second author thanks the US Government,and M r.Dy Le of ARL and M r.Nam Phan of NAVAIR,for their support.This work was completed during the second author’s pleasant and quiet stay at the Texas Institute for Advanced Study at the Texas A&M University;and thanks are expressed to John Junkins,Dimitris Lagoudas,and John Whitcomb for their hospitality and many illuminating conversations.

Atluri,S.N.(1975):On ‘Hybrid’Finite Element Models in Solid Mechanics,Advances in Computer Methods for Partial Differential Equations,R.Vichrevetsky,Ed.,AICA,pp.346-355.

Atluri,S.N.(1982):Path-Independent Integrals in Finite Elasticity and In elasticity with Body Forces,Inertia,and Arbitrary Crack-Face Conditions,Engineering Fracture Mechanics,vol.16,no.3,pp.341-364.

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.

Atluri,S.N.(2004):The Meshless Local Petrov-Galerkin(MLPG)Method for Domain&Boundary Discretizations,Tech Science Press,665 pages.

Atluri,S.N.;Han,Z.D.;Rajendran,A.M.(2004):A New Implementation of the Meshless Finite Volume Method,Through the MLPG “M ixed”Approach,CMES:Computer Modeling in Engineering&Sciences,vol.6,no.6,pp.491-514.

Atluri,S.N.;Zhu,T.(1998):A new meshless local Petrov-Galerkin(MLPG)approach in computational mechanics.Comput.Mech.,vol.22,pp.117-127.

Figure 3:Normalized relative displacement errors of a bar under uniform tension(regular sub-domains)

Figure 4:Normalized relative displacement errors of a bar under gravity load(regular sub-domains)

Figure 5:Normalized relative displacement errors of a bar under centrifugal force(regular sub-domains)

Figure 6:Normalized relative displacement errors of a bar under second order body force(regular sub-domains)

Figure 7:Normalized relative displacement errors of a bar under uniform tension(irregular sub-domains)

Figure 8:Normalized relative displacement errors of a bar under gravity load(irregular sub-domains)

Figure 9:Normalized relative displacement errors of a bar under centrifugal force(irregular sub-domains)

Figure 10:Normalized relative displacement errors of a bar under second order body force(irregular sub-domains)

Chadwick,P.(1975):Applications of an energy-momentum tensor in non-linear elastostatics,Journal of Elasticity,vol.5,nos.3-4,pp 249-258.

Eshelby,J.D.(1951):The Force on an Elastic Singularity,Phil.Trans.R.Soc.Lond.A,vol.244,pp.87-112.

Eshelby,J.D.(1957):The Determination of the Elastic Field of an Ellipsoidal Inclusion,and Related Problems,Phil.Trans.R.Soc.Lond.A,vol.241,pp.376-396.

Eshelby,J.D.(1959):The Elastic Field Outside an Ellipsoidal Inclusion,Phil.Trans.R.Soc.Lond.A,vol.252,pp.561-569.

Eshelby,J.D.(1975):The elastic energy-momentum tensor,Journal of Elasticity,vol.5,nos.3-4,pp.321-335.

Garcia-Vallejo,D.;M ikkola,A.M.;Escalona,J.L.(2007):A new locking-free shear deformable finite element based on absolute nodal coordinates,Nonlinear Dyn.,vol.50,pp.249-264.

Sugiyama,H.;Gerstmaya,J.;Shabana,A.A.(2006):Deformation modes in the fnite element absolute nodal coordinate formulation,Journal of Sound and Vibration,vol.298,pp.1129-1149.

Gurtin,M.E.;Podio-Guidugli,P.(1996):Con figurational forces and the basic laws for crack propagation,Journal of the Mechanics and Physics of Solids,vol.44,no.6,pp.905-927.

Gurtin,M.E.(2002):A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations,Journal of the Mechanics and Physics of Solids,vol.50,no.1,pp.5-32.

Han.Z.D.;Atluri,S.N.(2002):SGBEM(for Cracked Local Subdomain)–FEM(for uncracked global Structure)Alternating Method for Analyzing 3D Surface Cracks and Their Fatigue-Grow th.CMES:Computer Modeling in Engineering&Sciences,vol.3,no.6,pp 699-716.

Han,Z.D.;Atluri,S.N.(2003):On simple formulations of weakly-singular tBIE&dBIE,and Petrov-Galerkin approaches,CMES:Computer Modeling in Engineering&Sciences,vol.4,no.1,pp.5-20

Hill,R.(1986):Energy-Momentum Tensors in Elastostatics:Some re flections on the General Theory,J.Mech.Phys.Solids,vol.34,no.3,pp.305-317.

K ienler,R.;Herrmann,G.(2002):Fracture criteria based on local properties of the Eshelby tensor,Mechanics Research Communications,vol.29,pp.521-527.

Knowles,J.K.;Sternberg,E.(1972):On a Class of Conservation Laws in Lin-earized and Finite Elastostatics,Archive for Rational Mechanics and Analysis,vol.44,no.3,pp.187-211.

LSTC(2013):LS-DYNA Keyword User’s Manual,R7.0,Feb.2013,Vol.1.

Lubarda,V.A.;M arkenscoff,X.(2003):On conservation integrals in m icropolar elasticity,Philosophical Magazine,vol.83,no.11,pp.1365-1377.

Maugin,G.A.(1995):Material forces:concepts and applications,Appl.Mech.Rev.,vol.48,pp.213-245.

Mura,T.(1991):Micromechanics of defects in solids,2nded.,Dordrecht,Kluwer(1991).

Nikishkov,G.P.;Atluri,S.N.(1987):Calculation of Fracture Mechanics Parameters for an Arbitrary 3-D Crack by the EDI Method,International Journal for Numerical Methods and Engineering,vol.24,pp.1801-1822.

Nishioka,T.;Atluri,S.N.(1983):Path-Independent Integrals,Energy Release Rates,and General Solutions of Near-Tip Fields in M ixed-Mode Dynam ic Fracture Mechanics,Engineering Fracture Mechanics,vol.18,no.1,pp.1-22.

Noether,E.(1918):Invariante Variationsprobleme,G?ttinger Nachrichten,Mathematisch-physikalische Klasse,vol.2,p.235.(English translation by M.A.Tavel,Transport Theory and Statistical Physics,vol.1,p.183,1971)

Ogden,R.W.(1975):Dual reciprocal states in finite elasticity,Journal of Elasticity,vol.5,no.2,pp.149-153.

Okada,H.;Rajiyah,H.;Atluri,S.N.(1989)a:A Novel Displacement Gradient Boundary Element Method for Elastic Stress Analysis with High Accuracy,J.Applied Mech.,April 1989,pp.1-9.

Okada,H.;Rajiyah,H.;Atluri,S.N.(1989)b:Non-hyper-singular integral representations for velocity(displacement)gradients in elastic/plastic solids(small or finite deformations),Comp.Mech.,vol.4,pp.165-175.

Qian,Z.Y.;Han,Z.D.;Atluri,S.N.(2004):Directly derived non-hyper-singular boundary integral equations for acoustic problems,and their solution through Petrov-Calerkin schemes,CMES-Computer Modeling in Engineering&Sciences,vol.5,no.6,pp.541-562.

Shield,R.T.(1967):Inverse deformation results in finite elasticity,Zeitschrift für angewandte Mathematik und Physik ZAMP,vol.18,no.4,pp 490-500.

Zhao,Z.;Ren,G.(2012):A quaternion-based formulation of Euler–Bernoulli beam w ithout singularity,Nonlinear Dyn.,vol.67,pp.1825-1835.