Path independent integrals to identify localized plastic events in two dimensions
Mehdi Talamali, Viljo Petäjä, Damien Vandembroucq, Stéphane Roux
aa r X i v : . [ c ond - m a t . o t h e r] M a y Path independent integrals to identify localized plastic events in two dimensions
Mehdi Talamali , Viljo Pet¨aj¨a , Damien Vandembroucq ,
1: Unit´e Mixte CNRS/Saint-Gobain “Surface du Verre et Interfaces”39 Quai Lucien Lefranc,93303 Aubervilliers cedex, FRANCE2: Laboratoire PMMH, ESPCI/CNRS/Paris 6/Paris 710 rue Vauquelin, 75231 Paris cedex 05, France3: LMT-CachanENS de Cachan / CNRS-UMR 8535 / Universit´e Paris 6/PRES UniverSud61 avenue du Pr´esident Wilson,F-94235 Cachan Cedex, France.
We use a power expansion representation of plane elasticity complex potentials due to Kolossovand Muskhelishvili, to compute the elastic fields induced by a localized plastic deformation event.Far from its center, the dominant contributions correspond to first order singularities of quadrupolarand dipolar symmetry which can be associated respectively to pure deviatoric and pure volumetricplastic strain of an equivalent circular inclusion. Constructing holomorphic functions from thedisplacement field and its derivatives, it is possible to define path independent Cauchy integralswhich capture the amplitudes of these singularities. Analytical expressions and numerical testson simple finite element data are presented. The development of such numerical tools is of directinterest for the identification of local structural reorganizations which are believed to be the keymechanisms for plasticity of amorphous materials.
I. INTRODUCTION
Plasticity of amorphous materials has motivated an increasing amount of studies in recent years. In the absenceof underlying crystalline lattice in materials such as foams, suspensions or structural glasses, it is generally acceptedthat plastic deformation results from a succession of localized structural reorganizations[1, 2, 3, 4]. Such changes oflocal structure release part of the elastic strain to reach a more favorable conformation and induce long range elasticfields. The details of such local rearrangements and of the internal stress they induce obviously depend on the precisestructure of the material under study, and its local configuration . However, the important observation is that outsidethe zone of reorganization, a linear elastic behavior prevails . Therefore, elastic stresses can be decomposed onto amultipolar basis and independently of the material details, it is possible to extract singular, scale-free, dominant termswhich can be associated to a global pure deviatoric or pure volumetric local transformations of an equivalent circularinclusion. In particular, the elastic shear stress induced by a localized plastic shear exhibits a quadrupolar symmetry.This observation has motivated the development of statistical models of amorphous plasticity at mesoscopic scalebased upon the interaction of disorder and long range elastic interactions[5, 6, 7, 8]. In the same spirit statisticalmodels were also recently developed to describe the plasticity of poly-crystalline materials[9]. Several numericalstudies have been performed recently to identify these elementary localized plastic events in athermal or moleculardynamics simulations of model amorphous material under shear [10, 11].The question remains how to identify and analyze these transformation zones. In analogy with the path independentRice J-integral[12] developed to estimate the stress intensity factor associated to a crack tip stress singularity, we aimhere at capturing the stress singularity induced by the local plastic transformation which can be treated as anEshelby inclusion[13]. In two dimensions, we develop a simple approach based upon the Kolossov-Muskhelishvili(K&M) formalism of plane elasticity[14]. This is an appealing pathway to the solution since these zones will appearas poles for the potentials, and hence Cauchy integrals may easily lead to contour integral formulation which areindependent of the precise contour geometry, but rather rely on its topology with respect to the different poles whichare present.Although these techniques have been mostly used in the context of numerical simulations in order to estimate stressintensity factors from finite element simulations, they are now called for to estimate stress intensity factors fromexperimentally measured displacement field from e.g. digital image correlation techniques. In this case, interactionintegral techniques[15] or least square regression[16] techniques have been applied. Noise robust variants have alsobeen proposed[17]. These routes could also be followed in the present case.Though the present work is restricted to two dimensions due to the complex potential formulation, similar questionscan be addressed for the three dimensional version of this problem using the same strategy but a different methodology.In the following, we briefly recall the K&M formalism, and we give analytic expressions of contour integrals allowingto capture the singular elastic fields and we present a few numerical results based on a finite element simulationsupporting our analytical developments.In Section 2, we present the theoretical basis of our approach in terms of singular elastic fields, while Section3 introduce the contour integral formulation. In section 4, a numerical implementation based on finite elementsimulations is presented together with the results of the present approach. This application allows to evaluate theperformance and limitations of the contour integral procedure and check the detrimental effect of discreteness. Section5 presents the main conclusions of our study.
II. POTENTIAL FORMULATION
In two dimensions the Kolossov and Muskhelishvili potentials (K&M) can be used to write the elastic stress anddisplacement fields U and σ [14]. Using a complex formulation, we introduce the elastic displacement U = U x + iU y and the stress tensor field through two functions, the real trace S = σ xx + σ yy and the complex function S = σ yy − σ xx + 2 iσ xy . In the framework of linear elasticity, balance and compatibility equations can be rewritten as: S ,z − S ,z = 0 , (1) S ,zz = 0 , (2)where z = x + iy is the complex coordinate and the notation A ,x is used to represent the partial derivative of field A with respect to coordinate x . Note that we assumed zero surface density force and that equation (2) is here theclassical Beltrami equation which expresses the kinematic compatibility condition in terms of stress. The generalsolution to these equations can be obtained through the introduction of two holomorphic functions ϕ and ψ , calledthe K&M potentials. The displacement and the stress field can be written [14]2 µ U = κϕ ( z ) − zϕ ′ ( z ) − ψ ( z ) , (3) S = 2 h ϕ ′ ( z ) + ϕ ′ ( z ) i , (4) S = 2 [ zϕ ′′ ( z ) + ψ ′ ( z )] , (5)where µ is the elastic shear modulus and κ = (3 − ν ) for plane strain and κ = (3 − ν ) / (1 + ν ) for plane stress, ν being the Poisson’s ratio. III. PLASTIC INCLUSION AND SINGULARITY APPROACH IN 2DA. Singular terms associated to plastic inclusion
This K&M formalism can be applied to two-dimensional inclusion problems[18, 19]. Let us consider the case of asmall inclusion of area A experiencing plastic deformation and located at the origin of the coordinate system z = 0.It is assumed that the stress is a constant at infinity. Outside the inclusion, the K&M potentials can be expanded asa Laurent series as ϕ ( z ) = α out z + ∞ X n =1 ϕ n z n , ψ ( z ) = β out z + ∞ X n =1 ψ n z n (6)The linear terms can be easily identified as corresponding to uniform stresses while constant terms (omitted here)would lead to a rigid translation. It can be shown in addition that the dominant singular terms ϕ /z and ψ /z canbe associated to the elastic stress induced by the plastic deviatoric and volumetric of an equivalent circular inclusionof area A . Namely, considering a circular inclusion experiencing a plastic shear strain γ p and a plastic volumetricstrain δ p we have[18]: ϕ = 2 iµ A γ p π ( κ + 1) , ψ = − µ A δ p π ( κ + 1) . (7)In particular, for a pure shear plastic event we obtain a quadrupolar symmetry: σ xy = − γ p µκ + 1 A πr cos(4 θ ) (8)Note that we have in general to consider a complex value of γ p to include the angular dependence of the principalaxis. In contrast, the amplitude ψ is a real number (note that the imaginary part would correspond to a point-liketorque applied at the origin). B. Generic character of the expansion
Because of the well-known property of Eshelby circular inclusion, the above expansion limited to the ψ and φ terms only is the exact (outer) solution of a uniform plastic strain distributed in the inclusion, and vanishing stressat infinity. However, one should note that this result is much more general. Indeed, it is seen that the physical size ofthe inclusion does not enter into the solution but through the products A γ p and A δ p . Therefore, a smaller inclusionhaving a larger plastic strain may give rise to the very same field, provided the products remain constant. Therefore,one can consider the prolongation of the solution to a point-like inclusion (with a diverging plastic strain), as beingequivalent to the inclusion.Then from the superposition property of linear elasticity, a heterogeneous distribution of plastic strain γ p ( x ) in acompact domain, D , will give rise to such a singularity with an amplitude equal to[ A γ p ] eq = Z Z D γ p ( x ) dx (9)and the same property would hold separately for the volumetric part. As a particular case, one finds a uniform plasticstrain for an inclusion of arbitrary shape.This is the key property that allows to capture the equivalent plastic strain of an arbitrary complex configuration,for the above mentioned application to amorphous media. In fact, this is even the only proper way of defining theplastic strain for a discrete medium such as encountered in molecular dynamics simulations. The far-field behavior ofthe displacement and stress field can be accurately modelled, and without ambiguity by a continuum approach, andthus the above result will hold. In contrast, locally, the large scale displacement of several atoms may render difficultthe direct computation of the equivalent plastic strain experienced in such an elementary plastic event.Let us however stress one difficulty: As the above argument ignores the details of the action taking place withinthe “inclusion”, plasticity has to be postulated. However, a damaged inclusion, where the elastic moduli have beensoftened by some mechanism, or even a non-linear elastic inclusion at one level of loading would behave in a similarway as the above plastic inclusion. Obviously, to detect the most relevant physical description, one should haveadditional informations, say about unloading. If the above amplitudes remain constant during unloading, plasticitywould appear appropriate. If the amplitude decreases linearly with the loading, then damage is more suited. Finally,if the amplitudes varies reversibly with the loading, non-linear elasticity might be the best description. Thus, althoughone should be cautious in the interpretation, local damage detection from the far field may also be be tackled withthe same tools. C. Path independent contour integrals
In two-dimensions, this multipole expansion formalism in the complex plane suggests to resort to contour integralsto extract the singularities. However the displacement field is not a holomorphic function and cannot be used directlyfor that purpose. The strategy of identification of the singularities ϕ n and ψ n thus consists of expressing the potentialsfrom the displacement field and its derivatives in order to extract the singularities via Cauchy integrals. We now simplyexpress the displacement field and its derivatives:2 µU = κϕ ( z ) − zϕ ′ ( z ) − ψ ( z ) , (10)2 µU ,z = κϕ ′ ( z ) − ϕ ′ ( z ) , (11)2 µU ,z = − zϕ ′′ ( z ) − ψ ′ ( z ) , (12)2 µU ,zz = − ϕ ′′ ( z ) . (13)This gives immediately: ϕ ′ ( z ) = 2 µκ − (cid:2) κU ,z + U ,z (cid:3) , (14) ϕ ′′ ( z ) = − µU ,zz = − µ ∇ U , (15) ψ ′ ( z ) = − µ (cid:20) U ,z − z ∇ U (cid:21) . (16)Note that, except a multiplicative constant, the two last expressions are independent of materials properties. In lightof the expansion (6) of ϕ and ψ in Laurent series, if an anti-clockwise contour integration is considered along a path C , Cauchy residues can be formed as ϕ n = iµ πn ( n + 1) Z C z n +1 ∇ U dz , (17) ψ n = iµπn Z C z n (cid:20) U ,z − z ∇ U (cid:21) dz . (18)These expressions can thus be obtained from the sole knowledge of the displacement field, a quantity which can beaccessed from experiments, or from atomistic simulations. In the case of first order singularities (see Eq.(7)), theresidue term thus only depends on the local plastic deformation (size and amplitude of deformation) and on thePoisson’s ratio ν of the material.Reverting to Cartesian coordinates, where the contour is expressed as a function of the curvilinear abscissa s as( x ( s ) , y ( s )), the above expression can be written ϕ = µ π R C [ − xy ) + i ( x − y )][( U x,xx + U x,yy ) − i ( U y,xx + U y,yy )][ dxds + i dyds ] ds ,ψ = µ π R C [2[ ix − y ][( U x,x − U y,y ) − i ( U y,x + U x,y )] − [ i ( x + y )][ U x,xx + U x,yy − i ( U y,xx + U y,yy )] (cid:3) [ dxds + i dyds ] ds . (19) IV. NUMERICAL IMPLEMENTATION
The ultimate goal of a such a method would be to analyze numerical results obtained from molecular dynamicssimulations of amorphous plasticity where such local structural reorganizations are expected to take place. Thisobviously raises the question of a well defined method for writing the continuous displacement field from the data ofthe discrete displacements of particles [20] and more generally the question of the sensitivity to noise of the aboveexpressions. The first point is beyond the scope of the present work and we leave it for later studies. We thusfocus on the more restricted question of the numerical implementation and its efficiency in the case of artificiallynoise-corrupted displacement data.The method relies on contour integrations of derived fields of the displacement. The latter point induces a priori a strong sensitivity to noise. To limit such effects, first and second derivatives are extracted via an interpolation ofthe local displacement field by polynomial functions of the spatial coordinates. Moreover, the path independence ofthe contour integrals allows to perform spatial averages. We explore in the following the efficiency of this method onnoisy data.
A. Numerical generation of elastic fields induced by plastic inclusions
Displacement fields are computed numerically using a finite element code, with square elements and bi-linear shapefunctions { , x, y, xy } . Plane stress conditions of two-dimensional elasticity are used. The domain is a 150 ×
150 square.Stress free conditions are enforced all along the domain boundary. The Poisson’s ratio of the material is ν = 0 . A γ p and A δ p will differ slightly from the theoretical expectation.Nevertheless the path independence, (size, shape, center, ...) is expected to hold.We limited ourselves to such a description in order to mimic the difficulties one may face when having to deal withdiscrete element simulations. Indeed, the chosen finite element shape functions do not allow us to use this descriptionstrictly speaking in order to compute second order differential operators on the displacement, since the gradients ofthe latter are not continuous across element boundaries. Therefore, a regularization will be called for, as detailedbelow.The choice of a regular square lattice is obviously oversimplified compared with the case of the random latticesassociated with atomistic simulations. However, as x and y directions are obviously equivalent for square elementsand as linearity is preserved by the finite element formulation, this formulation should not introduce any breaking −0.4−0.3−0.2−0.100.10.20.30.4 AB C
FIG. 1: Map of the displacement field U y induced by a plastic shear strain of a central square element. The oriented pathsindicate contours for integration. R -0,2500,250,50,7511,25 ϕ , ψ Re φ square contourIm φ square contourRe ψ square contourIm ψ square contourRe φ cross contourIm φ cross contourRe ψ cross contourIm ψ cross contourRe ϕ rect. contourIm ϕ rect. contourRe ψ rect. contourIm ψ rect contour R -0,2500,250,50,7511,25 ϕ , ψ Re φ square contourIm φ square contourRe ψ square contourIm ψ square contourRe φ cross contourIm φ cross contourRe ψ cross contourIm ψ cross contourRe ϕ rect. contourIm ϕ rect. contourRe ψ rect. contourIm ψ rect. contour FIG. 2: Normalized values of numerical estimates ϕ and ψ obtained for three families of contours, respectively square, crossand rectangle shaped and of varying length L in the case of a displacement field induced by the plastic shear strain (left) orcontraction (right) experienced by the central element of a square lattice. Theoretical expectations are ϕ = i , ψ = 0 (left)and ϕ = 0 , ψ = 1 (right). of symmetry. More specifically, the displacement field induced by a quadrupole of principal direction off axis can beobtained by a linear superposition of x and y components of the displacement field induced by a quadrupole alignedwith the axis weighted by the sine and cosine of the quadrupole angle.Finally, the finite size of the system is also a specific difficulty encountered in practice, whereas the above argumentuses the assumption of an infinite domain. However, such a boundary condition should not induce additional poleswithin the domain, and can be considered as a common practical difficulty encountered for all practical use of thistool. All these arguments are possible causes of deviation from the theoretical expectation of path independence, andit thus motivates a detailed study of the method stability, robustness and accuracy.Two test cases are studied. I: a central inclusion experiencing a shear strain γ p = 1 (because of linearity, the actualamplitude is meaningless) along the x − axis; II: a central inclusion experiencing a volumetric contraction δ p = −
1. Amap of the displacement fields U y in case I is given on Fig. 1. -30 -20 -10 0 10 20 30 x -0.4-0.200.20.40.60.811.21.4 φ , ψ Re φ Im φ Re ψ Im ψ Re ψ theor. FIG. 3: Normalized values of numerical estimates ϕ obtained for square contours of center ( x ,
0) and of size M = 20. Theexpected behavior of ℜ ψ (unity when the inclusion lies within the contour, zero otherwise) is represented by the bold line.Other quantities are expected to be zero. B. Interpolating displacement data
The key ingredient is to go from a continuous but non-differentiable displacement field obtained from the finiteelement simulation to an evaluation of the second derivative at any point in the domain.The results which are presented below have been obtained using the following procedure. A quadratic fit is performedon a square centered on one node to extract the first and second order derivatives, the obtained values are used tocompute the integrals by quadrature. An alternative method has been tested: for an integration from ( x, y ) to( x + 1 , y ), a simple fit is performed of the 12 nodes ranging from ( x −
1) to ( x + 2), and from ( y −
1) to ( y + 1), by thetensor product of polynomials (1 , x, x , x ) and (1 , y, y ) (12 functions). Then the integral of all required quantitiescan be computed. Estimates of derivatives using of Fourier Series with and without low pass filtering have also beenperformed. All methods give similar results provided that the area of the region used for interpolation (or filtering)is comparable. C. Path independence
We first check the path independence of the contour integral in the cases I and II of isolated inclusions. For thatpurpose, we use two families of contours, square (A) and cross (B) and rectangular (C) shaped respectively as shownon Fig. 1. The size of these contours as well as their center can be varied. Figure 2 gives a summary of the results.For the three kinds of contours, we show the real and imaginary parts of the residues corresponding to Eq. (17). Notethat the numerical results have been normalized according to the theoretical expectations (7) so that the expectednumerical values are ϕ = i , ψ = 0 in case I (Fig. 2 left) and ϕ = 0 , ψ = 1 in case II (Fig. 2 right).These numerical results can be considered as rather satisfactory in terms of orientation and orthogonality betweenmodes ϕ and ψ : the measured values of quantities whose expected value is zero remain typically below 10 − . Whencompared to their theoretical values, ϕ s hear and ψ c ontraction exhibit relative differences of around 5-10% . Smallfluctuations (below 5%) can be found when changing shape and size of the contours. We already commented on thefact that the finite element simulation are performed with a single element for the inclusion, a procedure which isobviously not reliable in terms of accuracy, but which allows us to have a large ratio between element and system size.Another test of the dependence of the numerical procedure is on the sole topology, i.e. location of the inclusioninside or outside the contour, we show in addition the dependence of the measured values of ϕ and ψ on the locationof the contour center. Fig. 3 shows the singularity ψ measured from the integration of displacement field II along asquare contour centered along the x -axis. Results are normalized so that the expected value of ℜ ψ be unity whenthe inclusion is within the contour and zero elsewhere.The contour size is M = 20. We obtain the expected behavior:the values of ℜ ψ shifts form zero to unity depending on the inclusion is within or outside the contour. Significantfluctuations (10-20%) are however observed when the inclusion lies in the vicinity of the contour.Finally we test the dependence of the numerical method on the material properties. In the determination of ϕ and ψ (17) as residues, the need to resort to a second order derivative of the displacement field is balanced by the ν φ ψ theor. expectation FIG. 4: Numerical estimates of | ϕ | and ψ obtained for a central inclusion experiencing a unit shear and a unit contractionrespectively as a function of the Poisson ratio ν . The numerical results obtained in plane stress conditions for a system of size100 ×
100 with stress free boundary conditions are compared with the theoretical expectation (1 + ν ) / π . fact that the computation can be performed without any knowledge of the elastic properties of the material. Theindependence of the numerical procedure with the Young modulus is trivially obtained due to the linearity of the FEMcomputation. In Fig. 4 we show the dependence of the numerical results on the Poisson ratio. FEM computations havebeen performed on systems of size 100 ×
100 with stress free boundary conditions and a central inclusion experiencinga unit shear and a unit contraction respectively. Poisson ratios have been varied from 0.05 to 0.45 by step of 0.05.The results shown in the figure compare the numerical estimates obtained for a square contour of size 40 centered onthe inclusion with the theoretical expectation ψ = ϕ = (1 + ν ) / π . The numerical results show that the volumetricstrain is weakly dependent on Poisson’s ratio, but the elementary shear is more poorly estimated for high Poisson’sratio. V. DISCUSSION/CONCLUSION
The proposed approach is based on a an exact result and hence theoretically establish a parallel with other typesof elastic singularities (in particular for stress intensity factors which characterize crack loadings) where similar pathintegral are well known. When tested on direct numerical simulations, we could recover the main topological propertiesexpected in this context: path independence and detection of the absence/presence of a singularity within the contour.However, the quantitative results proved more disappointing: the method is rather unprecise on the determinationof the prefactor of the singularity and is more generally rather sensitive to noise. The main cause is presumablydue to the inconsistent regularity of the displacement field solution (simple continuity) with the need to resort toestimates of first and second order differentials. A piecewise high order polynomial interpolation is operational forintegrals over finite segments, however, from one segment to the next, first and second order differentials will displaya discontinuous character, which obviously affect the method and result. Moreover, being a path integral, the methoddoes not take advantage of the knowledge of the displacement field at all points of a domain. To make the methodmore robust with respect to noise, different approaches can be followed. One natural way is to average the result overdifferent contours, thus transforming the contour integral into a domain integral. An arbitrary weight average canalso be considered, and hence, one could optimize the weight in order to achieve the least noise sensitivity. Such amethod was explored successfully for cracks in ref. [17]. Note finally that extensions to three dimensions obviouslyrequire a different technique than Kolossov and Muskhelishvili potentials, and contour integrals, however still a linearextraction operator acting on the displacement field can be computed to provide similarly the equivalent plastic strain.
Acknowledgments
V. Pet¨aj¨a acknowledges the financial support of ANR program PlastiGlass NT05-4 41640 and of the academy ofFinland. [1] V. V. Bulatov and A. S. Argon, Modell. Simul. Mater. Sci. Eng. , 167 (1994).[2] V. V. Bulatov and A. S. Argon, Modell. Simul. Mater. Sci. Eng. , 185 (1994).[3] V. V. Bulatov and A. S. Argon, Modell. Simul. Mater. Sci. Eng. , 203 (1994).[4] M. L. Falk and J. S. Langer, Phys. Rev. E , 7192 (1998).[5] J.-C. Baret, D. Vandembroucq, and S. Roux, Phys. Rev. Lett. , 195506 (2002).[6] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, Phys. Rev. E , 010501(R) (2005).[7] A. Lemaˆıtre and C. Caroli, arxiv:cond-mat/0609689v1 (2006).[8] E. A. Jagla, Phys. Rev. E , 046119 (2007).[9] M. Zaiser and P. Moretti, J. Stat. Mech. P08004 , 79 (2005).[10] A. Tanguy, F. Leonforte, and J.-L. Barrat, Eur. Phys. J. E , 355 (2006).[11] C. E. Maloney and A. Lemaˆıtre, Phys. Rev. E , 016118 (2006).[12] J. Rice, J. Appl. Mech. , 379 (1968).[13] J. D. Eshelby, Proc. Roy. Soc. A , 376 (1957).[14] N. Muskhelishvili, Some basic problems of the mathematical theory of elasticity (Groningen, 1953).[15] J. R´ethor´e, A. Gravouil, F. Morestin, and A. Combescure, Int. J. Frac. , 65 (2005).[16] S. Roux and F. Hild, Int. J. Frac. , 141 (2006).[17] J. R´ethor´e, S. Roux, and F. Hild, Eng. Frac. Mech. in press (2008).[18] M. Jawson and R. Bhargava, Proc. Cam. Phil. Soc. , 669 (1961).[19] J. Mathiesen, I. Procaccia, and I. Regev, arxiv:0704.1206v1 (2007).[20] C. Goldenberg and I. Goldhirsch, Eur. Phys. J. E9