- Email: [email protected]
Viscous inclusions in anisotropic materials: Theoretical development and perspective applications
Viscous inclusions in anisotropic materials: Theoretical development and perspective applications Dazhi Jiang PII: DOI: Reference:
S0040-1951(16)30458-9 doi: 10.1016/j.tecto.2016.10.012 TECTO 127279
To appear in:
Tectonophysics
Received date: Revised date: Accepted date:
16 March 2016 7 October 2016 13 October 2016
Please cite this article as: Jiang, Dazhi, Viscous inclusions in anisotropic materials: Theoretical development and perspective applications, Tectonophysics (2016), doi: 10.1016/j.tecto.2016.10.012
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Viscous inclusions in anisotropic materials: theoretical development and perspective applications Dazhi Jiang1,2 Department of Earth Sciences, Western University, London, Ontario, N6A 5B7, Canada
PT
1
2
RI
([email protected])
State Key Laboratory of Continental Dynamics, Department of Geology, Northwest University,
SC
Xi’an, 710069, China
AC CE P
TE
D
MA
NU
Abstract: Theories and numerical solutions for a viscous ellipsoid in an infinite anisotropic viscous medium subjected to far-field homogeneous deformation lie at the heart of selfconsistent homogenization models and multiscale simulations of texture and fabric development in Earth’s lithosphere. There is considerable literature on ellipsoid inclusions, focused on anisotropic elastic materials, published in multi-disciplinary fields. To make this body of work more accessible as well as to advance viscous inclusion studies, an effort is made here to summarize recent advances and to further develop formally more explicit and, where possible, analytic solutions for incompressible viscous materials. The point-force concept and equivalent inclusion method of Eshelby are used together with the powerful Green function approach. This leads to generalized equations for ellipsoid inclusion behaviors in anisotropic materials. In the particular case of isotropic materials, the new mathematical development here enables the use of existing methods for elastic materials to get solutions for corresponding viscous inclusion problems efficiently and accurately. A 2D formulation is also presented for elliptic cylinders in plane-straining flows of anisotropic materials, using the same Green function method that is adopted for 3D inclusions. The 2D formulation is benchmarked with existing analytic solutions. A reconnaissance investigation to compare the behaviors of 2D elliptic inclusions and triaxial ellipsoids in a matrix of planar anisotropy undergoing far-field plane-straining flows suggests that conclusions drawn from 2D cannot be applied to 3D in anisotropic cases. The application of the viscous inclusion theory to the rheologically heterogeneous and non-linear lithosphere is discussed. By regarding a heterogeneous element as an ellipsoidal inclusion embedded in a hypothetical homogeneous equivalent matrix whose effective rheology is obtained through micromechanical homogenization and assuming that the conditions of scale separation are satisfied, the viscous inclusion theory forms the backbone of multiscale models for natural deformation and the accompanying texture and fabric development. Keywords: fabrics; viscous materials; multiscale modeling; inclusion problem; Eshelby solution; anisotropic viscosity; homogenization; rheology; flow partitioning
1
ACCEPTED MANUSCRIPT 1. Introduction
AC CE P
TE
D
MA
NU
SC
RI
PT
The instantaneous mechanical state and progressive evolution of a solitary inclusion in an infinite uniform material subjected to far-field homogeneous deformation, referred to as the “inclusion problems” below, have been studied extensively for a century and remain an active research subject today (see reviews by Mura, 1987, 1988; Ting, 1996; Zhou et al., 2013). Inclusion problems are significant in a wide variety of fields both theoretically and in practical applications. Micromechanics (e.g., Mura, 1987; Nemat-Nasser and Hori, 1993), the most recent development in applying continuum mechanics theories to real materials (Qu and Cherkaoui, 2006, p.2), stems from the seminal work of Eshelby (1957) on elastic inclusions. Earth scientists have applied and contributed to the advance of inclusion theories in the study of fabric development in the crust and strain and mixing of lithospheric heterogeneities in the mantle for half a century (e.g., Gay, 1968; Ramsay 1967, p.225; Hobbs et al., 1976, p.285; Passchier and Simpson, 1986; Passchier 1987; Spence et al., 1988; Spence and Wilmott, 1988; ten Brink and Passchier, 1995; Piazolo et al., 2002; Piazolo and Passchier, 2002; ten Grotenhuis et al., 2002; Treagus and Treagus, 2002; Ceriani et al., 2003; Mandal et al., 2003, 2005; Marques and Bose, 2004; Johnson et al., 2009; Li and Jiang, 2011; Mancktelow 2011, 2013; Xiang and Jiang, 2013; Chen et al., 2014). Among the considerable literature of inclusion studies, perhaps most influential to earth scientists are the work of Jeffery (1922) on the motion of a rigid ellipsoid in incompressible isotropic Newtonian fluids, its generalizations (Bretherton, 1962; Willis, 1977; Ferguson, 1979), and Eshelby’s (1957) work on isotropic elastic inclusions. The influence of the latter is largely through the works of Bilby et al. (1975) and Bilby and Kolbuszewski (1977) who applied Eshelby’s (1957) theory to isotropic incompressible viscous fluids. In the past two decades, Eshelby’s inclusion formalism has been extended to anisotropic viscous inclusions in nonlinear materials (Molinari et al., 1987; Lebensohn and Tomé, 1993). This has formed the backbone for the self-consistent viscoplastic model for texture development in poly-crystalline materials (Lebensohn and Tomé, 1993), the self-consistent MultiOrder Power Law Approach for multiscale fabric development in the heterogeneous lithosphere (Jiang and Bentley, 2012; Jiang 2014), and many self-consistent models of homogenization to obtain the bulk rheological properties from constituent properties (e.g., Molinari, 2002). Most works on inclusion problems in the past 5 decades are focused on elliptic inclusions in plane-strain (or plane-stress) and generalized plane deformation of anisotropic elastic materials (Jaswon and Bhargava 1961; Willis, 1964; Bhargava and Radhakrishna, 1964; Chen 1967; Yang and Chou, 1976; Hwa and Ting, 1989; Ting, 1996). These works adopt the complex variable formalism for plane problems of elasticity (Lekhnitskii 1950; Mushkelishvili 1953; Stroh 1958, 1962; Green and Zerna 1968). As will be shown in this paper, the work of Jaswon and Bhargave (1961) on elliptic inclusions in plane-strain (or plane-stress) deformation of isotropic elastic materials is directly applicable to corresponding 2D inclusions in isotropic viscous materials. However, it remains open how the works on 2D anisotropic elastic inclusions may be applied to corresponding viscous materials. The difficulty lies in the requirement of 2
ACCEPTED MANUSCRIPT
PT
considering incompressibility for viscous materials, which breaks the equivalence between linear elastic and incompressible Newtonian viscous problems. Molinari et al. (1987) first formally solve the incompressible 3D viscous inclusion problems in anisotropic viscous materials. Lebensohn et al. (1998) further develop this approach. Fletcher (2004, 2009) obtains analytic solutions for an elliptic isotropic viscous inclusion embedded in a viscous matrix with planar anisotropy, addressing incompressibility explicitly.
NU
SC
RI
Major progresses have also been made in developing numerical algorithms and programs to solve the equations of Jeffery (1922), Eshelby (1957, 1959), and more general ones for anisotropic materials (e.g., Freeman 1985, 1987; Ježek 1994; Jiang 2007a, 2007b, 2012; Meng, et al., 2012; Qu et al., 2016). Complete algorithms and software packages are now available for simulating single and multi-inclusion systems in general anisotropic viscous materials including non-linear power-law materials (Jiang 2012, 2013, 2014; Qu et al., 2016).
AC CE P
TE
D
MA
It is a bit challenging for one to follow the development on inclusion studies as the papers are published across a broad range of disciplinary fields and formulated in heavy and variable mathematical formalisms from one paper to another. Here I attempt to present a treatise on the subject. I review recent advances, develop more explicit formal and, where possible, analytical solutions for viscous inclusions in anisotropic materials in the narrative, and discuss the applications of the theoretical developments. The Green function method is used throughout the paper to render a more consistent mathematic notation. Many new equations governing the behavior of inclusions in anisotropic materials are established. New mathematical relations discovered in this paper enable efficient and accurate numerical solutions to viscous inclusion problems using existing equations for elastic inclusions. It is demonstrated that all existing results for ellipsoidal and elliptic inclusions can be derived readily from the general formulation as subsets or limiting cases. The more explicit formulation here greatly simplifies the efforts of numerical implementation of the inclusion theory in mathematic applications such as MATLAB, MATHCAD, and MAPLE.
2. Non-equivalence between compressible linear elasticity and incompressible linear viscosity It has been long held that linear elastic theories and linear viscous theories are formally equivalent (e.g., Spencer 1980, p.110-118). This implies that the former can be used, en masse, for the latter. However, elastic theories generally consider compressible solids whereas viscous materials are commonly assumed to be incompressible. This asymmetry does not cause undue difficulties for isotropic materials because their elastic stiffnesses are reduced to two scalar quantities (e.g., the bulk and shear moduli) which occur explicitly in elastic equations. One can address incompressibility easily by considering the limit case as the bulk modulus tends to infinity (or, equivalently, the Poisson’s ratio to 0.5). This was done by Bilby et al. (1975) and 3
ACCEPTED MANUSCRIPT
SC
RI
PT
Bilby and Kolbuszewski (1977). Later in this paper, it will be shown that Jaswon and Bhargava’s (1961) work on an elliptic inclusion in a 2D isotropic elastic medium can also be extended to viscous materials in a similar way. However, general anisotropic elastic materials have many non-vanishing stiffness terms and corresponding elastic equations are also more complex, generally expressed in convoluted and implicit forms. In 2D theories on anisotropic elastic solids, equations are also often expressed in heavy complex variable formalisms. This makes it non-trivial to consider incompressibility as a limiting case. It becomes more straightforward to regard incompressibility as an additional kinematic condition and solve the differential equation set for inclusion problems anew. The fundamental equations for linear anisotropic elasticity are listed below:
1 ui, j u j ,i 2
NU
ij
(1a)
el el ij Cijkl kl , ij M ijkl kl
constitutive law (generalized Hooke’s law)
(1b)
el Cijkl uk ,lj fi 0
equilibrium equations
(1c) (1d)
TE
el el el el s Cijpq M pqkl M ijpq C pqkl J ijkl
D
MA
compatible strain field
where ui , ij , ij , f i are respectively the displacement, elastic strain tensor, Cauchy stress
AC CE P
el el tensor, and body force, all referred to a fixed coordinate system xi i 1, 2, 3 ; Cijkl and M ijkl
are, respectively, 4th-order elastic stiffness and compliance tensors which are inverse of each s other in the sense of (1d), and J ijkl is the 4th-order symmetric identity tensor defined as
1 ik jl jkil ( ij being the Kronecker delta). The notation used in (1) and throughout 2 this paper follows Jiang (2014). A combination of Cartesian tensor notation and a standard tensor/matrix notation are used. Repeated indices in (1) and throughout the paper imply summation unless declared otherwise. A comma stands for differentiation with respect to the u Eulerian coordinate (e.g., ui , j i ). x j s J ijkl
To many readers, Eqs.1 are familiar. Eq.1a ensures that the elastic strain field is compatible. The generalized Hooke’s law (1b) relates the Cauchy stress tensor (symmetric, having 6 independent components) and the elastic strain tensor (also symmetric and having 6 el el independent components). Both Cijkl and M ijkl can be converted into 6 6 non-singular matrices (not the Voigt matrix notation) using a set of orthonormal second-order tensors so that their inversion (1d) calculations can be performed readily (see Jiang 2013 for more details).
4
ACCEPTED MANUSCRIPT In comparison, the corresponding fundamental equations for linear anisotropic viscosity are:
1 ui, j u j ,i 2
kk uk ,k 0
(2a)
PT
compatible strain rate field
RI
incompressibility
v v ijd Cijkl kl , ij M ijkl kld
rheology equations
(2c) (2d)
Cauchy stress field
v Cijkl uk ,lj p, j ij fi 0
equilibrium equations
NU
ij pij ijd
(2e) (2f)
MA
v v v v d Cijpq M pqkl M ijpq C pqkl J ijkl
(2b)
SC
ij
where, for easy comparison with (1), ui , ij , ij now stand for, respectively, the velocity,
D
viscous strain rate tensor, and the Cauchy stress tensor; ijd in Eq.2c stands for the deviatoric
AC CE P
TE
stress tensor; p is pressure. Eq.2b states the incompressibility condition. The rheology equations (2c) relate the deviatoric stresses (having 5 independent components) to the viscous strain rates v v (also having 5 independent components). The viscous stiffness ( Cijkl ) and compliance ( M ijkl ) are
1 d s J ijkl ij kl . inverse of each other in the deviatoric subspace (2f) where J ijkl 3 Note the viscous equilibrium equation (2e) differs from the elastic equilibrium equation (1c). Solving (2e) with the condition of (2b) for given boundary or initial conditions (a viscous problem) is different from solving (1c) for given boundary conditions (an elastic problem). This is why an elastic solution, in general, cannot be used directly for a corresponding viscous problem. In principle, one may still regard viscous incompressibility as a limiting case using a numerical approximation called ‘the penalty approach’ (e.g., Hutchinson 1976; Bilby et al. 1975; Lebensohn and Tomé 1993) which consists in using a modus operando ‘bulk viscosity’ ( b in the following expressions) to make an incompressible viscous material slightly compressible. Incompressible condition is approximated then by making b very large compared to other viscosity components. However, this method has serious drawbacks as explained below. A general penalty approach formulation is presented as follows. One can add a ‘slight compressibility term’ to express the viscous stiffness and compliance tensors as:
5
ACCEPTED MANUSCRIPT p p m v , M ijkl Cijkl 3b J ijkl Cijkl
1 m v J ijkl M ijkl 3b
PT
1 m v ij kl , and the superscript ‘p’ stands for the penalty approach. Unlike Cijkl where J ijkl and 3 v p p , the two tensors Cijkl and M ijkl are invertible in the Cauchy space (i.e., M ijkl
RI
p p p p s p p ). With Cijkl and M ijkl , equations for a viscous problem become Cijpq M pqkl M ijpq C pqkl J ijkl
SC
identical to those for a corresponding elastic problem and elastic solutions can be used. It would appear simple to get solutions for incompressible viscous problems by letting b tend to infinity.
NU
However, the difficulty lies in that fact that, except for isotropic materials, it is not possible to simply let b to get exact limit form expressions because of the implicit and complex nature of the solutions. This means that one has to use a large b numerically to obtain
MA
approximate solutions. Clearly a larger b would result in a better approximation of incompressibility. Unfortunately, because of the finite precision of the computer, too large a b v relative to the norm of Cijkl will deteriorate the numerical calculations due to ill conditioning of
TE
D
p and the inevitable accumulation of round-off errors (e.g., Pelletier et al., 1989). Cijkl
AC CE P
Because of this difficulty with the penalty approach, I shall consider incompressibility explicitly for general anisotropic viscous materials, advancing the work of Molinari et al. (1987) and Lebensohn et al (1998). This requires solving (2e) with condition (2b) afresh, which is the effort in the following section, through the powerful Green function method.
3. General solution for a single viscous inclusion in an infinite viscous matrix Two keystone ingredients in Eshelby’s (1957, 1959, 1961) [for a collection of Eshelby’s works, the reader is referred to Markenscoff and Gupta (2006)] treatment of the elastic inclusion problems are 1) the point force concept that allows an inclusion with eigenstrain to be treated as a distribution of body force and 2) the equivalent inclusion method that regards an inhomogeneous inclusion as a homogeneous one with proper eigenstrain. For a more detailed explanation of these two ingredients, the reader may refer to Jiang (2013, 2014). These ingredients are applied to solve the general three-dimensional viscous inclusion problems here, considering incompressibility explicitly. In the following, “viscous materials” shall always mean general “incompressible anisotropic viscous materials”. 3.1 Green functions for viscous materials The point force concept leads Eshelby (1957, 1959) to the use of elastic Green function for the solution of a general elastic inclusion problem. In linear elasticity, a single Green function 6
ACCEPTED MANUSCRIPT for displacement is used (e.g., Eshelby 1957). To address incompressibility explicitly, Molinari et al. (1987) and Lebensohn et al. (1998) use two separate Green functions: Gij x, x ' for velocity, and H i x, x ' for pressure. They relate, respectively, the velocity ui x and pressure
PT
p x at location x to the action of a unit point force at the source location x ' . With the Green
functions known, the velocity and pressure fields are then expressed by the convolution relations
RI
ui x Gij x, x ' f j x ' dV x ' and p x H i x, x ' fi x ' dV x ' , where f x ' is the body
kl* x ' x ' j
(e.g., Walker, 1993; Li and Wang, 2008, p.89;
NU
body force distribution fi x ' Cijkl
SC
force distribution in the deforming body and the volume containing all body force sources. The presence of an eigenstrain field kl* x ' in the medium is replaced by a layer of fictitious
where Cijkl is the elastic or viscous stiffness tensor of the material). With the equivalent
D
MA
inclusion method, a heterogeneous inclusion is then regarded as an equivalent homogeneous one with the appropriate eigenstrain (Eshelby 1957, see Jiang 2013, 2014 for more explanation). Therefore, the first important step toward solving the viscous inclusion problem is to find out expressions for the Green functions.
TE
Because in an infinite medium Gij x, x ' Gij x x ' and Hi x, x ' Hi x x ' , the Green functions are written as Gij x and H i x for simplicity below, setting the source point at the
AC CE P
origin. The Green functions for an infinite viscous material are solutions to the following auxiliary set of linear partial differential equations (Molinari et al. 1987, Eq.27 there): CijklGkm, jl x H m,i x im x 0
(3a)
Gkm,k x 0
(3b)
with boundary conditions
Gkm 0 , H m 0 at infinity
(3c)
where Cijkl is the viscous stiffness tensor of the material and x is the Dirac delta function. Applying Fourier transform to (3) converts them to the following algebraic system:
Aik K 2 gkm zi iKhm im
(4a)
zk K 2 gkm 0
(4b)
where 7
ACCEPTED MANUSCRIPT Aik Cijkl z j zl
(5a)
gkm (K )
G x exp iK x dx
(5b)
km
PT
hm (K )
H x exp iK x dx m
RI
(5c)
SC
The symmetric Aik is called the Christoffel stiffness tensor (e.g., Barnett, 1972; Walker, 1993),
gkm (K ) and hm (K ) are the Green functions in Fourier space, K Fourier wave vector, z unit
K ), and i 1 . The Green functions in Fourier space, gij and K hi respectively, can be solved from (4): i Aij ( z) , hi (K ) i 2 K K
MA
gij (K )
NU
vector parallel to K (i.e., z
(6)
ζ A T z
z 0
1
(7)
AC CE P
A T ζ
TE
D
where A , ζ , and , all functions of z , are submatrices of the following matrix:
We need Green functions in real space which are inverse Fourier transforms of (6). Analytical expressions for Gij x and H i x are known for isotropic materials (e.g., Bilby et al., 1975; Mura, 1987, p.22). Incompressible Green functions and their derivatives are unknown for anisotropic materials. In Appendix A, expressions of Gij , Gij ,l , and H i for viscous materials are obtained for the first time (A17) following the method of derivation of Barnett (1972) and Mura (1987) for elastic materials. The results are summarized below: Gij x
1
4 2 r 0
Aij d
(8a)
1 H i x 2 2 Fd i 4 r 0
(8b)
1 Gij ,l x 2 2 xl Aij zl M ij d 4 r 0
(8c)
8
ACCEPTED MANUSCRIPT where
r x , xi xi , Fi Ain mn m Aim i m xm , r
PT
M ij Aim A jn mn Aim j A jm i x m
RI
with C p q z p x q zq x p .
SC
The quantities Aij , Fi , and M ij in (8) are all functions of z which in the x z 0 plane
NU
can be expressed in terms of two mutually orthogonal unit vectors α and β (Synge, 1957; Barnett, 1972) shown in stereographic projection in Fig.1. The integrations of (8) are line integrals along a unit circle in the x z 0 plane.
Aij
1
ij
zi z j , Fi x i , and M ij
MA
For an isotropic viscous material with viscosity , it can be shown (Appendix A) that
1
z x i
j
z j x i . Explicit integration of (8) gives the
1 8 r
ij
xi x j
AC CE P
Gijiso x
TE
D
following well-known results (e.g., Bilby et al., 1975; Buryachenco, 2007, p.53, 56; Yin et al., 2014):
Gijiso,l x
H iiso x
1
8 r 2
x ij
l
jl
(9a)
x i il x j 3x i x j x l
xi 4 r 2
(9b)
(9c)
where the superscript “iso” stands for isotropic viscous materials. More explicit forms than (8) for Green functions and their derivatives are still not available for viscous materials. Some close-form expressions for elastic Gij and its derivatives have been found for transverse isotropic, orthotropic, and cubic materials, and on the symmetry plane of monoclinic materials (e.g., Ting and Lee, 1997 and references therein; Buryachenco, 2007). Unfortunately, because of the lack of full equivalence between compressible linear elasticity and incompressible linear viscous problems highlighted in Section 2, elastic Green function and their derivatives cannot be used readily for incompressible viscous materials. Nevertheless, the numerical evaluation of (8) is straightforward using simple quadrature techniques in standard mathematics applications such as MATHCAD or MATLAB. A MathCad worksheet to evaluate (8) for general anisotropic viscous materials is provided as online 9
ACCEPTED MANUSCRIPT Supplement Materials with this paper. The challenge remains in self-consistent simulations (below) where the integrations must be evaluated for many elements at every increment of deformation with a nontrivial amount of computation (8c, particularly) at each evaluation.
PT
3.2: Eshelby tensors and auxiliary quantities
RI
To present the solutions for the general viscous inclusion problems, two auxiliary tensors are defined: the second-order ij , related to the calculation of pressure (Lebensohn et al. 1998,
SC
their Eq.B14; Molinari et al., 1997, their Eq.31), and the fourth-order Tijkl (Molinari et al., 1987), sometimes referred to as Green’s interaction tensor (e.g., Lebensohn et al. 1998). They are respectively:
H x x 'n x ' dS x ' i
j
MA
NU
ij x H i , j x x ' dV x '
Tijkl x Gik , jl x x ' dV x '
G x x 'n x ' dS x ' ik ,l
j
(10a)
(10b)
TE
D
where the integrations are over the ellipsoid volume (first equal sign) or over the ellipsoid surface (the second equal sign) S x ' with n j being plane-normal unit vector component of the
AC CE P
area element dS on the ellipsoid surface. The second equal signs in (10a, b) are due to application of Gauss’s theorem (e.g., Korn and Korn, 1968, p.5.6.2). Two 4th order tensors, S and Π , referred to, respectively, as the symmetric and antisymmetric Eshelby tensors for viscous materials are defined as follows (Jiang 2013):
S J d : T : C and Π J a : T : C
(11)
d where C is the viscous stiffness (viscosity) tensor of the matrix, J ijkl is the same as in Eq.2f,
a J ijkl
1 ik jl jkil is the fourth-order anti-symmetric identity tensor as used in Jiang (2013, 2
2014). The operation “:” stands for double-index contraction between tensors. Eqs.11 are also equivalent to the following expressions used in Appendix B:
S Ts : C and Π Ta : C
10
ACCEPTED MANUSCRIPT s where Tijkl
1 1 a Tijkl T jikl Tijlk T jilk and Tijkl Tijkl T jikl Tijlk T jilk (Lebensohn et al., 4 4
PT
1998).
NU
SC
RI
An amazing result of Eshelby (1957) is that the mechanical fields inside the ellipsoid are uniform. This remarkable property, called the Eshelby uniformity (Zheng, et al., 2006), has played a key role in micromechanics. Note the Eshelby uniformity in ellipsoids (or ellipses in 2D) is a proven fact, not a conjecture, since Eshelby (1957). Eshelby (1961) did further conjecture ‘…Among closed surfaces, the ellipsoid alone has this convenient property…”. This conjecture of Eshelby’s has been proved for isotropic materials by Kang and Milton (2008) and Liu (2008). Eshelby’s conjecture remains open for anisotropic materials. As a consequence of the Eshelby uniformity, both ij and Tijkl are constant tensors for
MA
points within the ellipsoid (called interior points below). Accordingly, Sijkl and ijkl are also constant tensors for interior points. The mechanical fields in the matrix (exterior points) vary with space. To differentiate, the corresponding tensors for exterior points are denoted by ijE (x) ,
TE
D
E E E E , and ijkl (where superscript “E” stands for exterior field). Sijkl (x) , ijkl (x) , or simply ijE , Sijkl
AC CE P
3.3: Formal solutions for the interior and exterior fields of an incompressible linear viscous ellipsoid Finally, the mechanical fields inside the ellipsoid are related to the remote fields by the following set of equations: ε J d S1 : C1 : σ
(12a)
w Π : S1 : ε
(12b)
p Λ : C : S1 : ε
(12c)
1
where the mechanical field quantities inside the ellipsoid are represented by lowercase symbols ε , σ , w , and p , standing for, respectively, the strain rate tensor, the deviatoric stress tensor, the vorticity tensor, and the hydrostatic pressure; their uppercase counterparts Ε , Σ , W, and P shall stand for corresponding field quantities in the far-field matrix (Fig.2). The tilde over a quantity stands for the difference of the quantity between the field inside the ellipsoid and that in the far-field matrix, e.g., ε ε Ε . Eqs.12 collectively are called “solutions to an interior point Eshelby problem” or simply “the interior solutions”. 11
ACCEPTED MANUSCRIPT The exterior point fields are also expressed by a set of interaction equations as follows: (13a)
w(x) ΠE (x) : S1 : ε
(13b)
PT
ε(x) SE (x) : S1 : ε
(13c)
RI
p(x) Λ E (x) : C : S1 : ε
SC
This set is often referred to as “solutions to an exterior point Eshelby problem” or simply “the
NU
exterior solutions”.
MA
4. Viscous tensor identities
D
In this section we establish some previously unknown identities among viscous tensors, mostly for isotropic materials. These identities make the numerical evaluation of viscous inclusion problems more efficient and accurate.
TE
4.1: Interior field quantities
For the interior fields, ij and Tijkl , both constant tensors, can be expressed respectively
AC CE P
as (Lebensohn et al., 1998; Weinberger et al., 2005):
ij
a1a2a3 i z j 3 sin dd 2 0 0
(14a)
Tijkl
a1a2a3 Aik z j zl 3 sin dd 2 0 0
(14b)
where ai i 1, 2, 3 are the three semi-axes of the ellipsoid,
a1z1
2
a2 z2 a3 z3 , 2
cos sin and unit vector z is expressed in terms of spherical angles by z sin sin . These two cos integrals can be numerically evaluated efficiently using a product Gaussian quadrature (Jiang, 2013, 2014; Atkinson, 1982) or a Lebedev quadrature depending on the shape of the ellipsoid (see Qu et al., 2016).
12
2
ACCEPTED MANUSCRIPT From the definition of (7) it is easy to show that zi i 1 . Combining this with (14a), we a1a2a3 4
2
2
0 0
0 0
3 sin dd . Applying the identity
Michelitsch et al. (2003), we get the following identity:
kk 1
sin dd
4 due to a1a2a3
(15)
RI
(for general viscous materials)
3
PT
have kk
It is clear from (14a) that ij ji for anisotropic materials. For isotropic materials
SC
however, it can be shown (Appendix A) that i zi and the integrand in (14a) is odd with
ij 0 (if i j )
NU
respect to z if i j . This means: (for isotropic materials)
(16)
TE
D
MA
Since Bilby et al. (1975), it has been a common practice to use Eshelby tensor expressions for isotropic elastic materials to calculate the corresponding viscous Eshelby tensors by setting the Poisson’s ratio term ( in the following equations) to 0.5. In Appendix B, it is shown that the viscous and elastic symmetric Eshelby tensor components for isotropic materials are equal for off-diagonal components only. For diagonal components (i.e., components with indices like Siijj , including i j ), the following relations hold (see Appendix B for detail
ii
AC CE P
derivation):
1 el Sii11 Siiel22 Siiel33 0.5 3
el Siijj ii Siijj
0.5
(no sum)
(17a)
(no sum)
(17b)
where the superscript “el” stands for elastic. The antisymmetric Eshelby tensors for viscous and elastic isotropic materials are indeed equal: ( el ) ijkl ijkl
(18)
0.5
4.2: Exterior field quantities E For the exterior field solutions (13), the corresponding ijE , Tijkl , and Eshelby tensors E E , ijkl can be calculated numerically using integral expressions in (10) together with Green Sijkl
function expressions given in (8). The computation is however intensive even for isotropic materials (Eshelby 1959). Mura (1987) attempted to make the exterior solution for isotropic 13
ACCEPTED MANUSCRIPT
SC
RI
PT
elastic materials more explicit. He obtained rather long expressions (his Eqs.11.41 and 11.42) containing elliptic integrals (Mura 1987, p.85) and their derivatives that would still require considerable computational effort to evaluate. Thereafter, there have been two major advances for the exterior solution for isotropic elastic materials. First, Meng et al (2012) finally closed those unfinished derivatives in Mura’s expressions. They reduced the terms in Mura’s Eqs.11.41 and 11.42 that contain elliptic integrals and their derivatives to two simpler complete elliptic integrals of the first and second kind. This has greatly improved the efficiency and accuracy in numerical calculations of the exterior elastic field. The complete displacement field can be calculated using the far more explicit expression of Meng et al. (2012, Eq.21 there). Second, Ju and Sun (1999) developed a different and much more explicit approach to calculate S E ijkl [their
NU
G ijkl ] for isotropic elastic materials. Their approach allows the exterior strain and stress fields
(but not the complete displacement field) to be calculated even more efficiently. However, Ju and Sun’s (1999) approach cannot yield the complete displacement field because E ijkl is not
MA
calculated.
TE
D
The development below allows the complete viscous field quantities (deviatoric stress, strain rates, vorticity, and pressure) outside the inclusion to be computed efficiently and accurately. 4.2.1: Direct calculation of exterior field quantities
AC CE P
E To solve the exterior fields for viscous inclusions, we must obtain E ijkl , E ij , and Sijkl .
Let us investigate these quantities for incompressible isotropic viscous materials. First, because the Green functions for isotropic materials are available in analytic forms (9), it is straightforward, albeit a bit tedious, to get the following integral expressions: E Sijkl (x)
1 4
(19a)
xi x l jl x i x k ik x j x l il x j x k dV
(19b)
1 3 ij kl ik x j x l jk x i x l il x j x k jl x i x k 2 ij x k x l kl x i x j 15x i x j x k x l dV 3 2
r
E ijkl (x)
ijE (x)
3 8
1 4
1
r
1
r
3
where r x x ' , x i
3
jk
3x x i
j
ij dV
(19c)
xi x 'i , and the integration is over the volume of the ellipsoid x x'
x '12 x '2 2 x '32 2 2 1. a12 a2 a3
14
ACCEPTED MANUSCRIPT Second, using the two identities below:
xi x j x k 2 r
RI
xi x j x k xl 1 1 3 il x j x k jl x i x k kl x i x j 3 r 5 r x 'l
PT
x i x j 1 ij xi r3 3 r 3 x ' j r 2
SC
one can obtain the following relations with the Gauss theorem: xi x j 1 dV ij I 0 ij 3 r 3
(20a)
NU
xi x j x k x l 1 dV il jk jl ik kl ij ijkl 3 r 5
(20b)
MA
xi n dV x i x j x k nl , ij 2 j dS , and ijkl dS 3 r r r2
(20c)
TE
D
I0
AC CE P
Third, substituting (20) into (19) and then with simple variable substitutions to convert integrals over an ellipsoid surface to ones over a unit sphere, we get the following expressions:
1 4
E Sijkl (x)
E ijkl (x)
ijE (x)
1 ij kl 2 ik jl jk il il jk jl ik 3ijkl
1 jk il jlik ik jl il jk 8
(21a)
(21b)
1 ij 4
(21c)
where ij and ijkl are functions of x , which can be expressed explicitly as:
xi ai zi sin 2 2 2 0 x a z x a z x a z 2 2 2 3 3 3 1 11
2
ij (x) a1a2a3 0
2
ijkl (x) a1a2a3
aj
2
xi ai zi x j a j z j xk ak zk sin
d d
zl d d 2 2 2 2 a l x a z x a2 z2 x3 a3 z3 1 1 1 2
0 0
zj 3
5
15
(21d)
(21e)
ACCEPTED MANUSCRIPT and unit vector z is the same as in (14). With (21), solving the exterior fields eventually boils down to the evaluation of ij and
ijkl (21d and e), both of which can be performed efficiently and accurately with built-in
RI
PT
quadrature in Matlab, Maple, and Mathcad or with a combination of product Gaussian quadrature and Lebedev quadrature for extremely elongate or flattened ellipsoids (Qu, et al., 2016).
SC
4.2.2: Efficient calculation of exterior field quantities using Ju and Sun’s (1999) method The solution to exterior viscous inclusion problems in isotropic materials can be further
NU
E simplified. First, it is clear from Eq.19a that for incompressible isotropic viscous materials, Sijkl E E Sklij has the major symmetry, i.e., Sijkl . Whether this is also true for incompressible anisotropic
MA
materials remains open. Second, from (20a) it is clear that ij ji and kk 0 . Therefore ijE E 0 ) for isotropic materials. In summary, for must be symmetric ( ijE Eji ) and traceless ( kk
D
isotropic incompressible materials we have,:
(22)
TE
E E E Sijkl Sklij 0 , ij ji , kk 0 , ijE Eji , and kk
E From Eq.21a, it is easy to confirm that S kl 0 . This, together with the major symmetry
AC CE P
E of Sijkl , means that for incompressible materials, we have: E E E S kl Sij Sii 0 (sum applies to repeated Greek letter only)
(23)
Eqs.22 and 23 are now combined with Ju and Sun’s (1999) approach to greatly simplify the calculation of exterior field quantities. Ju and Sun (1999) develop an explicit approach to calculate the exterior-point symmetric Eshelby tensor for isotropic elastic materials [their G ijkl (x) ]. Comparing their Eq.4 (see also Li and Wang, 2008, p.105) with (19a) above, we get the following relation: E Sijkl (x) G ijkl (x)
0.5
ijE (x) kl
(24)
which can be regarded as a generalized relation of (B9) that we obtain for the interior field quantities in Appendix B.
16
ACCEPTED MANUSCRIPT E Applying the conditions of Eqs.23 and kk 0 to (24), we arrive at the following E E relations which allow all viscous tensors, ijE (x) , Sijkl (x) , and ijkl (x) , required for the solution
1 G ij (x ) , 3 0.5 E ij (x ) 1 G ij (x ) Gij ( x) , 0.5 3 0.5
ijE (x) kl ,
RI
(i j )
(25a)
except i j k l
i
j, i k , and j k
NU
E E E Siiii (x) Siijj (x) Siikk (x) ,
1 jk ilE jl ikE ik Ejl il Ejk 2
(25b) no sum
(25c) (25d)
MA
E ijkl (x)
(i j)
SC
E Sijkl (x) G ijkl (x)
PT
of exterior fields at a point, to be calculated simply from G ijkl (x) of Ju and Sun (1999):
TE
D
This greatly reduces the computation resources and improves the accuracy of the results. A MathCad worksheet to calculate exterior field quantities based on (25) is provided as online Supplement Materials with this paper.
AC CE P
5. Equations for rotational behaviors of ellipsoids in anisotropic viscous materials and principles of numerical implementation In this Section, we use the results of Section 4 to establish some new equations governing the motion of inclusions in anisotropic viscous materials and discuss numerical simulations of inclusion evolution. 5.1 Angular velocity of a viscous ellipsoid The rotation of an ellipsoidal inclusion in a viscous flow is determined by both the vorticity in the ellipsoid (12b), which determines the angular velocity of the material of the ellipsoid, and the migration of the principal axes through material due to ellipsoid’s shape change. The latter is described by a shear spin, denoted by w s , defined as (Goddard and Miller, 1967; Bilby and Kolbuszewski, 1977; Jiang 2012): ai2 a 2j 2 2 ij wijs ai a j w ij
ai a j
(no sum)
ai a j
17
(26)
ACCEPTED MANUSCRIPT The angular velocity tensor for an ellipsoid, written in ellipsoid’s coordinate system, ij , is then (e.g., Kocks, 1998; p. 395): ij wij wijs
PT
(27)
RI
Combining (12b) and (27), the angular velocity equation for the most general case of a deformable anisotropic inclusion in a dissimilar viscous anisotropic matrix is: (28)
SC
1 ij Wij ijmn Smnkl kl wijs
rigid in anisotropic
1 Wij ijmn Smnkl Ekl
MA
ij
NU
For a rigid inclusion, kl 0 (equivalent to kl kl ), (28) reduces to the following equation (29)
which is a generalization of Jeffery’s (1922) equation to an anisotropic viscous matrix. When the
D
1 matrix is isotropic, the term ijmn Smnkl Ekl becomes
TE
reduces to the well-known Jeffery’s result, ij
a I2 a J2 ij (Appendix C) and Eq.29 further a I2 a J2
rigid in isotropic
Wij
aI2 aJ2 ij (Jiang, 2012, Eq.15 aI2 aJ2
AC CE P
there). Here and below the notation adopted by Mura (1987) is used, whereby repeated lowercase indices are summed up from 1 to 3, while uppercase indices which take on the same numbers as the lowercase ones are not summed up. The behavior of an inviscid ellipsoid is another limiting case where the deviatoric stresses in the inclusion vanish (i.e., σ Σ ). Submitting this into (12a) and (12b), we get, after arranging: ε general J d S : Ε , 1
inviscid
ε general S1 J d : Ε 1
(30a)
inviscid
w general Π : J d S : Ε 1
(30b)
inviscid
ij
general inviscid
1 a2 a2 d d Wij ijkl I2 J2 J ijkl J S klmn Emn aI aJ
s
(30c)
5.2 Physical meaning of w , w , and justification for setting wijs wij when ai a j 18
ACCEPTED MANUSCRIPT
RI
PT
The vorticity in the inclusion w , which can be calculated from (12b), defines the rotation of the material element of the ellipsoid. w is always well defined whether the inclusion is a triaxial ellipsoid, spheroid, or sphere. Because the flow inside the ellipsoid is generally nonsteady, w can be decomposed to an internal vorticity and a spin in the sense of Astarita (1979) and Means et al (1980). The internal vorticity is the ‘objective’ (reference-frame invariant) part of vorticity that should be used to define the characteristics of flows (Jiang, 1999, 2010; Hobbs and Ord, 2015, p.77). Xiang and Jiang (2013) considered the decomposition of w in their investigation of the deformation of weak inclusions regarded as proxies for small-scale ductile s
NU
SC
shear zones. As pointed out earlier, the shear spin w is an entirely different quantity describing the phenomenon exhibited by a deforming ellipsoid whose principal axes swipe through the material as the ellipsoid changes shape (e.g., Goddard and Miller, 1967; Bilby and Kolbuszewski, s
MA
1977). It is called a ‘shear spin’ by Kocks (1998, p.385) in recognition that w arises from shearing between the instantaneous principal axes of the ellipsoid – wijs 0 if ij ( i j ) vanish. The shear spin must be considered for a deformable inclusion. However, the original definition of wijs is singular if the inclusion is a spheroid ( ai a j ). This occurs because a
TE
D
perfect circle does not have definable principal directions in the first place. Jiang (2007b) considers all possible singular cases that may be encountered in a numerical calculation of isotropic inclusion problems. Later, Jiang (2012) simplifies the singularity handling by setting wijs wij when ai a j which is the second equation in (26). Davies et al (2013) have
AC CE P
questioned the justification of this treatment. They instead suggested using a limit for wijs as ai a j based on the assumption that the orientation matrix Q is continuous at the ai a j
state. In the following, the second equation in (26) is further justified and it is shown that Q is intrinsically discontinuous at the ai a j state. The definition of wijs wij when ai a j first given in Jiang (2012) was based on the same line of reasoning as Jiang (2007b): Suppose an inclusion does go through a spheroid state at time t when two semiaxes are equal, say, a1 a2 (Fig.3). Call the partitioned flow field in the inclusion L= + w at time t (Fig.3). Immediately after t , at t t , the sectional shape of the inclusion must represent the infinitesimal strain ellipse in the inclusion. The two semiaxes ( a1 and a2 ) at t t are thus parallel to the 2 eigenvectors of ε , in the limit of t 0 . In an allisotropic system considered by Jiang (2007b), the eigenvectors of ε are the same as the eigenvectors of imposed far-field strain rate tensor. In general, the eigenvectors of ε do not coincide with that of the far-field strain rate tensor, but this does not affect the following argument. Immediately prior to the spheroid state, at t t , the inclusion’s sectional shape must be an ellipse which is reciprocal to the one at t t . This means that, at t t , the twosemiaxes are eigenvectors of ε (Fig.3). Since the eigenvalues of ε and ε are opposite in 19
ACCEPTED MANUSCRIPT sign but equal in magnitudes, it follows that the a1 - and a2 - axes at t t are parallel, respectively, to a2 - and a1 - at t t . In other words, the two semi-axes swap their identities but
PT
do not otherwise change orientations as the inclusion passes through the spheroid state (Fig.3). This behavior is fully captured by setting ij 0 which is equivalent to wijs wij . Numerically,
SC
RI
this treatment amounts to suspending the orientation update (achieved by ij 0 ), as the inclusion goes through a spheroid state, and updating the axial lengths only so that the inclusion’s new principal directions emerge. This fully justifies the second equation in (26) first introduced by Jiang (2012). The assumption that the orientation matrix Q is continuous at the ai a j stage is
NU
invalid. If is the angle between the major eigenvector of ε and one coordinate axis (Fig.3),
MA
cos sin sin cos then lim Qij t t = Qij t t and lim . As the left and t 0 x 0 sin cos cos sin right limits are not equal, Qij t is intrinsically discontinuous at the a1 a2 state.
D
One could make Q continuous if the ordering of a1 a2 is not imposed. In such a case,
TE
the two principal axes merely change relative lengths as the inclusion goes through the spheroid state with a1 a2 before t , a1 a2 at t , and a1 a2 after t . Even in this situation, we still have dQij
0 which is equivalent to setting ij 0 . A great deal of inconvenience arises if the a1 a2
AC CE P
dt
a1 a2 convention is removed. To use elliptic functions to calculate Eshelby tensors for
isotropic materials (Eshelby 1957; Mura 1987, p.78) for instance requires that the three semiaxes be ordered in terms of relative length. In a multi-inclusion system, imposing the a1 a2 a3 condition is the most straightforward way to keep track of the relative length of all inclusions. The intrinsic discontinuity of Qij if an inclusion does go through a spheroid state causes no trouble at all in numerical simulations if the complete definition (26) for shear spin is used. 5.3 Evolution of an inclusion through numerical calculations Given the formulation so far, the evolution of an ellipsoid inclusion in a general viscous material can be calculated numerically for any anisotropic viscous material in a given flow field. The procedure is as follows. The state of the ellipsoid at any time t is defined by its orientation matrix Q (Bunge, 1982; Jiang, 2007a, b) and the three semi-axes are expressed by a vector a1 a a2 . The evolution of the inclusion is found by solving the following two coupled a 3
differential equations (Jiang 2014): 20
ACCEPTED MANUSCRIPT dQ da ˆ , εa ΘQ dt dt
(31)
PT
where ˆij is a diagonal matrix made of the diagonal components of ij . Both ij and Θij are obtained from the above formulation. For a rigid ellipsoid, one only needs to solve the second
RI
equation.
NU
SC
Eqs.31 have been solved numerically by a fourth-order Runge-Kutta method (Jiang, 2007b, 2012), by using the Rodrigues rotation (e.g., Başar and Weichert, 2000, p.31) (Jiang 2013), and by a combination of the two (Jiang 2014). Clearly, it is the evaluation of associated Eshelby tensors in the process to obtain ˆij and Θij that takes the most computational time.
D
MA
Recently, Qu et al. (2016) have developed an optimal algorithm for calculating the Eshelby tensors of highly elongate and/or flat ellipsoids, combining the product Gaussian quadrature used in Jiang (2014) and the Lebedev quadrature.
TE
6. Isotropic systems
AC CE P
While for general anisotropic viscous materials, inclusion solutions must be expressed using integrals (12) and (13), it is shown in this section that, for all-isotropic (both inclusion and matrix are isotropic) systems, explicit analytic expressions can be derived directly from the general formulation. 6.1: Ellipsoid in 3D flows
From (12a and b), one can get the following explicit strain rate and vorticity partitioning equations (see Jiang (2013): ε J d S : C1 : Cinc J d : Ε
(32a)
w W Π : S1 : ε Ε
(32b)
1
inc
where C
stands for viscous stiffness of the inclusion. When the matrix and the inclusion are 1
inc
both isotropic, the C : C term in (32a) reduces to a scalar, r (not to be confused with r in Green functions), standing for the ratio of ellipsoid viscosity to the matrix viscosity and (32a) is simplified to (see Jiang (2013, Eq.11 there): 1
ε isotropic J d (r 1)S : Ε
(33)
21
ACCEPTED MANUSCRIPT All strain rate components in (33) can be solved explicitly (B7 and B8, Appendices B and C) to give:
1 2 r 1 Sijij
RI
, i j
(no sum)
(34b)
NU
where
ij
(34a)
SC
isotropic
ij isotropic
PT
1 1 (r 1) S2222 S2233 11 (r 1) S1122 S1133 22 1 22 isotropic (r 1) S2233 S2211 11 1 (r 1) S1111 S1133 22 33 11 22
11 isotropic
1 (r 1) S1111 S1133 1 (r 1) S2222 S2233 (r 1) 2 S1122 S1133 S2211 S2233 .
This
MA
reproduces Freeman’s (1987) equations with corrections by Jiang (2007b). All viscous Sijkl in (34) can be calculated using Eqs.17 and 18 from explicit expressions for elastic Eshelby tensors given by Eshelby (1957) and Mura (1987).
isotropic
aI2 aJ2 ij aI2 aJ2
(35)
TE
wij
D
In Appendix C, it is shown (Eq.C21) that:
wij
AC CE P
which combined with (34b) gives:
isotropic
ai2 a 2j 2 r 1 Sijij 2 ij (no sum) a a 2 1 2 r 1 S i j ijij
(36)
The shear spin is obtained by combing (26) with (34b) to get: s ij isotropic
w
ai2 a 2j ij , 2 2 a a 1 2 r 1 S j ijij i
(no sum)
(37)
Substituting (36) and (37) to (27) yields the following angular velocity equation for an isotropic ellipsoid embedded in an isotropic matrix: ai2 a 2j ai2 a 2j ij ij Wij 2 r 1 Sijij 2 2 2 2 a a a a 1 2 r 1 S isotropic j i j ijij i
(no sum)
(38)
Jeffery’s result is reproduced by letting r in (38). The angular velocity of an incompressible inviscid ellipsoid in an isotropic matrix is obtained by setting r 0 in (38) to get:
22
ACCEPTED MANUSCRIPT
ij
ai2 a 2j ai2 a 2j Wij 2Sijij 2 1 2Sijij ai2 a 2j ai a 2j ij
isotropic inviscid
(no sum)
(39)
PT
Mulchrone (2007) derives a 2D version of (37) for ellipses in isotropic plane flows.
ijd 2c ij , p 2m iI SiIjJ Jj
(40)
NU
1
SC
isotropic materials, the expressions are as follows:
RI
The deviatoric stresses inside the inclusion can be calculated from its strain rate tensor (obtained above) and the viscous stiffness. The pressure in the inclusion can be calculated using (12c). Because ij is diagonal (16) and Sijkl are zero except for terms with repeated indices for
MA
where the subscripts “c” and “m” stand for, respectively, the inclusion clast and the embedding matrix. As explained earlier, summation applies to repeated lowercase indices while uppercase indices take their lowercase values but are not summed.
TE
D
As already covered in Section 4.2, the mechanical fields outside the inclusion vary with space but can be easily obtained combining the formal solutions (13) with the relations (25) by which all relevant tensors can be calculated from the explicit approach of Ju and Sun (1999) efficiently and accurately.
AC CE P
6.2: Ellipses in 2D flows
An isotropic viscous ellipse in a 2D plane-straining flow can be regarded as an elliptic cylinder in 3D. Elastic Eshelby tensors for isotropic elliptic cylinders are already given in Mura (1987, p.80). With the Poisson’s ratio set to 0.5, the components are: el S1111
2R 1
R 1
2
el , S1122
1
R 1
2
el , S1133
el el el el S3311 0 , S3322 0 , S3333 0 , and S1212
where R
1 R R2 2R R2 el el el , S2211 , , S2233 , S 2222 2 2 R 1 R 1 R 1 R 1 R2 1
2 R 1
(41)
2
a ( a and b being the two semi-axes of the elliptic section with a b ). b
The corresponding non-vanishing viscous Eshelby tensor components can be calculated from (41) using (17) to get: S1111 S2222 S1122 S2211
R
R 1
2
23
(42a)
ACCEPTED MANUSCRIPT R2 1 2 R 1
11
(42b)
2
1 R , 22 R 1 R 1
(42c)
PT
S1212
RI
1 The non-vanishing components of Sijmn for elliptic cylinders under plane-straining flows 1 2D can be obtained using the identity Sijmn [where, as explained in more detail in Smnkl J ijkl
S
1 1212
R 1
1 1122
S
S
1 2211
R 1 4R
2
(43a)
2
2 R 2 1
(43b)
D
S
S
1 2222
MA
1 1111
1 1 ik jl jk il ij kl with i, j, k , l 1,2 ]. They are: 2 2
NU
2D deformations defined as: J ijkl
SC
2D Section 7, J ijkl is the fourth order deviatoric identity tensor for plane-straining incompressible
TE
With (42) and (43), we can express the interior strain rate and vorticity fields of the inclusion in terms of any given far-field plane-straining flow.
AC CE P
If the far-field plane-straining flow is a 2D general shear with vorticity parallel to x3 and the elliptic cylinder is parallel to the x3 axis as well, submitting (42) to (34) and (36) produces the following relations, equivalent to the results of Bilby and Kolbuszewski (1977):
2D 11
2D 22
R 1
2
R 2 2rR 1
11 ,
2D 12
R 1
2
rR 2 2 R r
2D 12
12 , w
1 r R 2 1 rR 2 2 R r
12
(44)
Submitting (42b) into (38) yields the angular velocity of the ellipse (Note for the elliptic cylinder 2D
2D
considered here the non-vanishing components are 12 21 ): 2D 12
isotropic
2 R 1 r R 1 2 R W12 12 2 R 1 rR 2 R r
(45)
which is equivalent to the result of Bilby and Kolbuszewski (1977, their Eq.15). The deviatoric stresses and pressure inside the inclusion can also be calculated similar to 3D cases (40). As ii and Sijkl are given in (42), the pressure in the ellipse, using (12c), is:
24
ACCEPTED MANUSCRIPT
p
2m11 1 r R 2 1
(46)
R 2 2rR 1
PT
which is equivalent to Schmid and Podladchikov (2003, their Eq.67).
NU
6.3: The work of Jaswon and Bhargava (1961)
SC
RI
The exterior fields can be obtained in a manner similar to the 3D case. Jin et al. (2011) give an explicit solution for elliptical cylinders as a special case of Ju and Sun (1999). Their solution can be combined with (25) to calculate the complete mechanical fields outside an inclusion. Another alternative to get the exterior fields is to apply the explicit solutions of Jaswon and Bhargava (1961) for plane-strain isotropic elastic inclusions to viscous conditions as described below.
MA
Analytic solutions for 2D elastic inclusion problems are found by Jaswon and Bhargava (1961). Although they use a complex variable formalism, most of their results are in explicit real-form expressions that can be readily extended to isotropic incompressible viscous materials.
D
The following expressions are from Eqs.19 and 26 of Jaswon and Bhargava (1961) recast in the current notation. They relate the constrained inclusion strain state ijc (Eshelby 1957) and
TE
the eigenstrains ij* :
1 R * * 11* 22 * 22 2 11 R 1 R 1
(47a)
c 22
1 R * * 11* 22 * 22 2 11 R 1 R 1
(47b)
12c
AC CE P
11c
R2 1
R 1
2
12*
(47c)
* Considering the incompressibility condition ( 11* 22 0 ) and using the relation ijc Sijkl kl* (c.f. 1 Eshelby, 1957), (47) immediately lead to all the expressions of Sijkl and Sijkl , given in (42) and
(43) respectively, which lead to the first two equations of (44) for the interior strain rates. Jaswon R 1 * and Bhargava (1961) also give an expression (their Eq.26) for inclusion rotation: w12c 12 R 1 2D R2 1 which readily leads to w12 2 12 and then the interior vorticity (the 3rd equation of 42) R 1 upon substitution of 12 12 E12 . This in turn leads to (45) upon consideration of the shear spin term. 25
ACCEPTED MANUSCRIPT Jaswon and Bhargava (1961) also discuss the limiting cases of a rigid inclusion and a void inclusion which lead to equations later rediscovered by Ghosh and Ramberg (1976) and Mulchrone (2007). The interior pressure can also be obtained from Eq.14a of Jaswon and Bhargava (1961)
PT
R 1
11* . This yields (46)
R 1 which follows from the 11 2
2R
incompressibility condition.
SC
1 1 upon substitution of 11* S1111 S1122 11
2m R 1
RI
which upon consideration of incompressibility becomes p
p ,
MA
NU
Solutions for the exterior fields are also included in Jaswon and Bhargava’s (1961). The expression for the exterior pressure field (their Eqs.13 and 21 combined, using the principle of superposition) is recast below: * 8m R sinh 2 sin 2 12* 1 11 2 R 1 cosh 2 cos 2 cosh 2 cos 2
D
where , are confocal elliptic coordinates which are related to the Cartesian coordinates by:
TE
x1 a 2 b2 cosh cos , and x2 a 2 b2 sinh sin . The eigenstrain terms are, for a 2
AC CE P
heterogeneous inclusion, replaced by an equivalent inclusion with
R 1 2 12 2
* 12
R 1
R 1 11 * 11
2R
and
to yield an explicit expression for the exterior pressure field:
p ,
8m R R 1 11 12 sinh 2 sin 2 1 R 1 cosh 2 cos2 2 R cosh 2 cos2 R 2 1
which can be further recast in the more familiar polar coordinates , , related to the Cartesian coordinates by x1 cos and x2 sin , as:
p , 8m
R R 1 11 12 B , 2 1 A , R 1 2R R 1
where
A ,
F 1 h cos 2 2F 26
(48)
ACCEPTED MANUSCRIPT F 1 h cos 2 2F
F 1 h 2h cos 2 , and h 2
a 2 b2 2
PT
B ,
B , takes positive sign when is in the first and 3 quadrants and negative sign when it is
RI
rd
NU
SC
in the second and fourth quadrants. 11 and 12 in (48) can also be expressed in terms of the far-field strain rates using (44). The expression of (48) is the most succinct real-form explicit expression available for the exterior pressure field.
AC CE P
TE
D
MA
Jaswon and Bhargava’s (1961) work has been overlooked by earth scientists who have, in their investigation of elliptic inclusions in 2D plane-straining flows of isotropic incompressible viscous materials, rediscovered equations available from Jaswon and Bhargava’s work. Ghosh and Ramberg (1976) obtained the angular velocity equation for a rigid ellipse in plane-straining flows of isotropic viscous materials, which is a limiting case explicitly discussed by Jaswon and Bhargava (1961, p. 677). Schmid and Podladchikov (2003) obtained close-form solutions (some expressed in complex variables) for the interior and exterior mechanical fields of a viscous ellipse. Their investigation of inclusions with a mantle layer of dissimilar viscosity is new contribution. Their other results overlap with Jaswon and Bhargava (1961). Mulchrone (2007) derives the angular velocity equation for an inviscid ellipse, which was also a limiting case explicitly considered by Jaswon and Bhargava (1961, p.677).
7. Two dimensional inclusion problems in anisotropic viscous materials 7.1 General statement
There is considerable literature on two-dimensional elastic inclusion problems since Jaswon and Bhargava (1961). Most of the work can be found in the books by Mura (1987) and Ting (1996). The work has been largely motivated by the attempt to avoid the ‘‘formidable” integrals (Yang and Chou, 1976) in the Green-function and point-force approach of Eshelby (Section 3). It has already been shown above that there is no longer any formidable integrals for isotropic viscous materials. The development in Section 4 together with the explicit approach of Ju and Sun (1999) and Jin et al (2011) has circumvented all formidable integrals. Only algebraic calculations and evaluation of two complete elliptic integrals of the first and second kind are necessary when one solves the interior and exterior fields of an ellipsoid of any shape. While it is straightforward to implement the general formulation in this paper numerically to solve any 3D inclusion problems in anisotropic materials, it is computationally intensive especially for the exterior field. Two-dimensional formulations are useful in at least three 27
ACCEPTED MANUSCRIPT
PT
respects. First, they lead to less intensive numerical calculations or even analytical solutions that allow inclusion behaviors to be investigated more thoroughly. Second, certain practical problems like deformation of thin elastic or viscous plates may be approximated by 2D plane-strain or plane-stress problems. Third, two-dimensional solutions, analytic or numerical, can serve as valuable benchmarks for more intensive and complex 3D inclusion systems.
AC CE P
TE
D
MA
NU
SC
RI
Theories for 2D anisotropic elasticity use a complex variable formalism pioneered by Lekhniskii (1950), Muskhelishvili (1953), and Green and Zerna (1968). Jaswon and Bhargava (1961) used the approach to study the 2D inclusion problem in isotropic elastic materials. Their work was extended to anisotropic elastic materials with cubic symmetry by Wallis (1964), with orthotropic symmetry by Bhargava and Radhakrishan (1964) and by Yang and Chou (1976) who gave real-variable expressions (see Mura 1987, p.142). Chen (1967) presented a solution procedure for monoclinic elastic materials based on the formulation of Green and Zerna (1960). Yang and Chou (1976), Hwa and Ting (1989), and Ting (1996) presented solution procedures for elastic inclusions in general anisotropic elastic materials under generalized plane deformation. As pointed out in the Introduction, it is not trivial to extend these works on 2D anisotropic elastic inclusions to incompressible viscous materials. The reason is that the “solutions” in these works are in fact menus, written in complex variable formalism, to be followed to get the solutions from the material’s anisotropic properties and boundary conditions, and it is often not trivial to execute the menu procedure. Take Chen’s (1967) “solution” for an elliptic inclusion in a monoclinic elastic material for example. According to the procedure, one must first solve 2 quartic equations (one for the matrix and one for the inclusion) of complex coefficients which are to be determined from the matrix and inclusion elastic compliances. Second, one solves two linear systems, each with four complex unknowns (one system is presented in Green and Zerna, 1968, p. 324 and the other in Eq.22-25 of Chen 1967). One must also consider possible degenerate cases if any of the quartic equations have repeated roots (Green and Zerna, 1968, p.205). The solutions in Hwa and Ting (1989) and Ting (1996, Chapter 10) require solving sextic equations of complex coefficients and manipulation of many matrices of complex components. Incompressibility may lead to new degenerate situations (Ting, 1996) that demand additional non-trivial treatment not available in the literature. In the following I use the same Green-function approach used in Section 3 for 3D inclusions to develop a 2D inclusion formulation for elliptic inclusions in anisotropic viscous materials. The formulation leads to integral expressions that can be evaluated efficiently in standard mathematic applications. 7.2 Generalized plane flows in anisotropic viscous materials Before proceeding, it is important to point out that a plane-straining flow in which the velocity field is confined to, say the x1 x2 plane [with velocity components being u1 u1 x1, x2 , u2 u2 x1 , x2 , and u3 0 ] is not possible in a general sense in anisotropic materials. This is
28
ACCEPTED MANUSCRIPT because the in-plane motion and anti-plane motion are not always separable in an anisotropic material (Yang and Chou, 1976; Ting, 1996). A flow field in which the anti-plane velocity component u3 does not vanish but nevertheless all velocity components depend on x1 and x2
PT
only [i.e., u1 u1 x1, x2 , u2 u2 x1 , x2 , and u3 u3 x1, x2 ] is called a generalized plane flow. In a generalized plane flow, the shear strain rates 13 and 23 do not vanish although 33 0 .
RI
Generalized plane flows are really special types of 3D flows.
SC
A plane-straining flow can exist in an anisotropic material only if the material has a symmetry plane parallel to the x1 x2 plane. Thus, triclinic materials can never have plane-straining
MA
NU
flows in them. If materials with monoclinic or higher symmetries are oriented such that no symmetry plane is parallel to the symmetry plane of the imposed velocity field, plane-straining flows will not occur in these materials. By inference, the occurrence of monoclinic flows (e.g., Jiang and Williams, 1998; Lin et al., 1998; Passchier, 1998) also requires that the flowing rock material have one symmetry plane parallel to the symmetry plane of the imposed velocity field.
D
The formulation below is valid for plane-straining flows only. Generalized plane flows can be treated using the formulation of Section 3.
TE
7.3 Formulation
The two Green functions, Gij2 D and H i2D respectively, (the superscript “2D” stands for
2D ij
G
AC CE P
the two-dimension case) are the inverse Fourier transforms of (6) in 2D space, respectively:
1
Aij ( z) 1 2 exp iK x dK and H i2 D 2 2 4 K 4
i i exp iK x dK K
Their derivatives, Gik2 D, jl and H i2, Dj , which are used in (10), can be obtained by simple differentiation and are expressed as: Gik2 D, jl
H
2D i, j
1
zzA
4 2
j l
1 4 2
ik
exp iK x dK
(49a)
z i
j
exp iK x dK
(49b)
Submitting (49) into (10), noting that the integrations are over the ellipse area for the 2D case, 2D and carrying out the explicit integration, we obtain the following expressions for ij2 D and Tijkl : 2D Tijkl x
ab
0
Aik z j zl B x, z d
(50a)
29
ACCEPTED MANUSCRIPT ab
0
i z j B x, z d
(50b)
i xz cos 1 1 where z , B r , z 2 2 sin 2 xz
stands for the real part of the complex function in the parentheses. Weinberger et al. (2005)
RI
with 2 a 2 cos2 b2 sin2 , and
PT
ij2 D x
SC
present a detailed procedure for getting (50a) from (49a). Using a similar procedure one obtains (50b) from (49b). For points inside the ellipse (the interior inclusion problem), 2 x z 0 and 1
2
NU
B x, z
2
2D , rendering constant ij2 D and Tijkl as expected. For exterior points (the exterior
TE
D
MA
inclusion problem), the two tensors vary in space. From (50) the corresponding Eshelby tensors for the interior and exterior points can be readily obtained which together with (12) and (13) solve the inclusion problem. The numerical evaluation of (50) can be done using built-in quadrature in any standard mathematics applications such as MAPLE, MATHCAD, or MATLAB, and the computation is efficient as only line-integrals are involved. 7.4 Application to materials with a planar anisotropy
AC CE P
Eqs.50 apply to any monoclinic material with symmetry plane parallel to the x1 x2 plane and subjected to plane-straining flows in the x1 x2 plane. As an illustration of using this 2D formulation, let us consider a viscous material of planar anisotropy specifically. Well-layered or foliated rocks (e.g., Johnson and Fletcher, 1994, p. 419) and minerals with a dominant slip system such as mica (Chen et al., 2014) have been treated as an idealized material with two distinct viscosities: a normal viscosity n that measures the resistance of the material to pure shearing along and perpendicular to the anisotropy plane and a shear viscosity s measuring the resistance to shearing parallel to the anisotropy. Fletcher (2004, 2009) obtains real-form analytic solutions for an isotropic viscous elliptic inclusion enclosed in such a planar anisotropic matrix undergoing plane-straining flows. Numerical calculations are compared with his analytic results. It should be pointed out that heterogeneous strain near the inclusion will lead to distortion of the planar anisotropy around the inclusion (Dabrowski and Schmid, 2011; Griera, et al., 2011, 2013). Such heterogeneity is not captured by the solution here or that of Fletcher (2004, 2009). As discussed above, in order for plane-straining flow to occur, the anisotropy plane (layering) must be parallel to the x3 axis. When the anisotropy plane is perpendicular to the x2 axis, the viscous stiffness tensor components can be written in the following block-matrix form:
30
ACCEPTED MANUSCRIPT 1 2D J 12 kl m 2D J 22 kl
(51)
PT
cijkl
2D J 11kl 2n 1 J 2 D 12 kl m
where the lowercase cijkl is used for 2D viscous stiffness tensor with all indices varying from 1 to
1 1 ik jl jk il ij kl ( i, j, k , l 1,2 ), m n is a measure of the strength of s 2 2
SC
2D J ijkl
RI
2D 2, J ijkl is the fourth order deviatoric identity tensor for plane-straining deformations defined as:
NU
anisotropy (Treagus, 2002; Fletcher 2004, 2009; Chen et al., 2014) which also corresponds to a similar concept in transverse isotropic elastic materials (Brace 1965). In (51) each component is itself a 2 2 matrix.
MA
In the ellipse’s coordinate system at angle relative to the anisotropy plane (Fig.4), the stiffness tensor components can be obtained by tensor transformation, c 'ijkl QiaQ jbQkcQld cabcd ,
TE
D
cos sin where Q . The non-vanishing components of c 'ijkl can be calculated and they sin cos are:
AC CE P
1 c '1111 = c '2222 = c '1122 =n cos2 2 sin 2 2 , m 1 c '1112 = c '1222 n 1 sin 2 cos 2 , m
(52)
1 c '1212 n sin 2 2 cos2 2 m 2D The corresponding compliances can also be obtained using the identity c 'ijmn s 'mnkl J ijkl . The
non-vanishing components are: s '1111 = s '2222 = s '1122 =
s '1112 = s '1222
s '1212
m 1 m 1 cos 4 , 8n
m 1 sin 4 ,
(53)
8n
m 1 m 1 cos 4 8n
31
ACCEPTED MANUSCRIPT which are equivalent to Fletcher (2009, Eq.5 there).
sin 2 2 m cos2 2 , A12 A21 n 1 m sin 4 2 , m 2m
n
m sin 2 2 cos2 2 m
RI
A22
n
(54)
SC
A11
PT
From (52), the corresponding Christoffel stiffness tensor components defined in (5a) can be calculated:
which in turn leads to the following expressions for Aik and i defined in (7) through matrix
n
ij
zi z j
m 1 cos m 1 cos(4 3 ) 2 m 1 sin m 1 sin(4 3 )
D
ζ
m
MA
Aij
NU
inversion:
TE
where m sin 2 2 cos2 2
1
(55a)
(55b)
(55c)
AC CE P
Substituting (55) into (50), the mechanical fields inside and outside the ellipse can be calculated. Below we consider the case in which the remote flow field is a plane-straining general shear with the shear plane parallel to the anisotropy and compare the numerical calculations based on the above with Fletcher’s (2009) analytic solutions.
7.5 Comparison with Fletcher (2004, 2009) Let the pure shearing strain rate and simple shear rate parallel to the anisotropy plane (the x1 x2 coordinate system in Fig.4) be respectively and . In the inclusion coordinate system (the x '1 x '2 coordinate system in Fig.4) at angle relative to anisotropy plane, the remote flow field strain rate and vorticity components are, expressed as functions of , , and :
11 , , 22 , , cos 2 sin 2 2
32
(56a)
ACCEPTED MANUSCRIPT
12 , , 21 , , sin 2 cos 2 2
(56c)
2
PT
W12 W21
(56b)
SC
m R 2 1
r m R 2 1 2 R
E12 , , , m, r
MA
12 , , , m, r
2 mR E11 , , , m, r R 2r mR 1 2
NU
11 , , , m, r
RI
The analytic expressions obtained by Fletcher (2009, his Eqs.10, 14, and 15) for the mechanical fields inside an isotropic elliptic inclusion can be recast into the following forms:
mR
11
2n R 2 1
11 , , r11 , , R 2r mR 1 m
TE
p
n R 2 1
D
m R 2 1 R2 1 w12 , , , m, r 2 12 E 12 , , , m, r R 1 r m R 2 1 2 R
2
(57a)
(57b)
(57c)
(57d)
as:
AC CE P
where r is the ratio of the inclusion viscosity to the matrix n and E ij , , , m, r are defined
E ij , , , m, r rEij , , Eij , , m One can easily confirm that (57) reduce to their isotropic versions in (44) and (46) when m 1 . Fig.5 presents the results for an inclusion ( R 5 , r 10 , and m 25 ) in plane-strain general shear ( 0.5 and 1 ) obtained from numerical calculations using the 2D formulation of this paper implemented in MATHCAD (solid lines) and those calculated directly using (57) (dashed lines). The results agree perfectly. The external mechanical fields can be calculated readily using the formulation here. Fig.6 shows the pressure fields normalized against the matrix deviatoric stress invariant around a strong inclusion ( R 3 , r 100 ). The inclusion long axis is horizontal and the far-field flow is pure shearing with horizontal stretching and vertical shortening. Fig.6a is for an isotropic matrix ( m 1 ). Figs.6b-e are when the matrix is anisotropic with m 10 . The anisotropy plane is horizontal in Fig.6b and increases by 30° increments clockwise until 90° in Fig.6e. It is seen clearly that the pressure perturbations are significantly subdued in an anisotropic matrix than in 33
ACCEPTED MANUSCRIPT an isotropic one (Fig.6a). This has implications for the study of tectonic pressure fluctuations during metamorphism.
PT
7.6: Comparison with 3D inclusion behaviors
NU
SC
RI
Although numerical calculations of 2D inclusion problems are far less intensive than 3D inclusion problems, most natural deformations involving heterogeneous inclusions cannot be treated as two-dimensional problems (see Li and Jiang, 2011 for more discussion). A question that naturally arises is: how well can 2D behaviors be applied to the understanding of 3D inclusion behaviors? For example, from Fletcher’s (2009) analytic solutions, one concludes that the rotational behavior of a rigid ellipse is unaffected by the planar anisotropy in the matrix. This is clear by considering the limit of (57c) as r . It reduces to the Ghosh and Ramberg (1976) relation for rigid inclusion in isotropic matrix. Is it also true for 3D ellipsoids? Is the rotation of a rigid ellipsoid unaffected by the presence of a planar anisotropy in the matrix?
MA
To answer this question, one can examine the general equation in (29) for a rigid ellipsoid 1 in anisotropic viscous matrix. As the term ijmn Smnkl in (29) must certainly depend on the
TE
D
anisotropy of the matrix, it cannot be true in general for rigid inclusion behaviors to be independent of matrix anisotropy. Is it true for the special case of planar anisotropy under the plane-straining flows? We can use the general formulation in Section 3 to answer this question.
Cijkl
AC CE P
In 3D, the viscous stiffnesses of a planar anisotropic material with the anisotropy plane normal to the x2 axis can be expressed in a block matrix form similar to (51):
d J 11kl 1 2n J 12d kl m J 13d kl
1 d J 12 kl m J 22d kl
1 d J 23kl m
J 13d kl 1 d J 23kl m J 33d kl
(58)
d where J ijkl are the 3D 4th-order deviatoric identity tensor as defined in (2).
Let us consider the rotational behavior of a rigid ellipsoid ( a1 : a2 : a3 ) with its a3 semi-axis parallel to the x3 axis embedded in a viscous matrix whose rheology is defined by (58). Suppose the far field flow is a plane-straining general shear like the 2D case considered above. Fig.7 presents the angular velocities (normalized against simple shearing) for spheroids ( a1 a2 ) in a plane-straining flow as functions of matrix anisotropy measured by m . The distinct
a3 axes of the spheroids are parallel to the vorticity of the flow (the x3 axis). The plane-straining flow in the x1 x2 plane is like the 2D case above with 0.5 and 1 along the anisotropy. A 34
ACCEPTED MANUSCRIPT circular cylinder ( 1: 1: ) rotates at a constant rate of half the shearing, independent of m , as (57c) predicts. The angular velocities depend on the relative lengths of a3 and m . A shorter a3
PT
and higher m imply a faster rotation rate. The angular velocity curves vary between the curve for a cylinder ( 1: 1: ) and the one for a disk with diminishing a3 ( 1: 1: 0 ). The same conclusions
RI
can be drawn for non-spherical objects. Fig.8a presents the angular velocities of a rigid triaxial ellipsoid ( 6 : 1: 5 , with the a3 axis parallel to the x3 axis) as a function of the orientation of the ellipsoid long axis with respect to the anisotropy plane ( x1 axis) for different matrix m . A
SC
greater m corresponds to a higher angular velocity at any given orientation, except near the stable orientation of the object. Fig.8b plots the angular velocity variations of the ellipsoid as its a3 axis varies from infinity to 0.01 for a fixed matrix m =10. A relatively shorter a3 generally
NU
corresponds to a higher angular velocity. Fig.9 shows the angular velocities of elliptic cylinder and ellipsoids, all with sectional aspect ratio of 2:1 on the x1 x2 plane but different a3 lengths, in a
MA
matrix with m =10 undergoing simple shearing. The angular velocity increases as a3 varies from infinity to nearly 0.
AC CE P
TE
D
The results above are based on a specific material with planar anisotropy and with the anisotropic plane parallel to the shearing. The behaviors of viscous inclusions in general anisotropic materials are terra incognita. There is indeed a great deal to be learnt using the formulation of this paper. The results from the simple example examined here already suggest that anisotropy affects the behavior of rigid and deformable inclusions significantly. In a 3D case, the ellipsoid can have any initial orientation with principal axes oblique to the far-field flow and the material anisotropy.
8. Applying the inclusion theory to lithosphere deformation So far, we have considered a single elliptic or ellipsoidal inclusion in an infinite homogeneous matrix. To apply the inclusion theory to the deformation of rocks in Earth’s lithosphere requires some assumptions. First, one must address the issue of non-linear rheology of the lithosphere as the inclusion theory is based on linear (viscous or elastic) constitutive behaviors. Second, rheologically distinct elements are irregularly shaped. How a theory based on ellipsoidal and elliptic inclusions may be applied to them also requires justification. Third, rock masses are heterogeneous poly-element aggregates. How and in what sense can they be treated as a “uniform matrix” needs to be clarified. Finally, assumptions like “an infinite matrix” and perfect bonding at inclusion interfaces used in the inclusion theory are not always justified in natural cases. We discuss these issues in more detail below. 8.1: Non-linear rheology and partitioning equations
35
ACCEPTED MANUSCRIPT
(59)
NU
ε(σ) M(sec) (σ) : σ , or σ(ε) C(sec) (ε) : ε
SC
RI
PT
There are no exact solutions like (12) and (13) for nonlinear materials such as the Earth’s ductile lithosphere which is often regarded as power law viscous (Kohlstedt et al., 1995; Tullis, 2002). However, it has been demonstrated that (13) can be used as approximate solutions for nonlinear materials using the linearization approach (Molinari et al., 1987; Lebensohn and Tomé, 1993; Ponte Castañeda 1996; Masson et al., 2000; Molinari 2004). The linearization approach consists in using a linear-form constitutive relation to approximate the actual nonlinear constitutive behavior in the vicinity of a given stress or strain (rate) state. For a power-law viscous material, the nonlinear constitutive equation is often written in a pseudo-linear form as (e.g., Hutchinson, 1976, p.104, his Eq.2.5):
where M(sec) (σ) and C(sec) (ε) are referred to, respectively, as the secant viscous compliance tensor and the secant viscous stiffness tensor (e.g., Lebensohn and Tomé, 1993). C(sec) is
D
MA
equivalent to the effective viscosity (e.g., Ranalli, 1995, p.78). M(sec) and C(sec) are no longer material constants but depend on the current state of deviatoric stress or strain rate. The linearization approach seeks the following linear-form relations: ε M : σ ε0 , or σ C : ε σ0
TE
(60)
to approximate (59). C and M in (60) are respectively the linearized stiffness and compliance
AC CE P
tensors, ε 0 and σ 0 are associated pre-strain rate and pre-stress terms. The two equations in (60) 1
are equivalent if C M , σ0 C : ε0 , and ε0 M : σ0 . A variety of linearization schemes have been proposed, differing in how M and ε 0 (or
C and σ 0 ) are constructed. For example, the tangent linearization uses the first-order Taylor expansion of (59) below:
ε(σ ) ε(σ 0 )
ε ε ε : σ σ0 : σ ε( σ 0 ) : σ0 σ σσ0 σ σσ0 σ σσ0
(61a)
σ( ε ) σ( ε 0 )
σ σ σ : ε ε0 : ε σ( ε 0 ) : ε0 ε ε ε0 ε ε ε0 ε ε ε0
(61b)
where M(tan)
ε σ and C(tan) are the tangent compliance and tangent stiffness σ σσ0 ε ε ε0
tensors respectively. For a power-law viscous material, (61) can be expressed as:
36
ACCEPTED MANUSCRIPT 1 ε(σ) nM(sec) (σ0 ) : σ ε0 , σ(ε) C(sec) (ε 0 ) : ε σ0 n
(62)
SC
RI
PT
1 where n is the power-law stress exponent, and ε0 1 n ε(σ0 ) , σ0 1 σ(ε 0 ) are the n back-extrapolated strain rate and stress terms (Hutchinson 1976; Molinari et al, 1987; Lebensohn and Tomé, 1993). If one uses a linearized stiffness or compliance, such as M(tan) , (12) and (13) are approximate solutions for the nonlinear viscous material. For simple presentation, (12a) is commonly recast in the following form (e.g., Jiang, 2014) with linearized stiffnesses and compliances used:
NU
ε H : σ , σ H : ε
(63)
where H S1 J s : M and H C : S1 J s .
MA
1
The following partitioning equations are obtainable from (63) (c.f. Jiang, 2014, Eq.10
D
there):
ε( k ) A( k ) : E ; σ( k ) B( k ) : Σ
TE
(64)
where A( k ) J d S( k ) : nmC(sec) 1 : C( k )(sec) J d : J d nm 1 S( k ) is the strain-rate 1
AC CE P
partitioning tensor and B( k ) nmC( k )(sec) : A( k ) : C(sec) 1 is the stress partitioning tensor, and nm is the power-law stress exponent of the matrix. The super- or sub- script “k” stands for the k-th element. Alternatively, the partitioning equations are expressed in terms of linearized stiffnesses or compliances (Jiang 2014): ε( k ) A
where A
(k )
(k )
(k )
H
(k )
( k ) 1
H( k ) C
partitioning tensors, A
and β
σ( k ) B
:Eα;
M
(k )
( k ) 1
(k )
(k )
: Σβ
: H( k ) C and α
(k )
(65)
( k ) 1
H( k ) C
being 4th-order and α second order; B
(k )
: Σ0( k ) σ0( k ) are strain rate
H
(k )
: E0( k ) ε0( k ) are stress partitioning tensors, B
37
M
(k )
: H
( k ) 1
(k )
M
being 4th-order and
ACCEPTED MANUSCRIPT β
(k )
second order. Note that the secant-based partitioning tensors ( A( k ) , B( k ) ) are distinct from (k )
(k )
,α ,B
(k )
,β
(k )
).
PT
the linearization-based partitioning tensors ( A
The partitioning equations relate the local deformation in a rheologically distinct element
RI
(RDE, regarded as an ellipsoidal inclusion) to the bulk deformation of the matrix.
SC
8.2: Irregular shape of rheological heterogeneities
AC CE P
TE
D
MA
NU
The simple reason for the common practice of regarding heterogeneous elements as Eshelby ellipsoids in micromechanics (e.g., Mura, 1987; Nemat-Nasser and Hori, 1999; Li and Wang, 2008) is because exact solutions for inclusion problems are limited to ellipsoids and cuboids in an infinite matrix and the Eshelby uniformity is limited only to ellipsoids (Eshelby, 1961; Mura, 1987). Where an inclusion is non-ellipsoidal, the Eshelby solutions (12) have been taken as an approximation for the average fields in the inclusion. It has been hard to justify this practice quantitatively (see Mancktelow, 2011; Jiang 2013 and references therein for some discussions). Budiansky and Mangasarian (1960) appear to be the first to treat individual grains in a polycrystalline material as ellipsoidal inhomogeneities embedded in the polycrystalline medium which is taken to be mechanically homogeneous (Qu and Cherkaoui, 2006). Their treatment has been followed in subsequent wide applications of the Eshelby formalism in micromechanics to obtain the macroscopic effective mechanical properties of a material from the properties of its constituent elements (Mura, 1987, p.421-433; Nemat-Nasser and Hori, 1999, Part I; Molinari, 2002). The practical success of this approach gives confidence on its use. Dislocations, stacking faults, cracks, weakened zones, and other discontinuities have also been treated as Eshelby “ellipsoids” (e.g., Mura, 1987, p.15-20, 240-379; Rudnicki, 1977; Healy et al., 2006; Exner and Dabrowski, 2010). In geology literature, rigid or deformable clasts have commonly been treated as ellipsoidal objects (e.g., Ramsay, 1967, p.209-221; Dunnet, 1969; Gay, 1968; Ghosh and Ramberg, 1976; Ježek et al., 2002). Many physical experiments (e.g., Ferguson, 1979; Arbaret, et al., 2001; Ghosh and Ramberg, 1976) and some theoretical considerations (Bretherton, 1962; Wallis, 1977) find that the rotational behavior of rigid non-ellipsoidal objects is very close to that of their best-fit ellipsoids. The assumption allows modern micromechanics to be applied to the deformation of heterogeneous materials like Earth’s lithosphere. It is also worth pointing out that the shape of an ellipsoid can be varied to fit a wide range of rheologicallydistinct elements. By changing the relative length of its three semi-axes, an ellipsoid varies from rod-like ( a1 , a2 a3 0 ), to spherical ( a1 a2 a3 ), to plane-like ( a1 a2 a3 0 ), and many other shapes. In a non-ellipsoid inclusion, the Eshelby tensor is no longer a constant but an Eshelby Tensor Field (ETF). Zou et al (2010) recently obtain explicit expressions of the ETF and its average for a wide variety of non-elliptical inclusions including arbitrary polygonal ones and 38
ACCEPTED MANUSCRIPT
8.3: Heterogeneous matrix and homogenization
SC
RI
PT
those characterized by the finite Laurent series. They conclude that the ETF inside a nonelliptical inclusion can be quite non-uniform. However, using the elliptical approximation as the average of the ETF is valid for convex non-elliptical inclusions. The approximation becomes less acceptable for non-convex non-elliptical inclusions. Where the shape effect is important, one must take the numerical challenge to address the so-called Eshelby’s problem of non-ellipsoidal inclusions for which there has been progress in isotropic elastic materials (e.g., Zou et al., 2010). In structural geology and tectonics studies, mechanically more competent elements tend to exhibit convex bodies. It is therefore more acceptable to regard Eshelby solutions as the average fields in these elements.
AC CE P
TE
D
MA
NU
In micromechanics, the heterogeneous rock mass composed of rheologically distinctive elements (RDEs hereafter) is replaced by a hypothetical homogeneous equivalent (or effective) medium (HEM). The mechanical interaction between a specific RDE (regarded as an Eshelby ellipsoid) and the surrounding matrix is then replaced by a heterogeneous ellipsoid with a rheologically uniform HEM. The rheological properties of the HEM at a point are taken to be the overall effective properties of all RDEs contained in a representative volume element (RVE) centered at that point. The procedure of getting HEM properties from the properties of constituent RDEs is called homogenization. Most homogenization methods are described in the books by Mura (1987), Nemat-Nasser and Hori (1999), and reviews by Zaoui (2002), Molinari (2002), and Mercier and Molinari (2004). Jiang (2014) summarizes the self-consistent homogenization method for poly-element continua based on Lebensohn et al (2003, 2004). A more complete summary of self-consistent homogenization is given below briefly. When the deforming rock mass is made of a more or less homogeneous matrix containing dilute heterogeneous inclusions, the contribution of the inclusions to the bulk rheology is negligible because the interactions among inclusions are insignificant. The HEM property is essentially the matrix property and the Eshelby solutions (63) and (64) can be used directly to find the mechanical fields inside individual inclusions (using the interior solutions) or to get the mechanical fields surrounding an inclusion (using the exterior solutions). As the inclusion volume fraction increases, interactions among inclusions become significant (e.g., Ildefonse et al., 1992, Mandal et al., 2003). A self-consistent homogenization is necessary to obtain the effective properties of the HEM, which is conducted by considering a representative volume element (RVE) of the material sufficiently large so that it statistically represents the material. It is also tacitly assumed that the constituent RDEs in the RVE are uniformly mixed and have a fractal-like structure in shape and orientations regardless of sizes. A salient reference medium is defined. Where the rock mass has a connected matrix containing heterogeneous inclusions, the Mori-Tanaka (1973) approach can be adopted. For a heterogeneous poly-element aggregate without a connected matrix, the reference medium is the HEM. In either case, each element is regarded as an ellipsoidal inclusion surrounded by the 39
ACCEPTED MANUSCRIPT reference medium and subjected to remote loading ( Em for strain rate and Σm for stress) which are, respectively, the average strain rate and associated deviatoric stress tensors in the reference material. The self-consistent homogenization works as follows. Consider a representative volume element (RVE) made of N different constituent elements and label them by k (k=1, 2, … N). The
where H
: σ( k ) Σm
(66)
RI
(k )
(k )
SC
ε( k ) Em H
PT
average strain rate ε ( k ) and stress σ ( k ) in a given element is given by (63):
is evaluated at Em .
NU
The strain rate partitioning equations follow from (64) and (65) as:
ε( k ) A
(k )
: Em α
(67a)
MA
ε( k ) A( k ) : Em (k )
(67b) means volume average over
D
Imposing the consistency condition ε ( k ) E on (65), where
1
:E
AC CE P
Em A( k )
TE
the RVE, gives:
Em A
(k )
1
:E A
(k )
1
: α
(68a)
(k )
(68b)
Combining the nonlinear response of each element (59), σ( k ) C( k )(sec) : ε( k ) , with (67a) and (68a) yields σ( k ) C( k )(sec) : A( k ) : A( k )
1
: E . Imposing the stress consistency condition
σ( k ) Σ on this equation gives the following expression for the homogenized secant viscosity: (sec)
C
C( k )(sec) : A( k ) : A( k )
1
(69) (k )
Combining the linearized response of each element (48), σ( k ) C
: ε( k ) σ0( k ) , with
(67b) and (68b) and then imposing the stress consistency condition σ( k ) Σ gives the following expressions for the linear-form of homogenized stiffness and back-extrapolated term: C C
(k )
Σ0 C
:A
(k )
(k )
:α
(k )
: A σ0
(k )
(k )
1
C
(k )
:A
(k )
: A
40
(k )
(70)
1
: α
(k )
ACCEPTED MANUSCRIPT A similar procedure can be followed with the stress partitioning equations from (64) and (65) to get the following homogenized compliances and the back-extrapolated term: M
M M
(k )
:B :β
(k )
(k )
: B( k ) : B( k )
: B
(k )
1
(71a)
1
ε 0( k ) M
(k )
:B
(k )
: B
(k )
1
: β
(71b)
(k )
SC
E0 M
(k )
( k )(sec)
PT
(sec)
RI
M
NU
Eqs.71b are the homogenization equations used in the self-consistent MOPLA of Jiang (2014). An initial set of guesses for the macroscale stiffness and compliance quantities are made to calculate them iteratively (Jiang 2014).
relations [like M nM
(sec)
MA
For a poly-element material, its secant properties, C
(sec)
and M
(sec)
, do not have simple
] with the linearized properties, C , M , Σ 0 , and
E0 that one may
AC CE P
TE
D
expect of mono-phase polycrystals (Hutchinson, 1976; Lebensohn and Tomé 1993). The rheological response of a general deforming heterogeneous material varies with time because fabric develops in the material as strain increases. Such a material is not expected to have a constant stress exponent. The common use of a steady-state, always-isotropic power-law rheology to represent Earth’s lithosphere in geodynamic models is likely unrealistic. Cook et al. (2015) arrived at a similar conclusion . 8.4: Inclusions in a finite space
To apply the inclusion theory to lithospheric deformation, the conditions of scale separation (Zaoui 2002; Jiang 2014) must be satisfied. To state simply, this requires that the size of an average RDE be small compared to the size of the RVE which itself must be small compared to the scale of the deformation. The scale separation conditions justify the use of the inclusion theory developed for an infinite body. Where the size of the inclusion phase becomes appreciable compared to its distance from the deformation boundary, boundary effects have to be considered. Both translational and angular velocities of the inclusion are affected (e.g., Goldsmith and Mason, 1967). This constitutes inclusion problems in a finite space that has been considered by Kinoshita and Mura (1984). Existing work is still limited to the spherical shape. Li et al. (2007) obtained closed-form solutions for the stress fields of a spherical inclusion concentrically placed in a finite spherical space. Zou et al. (2012) proposed a superposition framework to study inclusions in a finite space. Further progress to these finite space problems would enable better estimation of the mechanical properties for actual materials.
41
ACCEPTED MANUSCRIPT On the other hand, analogue experiments and numerical modeling have improved understanding of the rotational behaviors and surrounding flow patterns of rigid inclusions in confined simple shear (e.g., Marques et al., 2005a, b).
PT
8.5: Interface properties
AC CE P
TE
D
MA
NU
SC
RI
It is well known that the mechanical interactions between the inclusion and the embedding matrix and therefore the overall effective properties of heterogeneous materials are significantly influenced by the interface properties between the constituents. The interfaces may represent weak materials due to one or a combination of changes in mineral phase and texture, metamorphic reaction zones, and fluids. Consideration of interface properties constitutes problems of coated inclusions (Buryachenko, 2007, p.78-79). A weak interface may be due to the presence of a low-viscosity rim (Schmid and Podlachikov, 2003), creation of Somigliana’s dislocations (Mura, 1987, p.484-492) or other mechanisms of interface failure (Hashin, 1991; Huang et al., 1993; Lubarda and Markenscoff 1998). In any of these cases, the shear stress on the inclusion interface is relaxed. The simplest end-member situation is a perfectly-slippery interface which supports no shear stress. This implies that the inclusion, rigid or deformable, supports no deviatoric stresses and must behave like an inviscid inclusion (Xu et al., 2011). Eq.30c describes the instantaneous angular velocity of such an inclusion. For a perfectly-slippery inclusion, it is important to consider the shear spin effect, as in (30c), although the inclusion is ‘rigid’ in the sense that it does not change shape (because of vanishing deviatoric stresses). This is because a freely-slipping interface by definition requires that the hoop particles in contact with the inclusion slide freely (Fig.10) which is equivalent to allowing the principal axes of the inclusion to swipe away from their current positions. Failure to consider the shear spin contradicts the slippery interface assumption. There is good evidence from analogue experiments and natural observations that slippery interface may be significant in natural deformations (Ildeonse and Mancktelow, 1993; Arbaret et al., 2001; Marques and Coelho, 2001; Pennacchioni et al., 2001; Ceriani et al., 2003; Marques and Bose, 2004; Schmid and Podladchikov, 2004, 2005; Marques, et al., 2005). Eq.30c can be used to describe the angular velocity of a perfectly-slippery ellipsoid in a general anisotropic viscous material and Eq.39 is for the special situation where the matrix is isotropic. Mulchrone (2007) derived a 2D version of (39). Johnson et al. (2009) have validated the 2D equation numerically. The work of Johnson et al. (2009) confirms the importance of shear spin as discussed above.
9. Discussion 9.1 Mechanical interactions is key to natural deformation
42
ACCEPTED MANUSCRIPT
TE
D
MA
NU
SC
RI
PT
In essence, the inclusion theory is about the mechanical interactions between a heterogeneous element (an Eshelby inhomogeneity) and the embedding material (the hypothetical uniform HEM). In fact, the Eshelby-type interior solutions (12) have often been referred to as “interaction laws” (Molinari, 2002; Molinari and Mercier, 2004). To be complete, the interaction laws should also include the exterior solutions (13). Such mechanical interactions are indeed the key to understanding natural deformation and the accompanying fabric development. At the structural geology scale, the behaviors of deformable and rigid elements, which give rise to the development of shape fabrics, are manifestation of the mechanical interaction of the elements with the embedding material. The resulting interior fields in an element are responsible for the development of the fabrics inside the element. The development of tails and other geometric patterns around a heterogeneous element (e.g., Passchier and Trouw, 1996) results from the exterior flow fields and can be investigated using exterior solutions (13). The kinematics and mechanics in a Ramsay and Graham (1970) type ductile shear zone may be treated as interior problems of a flattened weak ellipsoidal inclusion (e.g., Xiang and Jiang, 2013) embedded in a country rock mass. The mechanical interactions governed by the interaction laws (12) and (13) determine the deformation history of the shear zone (interior) as well as the manner in which the active zone may increase or decrease in volume (White et al., 1980) during deformation. It is not appropriate to separate the deformation within the zone from that of the “country rock”.
AC CE P
The significance of mechanical interactions extends to all scales of natural deformation. A regional geological unit such as a plutonic body or a continental block may be regarded as a rheologically heterogeneous inclusion in a broader orogenic belt undergoing large-scale deformation modulated by plate tectonics. The deformation and related fabric development within the body constitute an interior problem. The deformation surrounding the body forms an exterior problem. In this context, one considers the orogenic system as a whole subjected to a far-field loading or displacement condition rather than treating the deformation of a geological body in isolation with ad hoc imposed boundary conditions. In such a whole system approach, the deformation of an individual element is the consequence of partitioned deformation due to rheological interactions between the element and the whole deforming system. The partitioning varies with time as the configuration of the deforming system varies continually. To capture this, partitioning (interaction laws) and homogenization equations must be considered simultaneously in a self-consistent way (Jiang 2014). Treating the deformation of an element in isolation, such as a ductile shear zone under a prescribed constant boundary condition, fails to consider the dynamic evolution of the whole deforming system and mechanical interactions among rheologically distinct geological elements. 9.2 A full mechanical and multiscale approach to natural deformation Studies of natural deformation involving large finite strains have been an important subject in structural geology and tectonics (e.g., Ramsay 1967; Hobbs et al., 1976; Passchier and Trouw, 1996). Traditionally the focus has been based on kinematic analysis. In short, the 43
ACCEPTED MANUSCRIPT
MA
NU
SC
RI
PT
kinematic approach investigates how lines, planes, and 3D elements rotate and change shape in a given flow field and, how all this relates to fabric development. The choice of the flow field is left to the intuition of the investigator. While on the plate scale, one may reasonably define the movement (the velocity field) across plate boundary regions because relative plate motions can be steady for many million years (e.g., DeMets and Wilson, 1997; DeMets and Traylen, 2000; Bird, 2003), how the encompassing rheologically heterogeneous lithosphere (made of multiscale elements) actually accommodates the imposed boundary motion, on progressively smaller scales, is governed by mechanical interactions among the constituent elements. It is impractical for one to assign a realistic flow field for the deformation within a rheological element as the latter is a partitioned field due to the interactions. The partitioning varies with time because of the continual change of the configurations of the elements and hence their mechanical interactions. The inclusion theory and associated micromechanics provide a powerful and rigorous mean to tackle such mechanical interactions and continual partitioning. It allows the kinematic aspects (strain rates, vorticity, finite strains) and dynamics of deformation (deviatoric stresses and pressure) to be derived completely from the solutions (12) and (13) rather than from simplistic kinematic models with ad hoc boundary conditions.
AC CE P
TE
D
Furthermore, the partitioning and homogenization equations are tools to bridge the scale gaps. The inclusion theory can be applied on multiple levels. Lower-order (usually larger) inclusions can contain higher-order (usually smaller in scale) inclusions. The partitioned mechanical fields in a lower-order inclusion can be further partitioned to give mechanical fields in higher-order inclusions contained in the lower-order one. This approach can thus address multi-scale deformations (Jiang and Bentley, 2012; Jiang 2014). In the self-consistent MultiOrder Power-Law Approach (MOPLA, Jiang and Bentley 2012; Jiang 2014), two levels of partitioning were considered. A macroscale deformation is first partitioned self-consistently into the constituent rheologically distinct elements making up the macroscale deforming rock masses. The partitioned flow fields in individual heterogeneous elements are used to investigate fabric development, on a smaller scale, in the elements. Through this approach, observations made on smaller scales are related to deformation boundary conditions at larger scales.
10. Summary and conclusions The long history, continued interest, and recent advances in inclusion studies are testament of the enormous significance of the problem. The theory of viscous inclusions is the backbone for the visco-plastic self-consistent model for texture development in polycrystalline materials (Lebensohn and Tomé, 1993), the self-consistent Multi-Order Power Law Approach (MOPLA) for ductile deformation of Earth’s heterogeneous lithosphere (Jiang, 2014), and many self-consistent homogenization models to estimate the overall rheological properties of heterogeneous materials (Molinari, 2002). The considerable literature over many disciplinary fields on the problem is, however, a challenge to earth scientists. As pointed out in this paper, 44
ACCEPTED MANUSCRIPT
RI
PT
duplicate investigations have occurred. In an effort to summarize and advance the current status on viscous inclusion studies, the powerful Green function method has been used in this paper to unify the mathematic notations and presentation. In order to address viscous incompressibility explicitly, two Green functions are used, following Molinari et al. (1987). This necessitates deriving expressions for 3D Green functions and their derivatives for incompressible materials using Barnett’s (1973) method. More explicit solutions than Molinari et al. (1987) are obtained for both the interior fields, expressed in (12), and exterior-point fields, expressed in (13). The general solutions are expressed in terms of Eshelby tensors and an auxiliary tensor
NU
SC
for the pressure field (the interior S , Π , Λ and exterior S E , Π E , Λ E ). It is straightforward to develop numerical algorithms to evaluate these quantities using standard mathematics applications. For isotropic viscous materials, evaluations of all these tensors are simpler, more
MA
efficient and more accurate with the advances made in this study. The evaluations of S E and Π E remain computationally intensive for general anisotropic materials. Fortunately, for the selfconsistent MOPLA and homogenization models, only interior solutions are needed.
TE
D
For isotropic viscous inclusion problems, many mathematic identities among viscous tensor quantities are discovered in this paper that greatly simplify the solutions and their numerical evaluations. These identities enable the use of existing analytic expressions for isotropic elastic interior-point Eshelby tensors to calculate the interior-point viscous Eshelby tensors ( S , Π ) and Λ . The identities also enable the use of the explicit method of Ju and Sun
AC CE P
(1999) for exterior-point elastic Eshelby tensor to calculate viscous exterior-point S E , Π E , and Λ E . As a result, numerical solutions for isotropic inclusion problems can be obtained much more efficiently and accurately than before. The most general equation is established for the rotation of a deformable anisotropic 1 kl wijs , Eq.28). ellipsoid in a rheologically dissimilar anisotropic matrix ( ij Wij ijmn Smnkl 1 The particularization of this equation for rigid ellipsoids leads to ( ij Wij ijmn Smnkl Ekl , Eq.29)
which is the generalization of Jeffery’s (1922) equation to an arbitrary anisotropic viscous matrix. The rotational behavior of a slippery ellipsoid is also obtained 1 a 2 aJ2 d d [ ij Wij ijkl I2 J ijkl J S mn , Eq.30c]. 2 klmn aI aJ
The shear spin of a deformable ellipsoid describes the swiping-through-material phenomenon of the ellipsoid’s principal axes. This phenomenon has nothing to do with material rotation of the ellipsoid which is described by the vorticity in the ellipsoid (the partitioned field). This helps to clarify and justify Jiang’s (2012) treatment for shear spin in the singular event where a deformable ellipsoid passes through a spheroid or sphere state.
45
ACCEPTED MANUSCRIPT
PT
In applying the general solutions (12) and (13) to isotropic inclusions, it is shown that all existing results that have been developed using different approaches in the literature can be derived directly from the general formulation without undue difficulty. This again simplifies inclusion research. Isotropic solutions can serve as benchmarks for numerical anisotropic solutions, and together they can be used to benchmarks more complex numerical models.
MA
NU
SC
RI
There has been significant progress on 2D inclusion studies in anisotropic elasticity. The general complete solutions are now available for generalized plane deformation (Ting, 1996). This body of work adopts the complex variable formalism. Because of the nature of the formalism and resulting generally implicit expressions, it remains open how or whether the results can be applied to incompressible viscous anisotropic materials. Considering this and the fact that numerical evaluation of 3D inclusion problems in anisotropic materials is computationally intensive, a 2D formulation is developed, using the same Green function method that was used in 3D, specifically for anisotropic viscous elliptic inclusions in a rheologically dissimilar anisotropic viscous matrix that can host plane-straining flows. Numerical calculations based on this 2D formulation agree perfectly with the analytic solutions of Fletcher (2009).
AC CE P
TE
D
A reconnaissance investigation on elliptic cylinders and ellipsoids in a matrix of planar anisotropy undergoing plane-straining flows shows that elliptic cylinders behave quite differently from 3D ellipsoids. This implies that conclusions drawn from 2D inclusion studies may not apply to 3D situations for anisotropic materials. As rocks are generally anisotropic, one should always consider the 3D effect of deformations where heterogeneous inclusions are investigated. More detailed investigation on this is under way. The application of the viscous inclusion theory based on linear rheology to the rheologically non-linear lithosphere is through the linearization approach. The heterogeneous lithosphere is approximated by a homogeneous equivalent matrix whose effective rheology is obtained through a micromechanical homogenization method. A heterogeneous rheological element is regarded as an ellipsoidal inclusion which interacts with the homogeneous equivalent matrix. The conditions of scale separation must also be satisfied in order for a finite volume deformation problem to be regarded as one in an infinite body. The interface between the inclusion and the matrix is assumed to be coherent. These conditions and assumptions may be justified from certain problems under given circumstances. Where they cannot be justified, modification of the theory is necessary to account for, such as irregular inclusion shapes, incoherent interface, and boundary effects of deformation. The viscous inclusion theory with its extension to rheologically heterogeneous and nonlinear materials forms the backbone theory to address the ubiquitous mechanical interactions on all scales in natural rock deformations. It enables a full mechanical and multiscale approach toward lithospheric deformation and the accompanying texture and fabric development.
46
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Acknowledgements: The author thanks the Journal editor and three reviewers for their critical comments and acknowledges financial support for the research from Natural Science and Engineering Research Council of Canada (NSERC) through a Discovery Grant, the China National Natural Science Foundation (grant NSFC41472184), a Distinguished Visiting Professorship from China Shaanxi Province and Northwest University as well as internal funding from the State Key Laboratory for Continental Dynamics at Northwest University.
47
ACCEPTED MANUSCRIPT Appendices A: Expressions of Gij , Gij ,l , and H i
PT
Integral expressions for Gij and Gij ,l for anisotropic elastic materials are given by
RI
Barnett (1972) and Mura (1987, p.33-34). Barnet’s derivations are repeated here in order to get expressions of Gij , Gij ,l , and H i for incompressible viscous materials.
SC
The two Greens functions in real space are inverse Fourier transforms of the expressions in (6): Aij (z) exp iK x dK 8 K 2
Hi x
i i exp iK x dK 8 K
3
3
MA
1
(A1)
(A2)
D
K and the integrations are over infinite 3D space. K
TE
where z
1
NU
Gij x
With Euler’s formula, the expressions become:
1
Aij (z) cos Kz x dK 8 K 2
AC CE P
Gij x
Hi x
3
1
8
3
(A3)
i
K sin Kz x dK
(A4)
Here only the real parts of the integrals are concerned as Gij and H i are real. The expression for Gij , l can be obtained by differentiating (A3) with respect to xl :
Gij ,l x
1
zl Aij (z) sin Kz x dK 8 K 3
(A5)
With the following set of variable substitution: x
x k dk , K , and therefore dK 3 r r r
(A6)
(A3-5) are reduced to:
48
ACCEPTED MANUSCRIPT
1 8 r
i
3 2
Gij ,l x
Aij (z) cos kz x dk k2
r
k
sin kz x dk
1
(A7a)
PT
8
3
zl Aij (z) sin kz x dk 8 r k
3 2
RI
Hi x
1
SC
Gij x
(A7b)
(A7c)
A spherical coordinate system aligned with x can be used to evaluate these integrals:
NU
d 3k k 2 sin dkd d , z x cos
(A8)
MA
where is a polar angle in the plane perpendicular to x (i.e., z x 0 ). With (A8), the three integrals in (A7) become:
Gij x 3 Aij (z)sin cos k cos dk d d 8 r 0 0 0
(A9a)
H i x 3 2 i sin k sin k cos dk d d 8 r 0 0 0
(A9b)
Gij ,l x 3 2 zl Aij (z)sin k sin k cos dk d d 8 r 0 0 0
(A9c)
2
TE
D
1
2
AC CE P
1
1
2
Making use of the following identities (e.g., Jones, 1966): cos k cos dk cos , 0
k sin k cos dk sin cos , and
f ( x ) ( n ) x dx ( 1)n f ( n ) (0) , we get from (A9):
0
Gij x
2
1 8 2 r
Hi x
A (z ) ij
0
1 8 r
2 2
2
i
0
d
(A10a)
2
d
(A10b)
2
49
ACCEPTED MANUSCRIPT
z A ij
l
0
More explicit expressions for
i
d
(A10c)
2
and
zl Aij
2
2
SC
First, z x cos implies (Barnett, 1972):
zi xi
2
Aij Aij z l Aij zl x l Aij zl
(A12)
MA
zl Aij
(A11)
NU
2
thus
are derived below.
PT
8 r
2 2
2
RI
Gij ,l x
1
2
D
i The expressions for Aij and are obtained by differentiating the identity
AC CE P
TE
A z A ζ I and reorganizing: T T z 0ζ
Aij A Aim A jn mn Aim j A jm i x m
(A13a)
i A Ain mn m Aim i m x m
(A13b)
Using the definition of the Christoffel stiffness matrix in (5a), we have A
2
z p zq C p q z p zq C p q z p x q zq x p
(A14)
2
Calling C p q z p x q zq x p Fi
A
and inserting it into (A13), we get:
i Ain mn m Aim i m x m
(A15a)
50
ACCEPTED MANUSCRIPT M ij
Aij Aim A jn mn Aim j A jm i x m
(A15b)
8 2 r 2
F d i
0
2
1 8 2r 2
x A
ij
l
zl M ij d
0
RI
Gij ,l (x)
2
1
SC
H i ( x)
PT
With (A15), (A10b) and (A10c) become respectively: (A16a)
(A16b)
NU
Finally, from (5a), it is clear that Aij z Aij z . From (7) it is easy to show that the following relations are true:
MA
1 1 zi i 1 , Aij Aim jm j zm Ajm im i zm , and Aij i z j
(A17)
D
Therefore, we have Aij (z) Aij (z) , (z) (z) , and i (z) i (z) . Considering this in
TE
(A15), it is clear that the integrands in (A10a) and (A16) are all even with respect to z. Thus Gij ,
AC CE P
Gij , l , and H i can be written as:
Gij (x)
1
4 2 r 0
Aij d
(A18a)
1 H i (x) 2 2 Fd i 4 r 0 Gij ,l (x)
1
4 2r 2 0
x A l
ij
(A18b)
zl M ij d
(A18c)
For isotropic materials, (A17) becomes:
i zi , Aij
1 z z , and ij i j
(A19)
B: Relations among tensor quantities for isotropic incompressible materials 51
ACCEPTED MANUSCRIPT Consider the integral of (14b) for isotropic incompressible viscous materials with
1
PT
s Cijkl 2 J ijkl , where is the viscosity. It can be shown that Aik ik zi zk , which turns
(14b) to:
ik
zi zk z j zl 3dS
RI
a1a2 a3 4
z 1
(B1)
SC
Tijkl
2 ij kl , ( and being respectively 1 2
NU
where is the same as defined in (14).
MA
el s For isotropic elastic materials, Cijkl 2 J ijkl +
the shear modulus and Poisson’s ratio, e.g., Mura 1987). We have Aij ij
el Tijkl
a1a2 a3 4
TE
D
1 1 zi z j . The interaction tensor is: ij 2 1 1 zi zk z j zl 3dS ik 2 1 z 1
AC CE P
Aij1
1 zi z j and 1 2
s The symmetric part of (B1) and (B2), Tijkl
(s) Tijkl
a1a2 a3 4
el(s) Tijkl
a1a2 a3 4
(B2)
1 Tijkl Tjikl Tijlk Tjilk , are respectively: 4
1 3 4 ik z j zl jk zi zl il z j zk jl zi zk zi z j zk zl dS z 1
(B3a)
1 1 zi z j zk zl 3dS (B3b) ik z j zl jk zi zl il z j zk jl zi zk 2 1 z 1 4
Therefore, we have from (B3):
Tijkl( s ) Tijklel(s)
(B4)
=0.5
The viscous Eshelby tensor and the elastic Eshelby tensor are respectively:
52
ACCEPTED MANUSCRIPT (s) Sijkl 2Tijkl
(B5) 2 el(s) Tij kl 1 2
(B6)
PT
el el(s) Sijkl 2Tijkl
Tijel(s)
a1a2 a3 1 2
zz 8 1 i
3
j
dS
1
el_s ij kl 2Tijkl
which yields: 0.5
el(s) ij kl 2 Tijkl
0.5
(B8)
D
el Sijkl
MA
NU
Submitting (B7) to (B6) reduces it to: el Sijkl
(B7)
SC
ξ 1
1 2 ij 2 1
RI
From (B3b) and (11), one gets:
TE
Considering (B4) and replacing the shear modulus by viscosity in (B8), the second term on the
AC CE P
right of (B8) becomes the viscous Sijkl because of (B5) and we have: el Sijkl ij kl Sijkl
(B9)
0.5
el Furthermore, because Sijkl are zero except for terms with repeated indices and because
Sii11 Sii 22 Sii 33 0 due to incompressibility (Jiang 2013), we get the following two relations
from (B9): ii
1 el Sii11 Siiel22 Siiel33 0.5 3
el Siijj ii Siijj
0.5
(no sum)
(no sum)
(B10) (B11)
which are respectively (17a) and (17b) in the text.
53
ACCEPTED MANUSCRIPT a Taking the anti-symmetric part of (B1) and (B2), Tijkl
1 Tijkl Tjikl Tijlk Tjilk , and 4
a1a2 a3 16
ik
z j zl jk zi zl il z j zk jl zi zk 3dS
z 1
RI
(a) ( a _ el ) Tijkl Tijkl
PT
a then contract Tijkl with the corresponding stiffness, it is easy to demonstrate that
ik
z j zl jk zi zl il z j zk jl zi zk 3dS
z 1
(B12)
NU
a1a2 a3 8
( el ) ijkl ijkl
SC
which is equivalent to:
MA
This is (18) in the text.
C: Derivation of equations for incompressible isotropic materials
ai2 a 2j ai2 a 2j
S ijij
(no sum)
(C1)
TE
ijij
D
First, we prove that the following relation holds for isotropic materials:
Siiii
AC CE P
Proof: Components of Eshelby tensors for isotropic incompressible materials are zero except for the following ones with repeated indices [combining Eshelby, 1957, (17) and (18) in the text]: 3 2 ai 4
ii
ii , Siijj
ijij ijji
j
8
i
3 2 aj 4
ij
ii , Sijij
3 ai 2 a j 2 8
, jiij jiji ijij
ij
,
(C2a)
(C2b)
where, as in the rest of this Appendix, repeated indices do not imply summation unless specified otherwise and the
terms are defined by the following integrals (Jeffery, 1922; Eshelby, 1957):
i
du 2 a1a2 a3 2 , 0 ai u
ii
2 a1a2 a3 0
a
2 i
du
u 2
, and
ij
2 du a1a2 a3 2 , ( i j) 2 3 0 ai u a j u
54
(C3)
ACCEPTED MANUSCRIPT a
where
2 1
u a22 u a32 u .
1 j i , ( i j , ai a j ) 3 ai2 a 2j
SC
Combining (C4) and the 3rd expression in (C2a) gives (C1).
(C4)
RI
ij
PT
The following relation was given by Eshelby (1957, Eq.3.13 there) but also follows directly from C3:
1 4Sijij
MA
1 Sijij
NU
Second, we prove the following relation which will be used in the derivation of many expressions:
d Proof: Because there are only 4 non-vanishing components of J ijkl if
(C5)
i j and they are
sum on lowercase indices only
TE
1 d SIJkl S klIJ J IJIJ 0.5
D
d d J ijij J djiij J ijji J djiji 0.5 . We therefore have by definition:
(C6)
AC CE P
Spelling out all nonzero terms in the contraction, the left side of (C6) is: 1 1 1 1 SIJkl SklIJ SIJIJ SIJIJ SIJJI S JIIJ 2SIJIJ SIJIJ
(no sum)
(C7)
Combining (C6) and (C7) gives (C5). Third, we derive (32) in the text. The tensor form equation (31) is rewritten below omitting the label: 1
ε J d (r 1)S : Ε
(C8)
For the three diagonal components of the strain rates ii (i 1, 2,3) , we expand (C8):
II (r 1)SIIkl kl II
(sum on repeated lowercase indices only)
Expanding (C9) and using the incompressible condition lead to the following linear set:
55
(C9)
ACCEPTED MANUSCRIPT 11 (r 1) S111111 S1122 22 S1133 33 11 22 (r 1) S221111 S2222 22 S2233 33 22
(C10)
PT
33 (r 1) S331111 S3322 22 S3333 33 33 11 22 33 11 22 33 0
RI
which can be reduced to two equations with two unknowns ( 11 and 22 ) written in the following
SC
matrix form:
(C11)
NU
1 (r 1) S1111 S1133 (r 1) S1122 S1133 11 11 (r 1) S2211 S2233 1 (r 1) S2222 S2233 22 22
MA
The solution to this system gives the following expressions (see Freeman, 1987 with errors corrected in Jiang 2007):
(C12)
AC CE P
TE
D
1 1 (r 1) S2222 S2233 11 (r 1) S1122 S1133 22 1 22 (r 1) S2233 S2211 11 1 (r 1) S1111 S1133 22 33 11 22
11
where 1 (r 1) S1111 S1133 1 (r 1) S2222 S2233 (r 1)2 S1122 S1133 S2211 S2233 . This is (32a) in the text.
To get (32b) for the off-diagonal terms of ij ( i j ), we proceed as follows. First by denoting Ξ J d (r 1)S , (C8) is rewritten as:
ε Ξ1 : Ε
(C13)
It follows readily that:
ijij J ijijd (r 1)Sijij 0.5 (r 1)Sijij , i j
(C14)
Applying the same argument and procedure that (C5) was proven to (C14), we get: 1 ijij
1 4 0.5 (r 1) Sijij
(C15)
Expanding (C13) and making use of (C15), we get: 56
ACCEPTED MANUSCRIPT 1 1 1 1 IJ IJkl kl IJIJ IJ IJJI JI 2 IJIJ IJ
IJ , IJ 1 2(r 1) S IJIJ
(C16)
PT
which is (32b) in the text. Fourth, we derive (33) below.
IJIJ
1 aI2 aJ2 4 aI2 aJ2
SC
S
1 IJIJ
RI
Combining (C1) and (C5) leads to the following relation:
(C18)
MA
1 wij ijkl Sklmn mn
NU
Now, write (12b) in index form:
(C17)
TE
to the kl indices, (C18) becomes:
D
Because there are only 2 nonzero terms for ijkl (i.e., ijij ijji ), after contraction with respect
AC CE P
1 wIJ 2 IJIJ SIJmn mn (sum on repeated lowercase indices only)
(C19)
Because there are only 2 nonzero terms of S IJmn for isotropic materials ( S IJIJ S IJJI ), after contraction with respect to the mn indices (C19) reduces to: 1 wIJ 4 IJIJ SIJIJ IJ
(no sum)
(C20)
Applying (C17) to (C20) yields: w IJ
aI2 aJ2 IJ aI2 aJ2
(C21)
which is (33) in the text. D: Orthonormal base tensors for plane-straining analysis To perform inversion of fourth-order symmetric tensors in 2D space, we can follow the approach of Lebensohn et al. (1998) and Jiang (2013) by introducing a set of 3 orthonormal 2nd order tensors: 57
ACCEPTED MANUSCRIPT b1
1 1 0 1 0 1 1 1 0 2 3 , b , b 2 0 1 2 1 0 2 0 1
(D1)
This set of tensors are orthonormal, satisfying b : b . A symmetric fourth-order tensor in
PT
2D, such as the symmetric Eshelby tensor S, can be reshaped into a 3 by 3 matrix by: ( , 1..3 )
(D2)
RI
Sλη = bij bkl Sijkl
λη
This inversion ensures that 1 ik jl jkil 2
MA
1 1 Sijmn Smnkl Sijmn Smnkl
(D3)
NU
-1
1 Sijkl bij bkl S
SC
The matrix can be inverted and then reshaped back to a fourth-order tensor through:
i, j, k , l 1,2
(D4)
AC CE P
TE
D
For incompressible materials in 2D, the deviatoric stress and strain rate tensors only have 2 independent components. The corresponding 4th order stiffness tensor is singular in the complete Cauchy stress space. The inversion of such fourth-order tensors can be constrained in the deviatoric space with , 1,2 (compare with 3D case illustrated in Jiang 2013). Such inversion ensures that: 1 1 2D Sijmn Smnkl Sijmn Smnkl J ijkl
1 1 ik jl jk il ij kl 2 2
58
i, j, k , l 1,2
(D5)
ACCEPTED MANUSCRIPT REFERENCES Arbaret, L., Mancktelow, N.S., Burg, J.P., 2001. Effect of shape and orientation on rigid particle rotation and matrix deformation in simple shear flow. Journal of Structural Geology 23,
PT
113-125.
RI
Astarita, G. 1979. Objective and generally applicable criteria for flow classification, J. NonNewtonian Fluid Mechanics 6, 69–76.
SC
Barnett, D.M., 1972. The precise evaluation of derivatives of the anisotropic elastic Green’s functions. Physica Status Solidi (b), 49, 741-748.
NU
Başar, Y., and Weichert, D., 2000. Nonlinear Continuum Mechanics of Solids. Springer. Bhargava, R.D., Radhakrishna, H..C. 1964. Elliptic inclusion in orthotropic medium. Journal of
MA
the physical society of Japan, 19(3), 396-405.
Bilby, B.A., Eshelby, J.D., and Kundu, A.K., 1975, The change of shape of a viscous ellipsoidal region embedded in a slowly deforming matrix having a different viscosity:
D
Tectonophysics, v. 28, p. 265-274.
TE
Bilby, B.A., Kolbuszewski, M.L., 1977. The finite deformation of an inhomogeneity in twodimensional slow viscous incompressible flow. Proceedings of the Royal Society of
AC CE P
London, A355, 335-353.
Bird, P. 2003. An updated digital model of plate boundaries. Geochemistry Geophysics Geosystems, 4(3), 1027, doi:10.1029/2001GC000252., Brace, W.F., 1965. Relation of elastic properties of rocks to fabric. Journal of Geophysical Research, 70, 22-35.
Bretherton, F.P., 1962. The motion of rigid particles in a shear flow at low Reynolds number. Journal of Fluid Mechanics, 14, 284-304. Budiansky, B., Mangasarian, O., 1960. Plastic stress concentration at circular hole in infinite sheet subjected to equal biaxial tension. Journal of Applied Mechanics, 27, 59-64. Bunge, H.J. 1982. Texture analysis in materials science, mathematical methods. (English edition translated by Morris, P.A.), Butterworths & Co. Buryachenko, V.A., 2007. Micromechanics of heterogeneous materials. Springer. Ceriani, S., Mancktelow, N.S., Pennacchioni, G., 2003. Analogue modelling of the influence of shape and particle/matrix interface lubrication on the rotational behaviour of rigid particles in simple shear. Journal of Structural Geology 25, 2005-2021. 59
ACCEPTED MANUSCRIPT Chen, T., 1967. On an elliptic elastic inclusion in an anisotropic medium. Quarterly Journal of Mechanics and Applied Mathematics. 20, 307-313. Chen, Y., Jiang, D., Zhu, G., Xiang, B. 2014. The formation of micafish: a modeling
PT
investigation based on micromechanics. Journal of Structural Geology 68, 300-315. Cook, A.C., Vel, S.S., Gerbi, C., Johnson, S.E., 2014. Computational analysis of nonlinear creep
RI
of polyphaser aggregates: Influence of phase morphology. Journal of Geophysical
SC
Research: Solid Earth, 119, 6877-6906. doi:10.1002/2014JB011197. Davis, J.R., Titus, S.J., Horsman, E., 2013. Non-steady homogeneous deformations:
NU
Computational techniques using Lie theory, and application to ellipsoidal markers in naturally deformed rocks. Journal of Structural Geology, 56, 142-155.
MA
Dabrowski, M., Schmid, D.W. 2011. A rigid circular inclusion in an anisotropic host subject to simple shear. Journal of Structural Geology 33 (7), 1169-1177.
D
Dabrowski, M., Schmid, D.W., Podladchikov, Y.Y. 2012. A two-phase composite in simple
TE
shear: effective mechanical anisotropy development and localization potential. Journal of Geophysical Research, 117, B08406, doi:10.1029/2012JB009183.
AC CE P
DeMets, C., Traylen, S., 2000. Motion of the Rivera plate since 10 Ma relative to the Pacific and North American plates and the mantle. Tectonophyics, 318, 119–159. DeMets, C., Wilson, D.S., 1997. Relative motions of the Pacific, Rivera, North America, and Cocos plates since 0.78 Ma. Journal of Geophysical Research, 102, 2789–2807. Eshelby, J.D., 1957, The determination of the elastic field of an ellipsoidal inclusion, and related problems: Proceedings of the Royal Society of London, v. A241, p. 376-396. Eshelby, J.D., 1959. The elastic field outside an ellipsoidal inclusion. Proceedings of the Royal society of London, A 252, 561-569. Eshelby, J.D., 1961. Elastic inclusions and inhomogeneities. In: Snedden, N.I., Hill, R. (eds.) Progress in solid mechanics, Vol. 2, 89-140. Eshelby, J. D., Read, W.T., Shockley, W. 1953. Anisotropic elasticity with applications to dislocation theory. Acta Metallurgica. 1, 251-259. Exner, U. and Dabrowski, M., 2010. Monoclinic and triclinic 3D flanking structures around elliptical cracks. Journal of Structural Geology, 32(12): 2009-2021.
60
ACCEPTED MANUSCRIPT Ferguson, C.C., 1979. Rotations of elongate rod of particles in slow non-Newtonian flows. Tectonophysics 60, 247-262. Fletcher, R.C. 2004. Anisotropic viscosity of a dispersion of aligned elliptical inclusions in a
PT
homogeneous incompressible anisotropic viscous fluid. Journal of Structural Geology, 26, 1977-1987.
RI
Fletcher, R.C. 2009. Deformable, rigid, and inviscid elliptical inclusions in a homogeneous
SC
incompressible anisotropic viscous fluid. Journal of Structural Geology, 31, 382-387. Fletcher, R.C., Pollard, D.D. 1999. Can we understand structural and tectonic processes and their
NU
products without a complete mechanics? Journal of Structural Geology, 21, 1071-1088. Freeman, B. 1985. The motion of rigid ellipsoidal particles in three-dimensional slow flows:
MA
implications for geological strain analysis. Tectonophysics 113, 163-183. Freeman, B. 1987. The behavior of deformable ellipsoidal particles in three-dimensional slow flows: implications for geological strain analysis. Tectonophysics, 132, 297-309.
D
Gao, X.-L., Ma, H.M., 2010. Strain gradient solution for Eshelby’s ellipsoidal inclusion problem.
TE
Proceedoings of the Royal Society A, 466, 2425–2446. doi:10.1098/rspa.2009.0631 Gay, N.C., 1968. The motion of rigid particles embedded in a viscous fluid during pure shear
AC CE P
deformation of the fluid. Tectonophysics 5, 81-88. Ghosh, S.K., Ramberg, H., 1976. Reorientation of inclusions by combination of pure and simple shear. Tectonophysics, 34, 1-70. Goddard, J.D., Miller, C., 1967. Nonlinear effects in the rheology of dilute suspensions. Journal of Fluid Mechanics 28, 657-673. Goldsmith, H.L., Mason, S.G., 1967. The microrheology of dispersions. Rheology: theory and applications. Vol 4, 86-250. Griera A., Bons P.D., Jessell M.W., Lebensohn R., Evans L. and Gomez-Rivas E. 2011. Strain localization and porphyroclast rotation. Geology, 39, 275-278. Griera, A., Llorens, M-.G., Gomez-Rivas, E., Bons, P.D., Jessell, M.W., Evans, L.A. and Lebensohn, R. 2013. Numerical modelling of porphyroclast and porphyroblast rotation in anisotropic rocks. Tectonophysics, 587, 4–29.
61
ACCEPTED MANUSCRIPT Hashin, Z., 1991. Thermoelastic properties of particular composites with imperfect interface. Journal of Mechanics and Physics of Solids, 39, 745–762. Hill, R. 1961. Discontinuity relations in mechanics of solids. In: Snedden, N.I., Hill, R. (eds.)
PT
Progress in solid mechanics, Vol. 2, 245-276.
Hill, R. 1983. Interfacial operators in the mechanics of composite media. Journal of Mechanics
RI
and Physics of Solids, 31, 347–357
SC
Hobbs, B.E., Means, W.D., and Williams, P.F., 1976. An Outline of Structural Geology. John Wiley & Sons.
NU
Hobbs, B.E., Ord, A., 2015. Structural geology: the mechanics of deforming metamorphic rocks. Volume 1: principles. Elsevier.
MA
Huang, J.H., Furuhashi, R., Mura, T., 1993. Frictional sliding inclusions. Journal of Mechanics and Physics of Solids, 41, 247–265.
Hutchinson, J.W., 1976. Bounds and Self-Consistent Estimates for Creep of Polycrystalline
D
Materials. Proceedings of the Royal Society of London. A. Mathematical and Physical
TE
Sciences, 348(1652): 101-127.
Hwa, C., Ting, T.C.T., 1989. Two-dimensional problems of the anisotropic elastic solid with an
AC CE P
elliptic inclusion. Quarterly Journal of Mechanics and Applied Mathematics, 42, 553-572. Ildefonse, B., Mancktelow, N.S., 1993. Deformation around rigid particles: the influence of slip at the particle/matrix interface. Tectonophysics 221, 345-359. Ildefonse, B., Launeau, P., Fernandez, A., Bouchez, J.L., 1992. Effects of mechanical interactions on development of shape preferred orientations: a two-dimensional experimental approach. Journal of Structural Geology 14, 73-83. Jaswon, M.A., Bhargava, R.D., 1961. Two-dimensional elastic inclusion problems. Proc. Camb. Phil. Soc., 57, 669-680. Jeffery, G.B., 1922, The motion of ellipsoidal particles immersed in a viscous fluid: Proceedings of the Royal Society of London, v. A102, p. 161-179. Ježek, J., Melka, R., Schulmann, K., Venera, Z., 1994. The behaviour of rigid triaxial ellipsoidal particles in viscous flows -- modeling of fabric evolution in a multi-particle system. Tectonophysics 229, 165-180. Ježek, J., Schulmann, K., Segeth, K., 1996. Fabric evolution of rigid inclusions during mixed coaxial and simple shear flows. Tectonophysics 257, 203-221. 62
ACCEPTED MANUSCRIPT Jiang, D. 1999. Vorticity decomposition and its application to sectional flow characterization. Tectonophysics 301, 243-259. Jiang, D., 2007a. Numerical modeling of the motion of rigid ellipsoidal objects in slow viscous
PT
flows: A new approach. Journal of Structural Geology, 29(2): 189-200. Jiang, D., 2007b. Numerical modeling of the motion of deformable ellipsoidal objects in slow
RI
viscous flows. Journal of Structural Geology, 29(3): 435-452.
SC
Jiang, D., 2010. Flow and finite deformation of surface elements in three dimensional homogeneous progressive deformations. Tectonophysics, 487(1-4): 85-99.
NU
Jiang, D., 2012. A general approach for modeling the motion of rigid and deformable objects in ductile flows. Computers & Geosciences 38, 52-61.
MA
Jiang, D., 2013. The motion of deformable ellipsoids in power-law viscous materials: Formulation and numerical implementation of a micromechanical approach applicable to flow partitioning and heterogeneous deformation in Earth’s lithosphere. Journal of
D
Structural Geology 50, 22-34.
TE
Jiang, D. 2014. Structural geology meets micromechanics: a self-consistent model for the multiscale deformation and fabric development in Earth’s ductile lithosphere. Journal of
AC CE P
Structural Geology 68, 247-272.
Jiang, D. and Bentley, C., 2012. A micromechanical approach for simulating multiscale fabrics in large-scale high-strain zones: Theory and application. Journal of Geophysical Research 117: B12201, doi:12210.11029/12012JB009327. Jiang, D., Williams, P.F., 1998, High-strain zones: a unified model. Journal of Structural Geology, v. 20, p. 1105-1120. Jin, X., Keer, L.M., Wang, Q.J., 2011. A closed-form solution for the Eshelby tensor and the elastic field outside an elliptic cylindrical Inclusion. Journal of Applied Mechanics, Transactions ASME 78, 031009. Johnson, A.M., Fletcher, R.C., 1994. Folding of viscous layers. Columbia University Press, New York. Johnson, S.E., Lenferink, H.J., Price, N.A., Marsh, J.H., Koons, P.O., West Jr., D.P., Beane, R., 2009. Clast-based kinematic vorticity gauges: the effects of slip at matrix/clast interfaces. Journal of Structural Geology 31, 1322-1339. Jones, D.S., 1966. Generalized functions. McGraw-Hill. 63
ACCEPTED MANUSCRIPT Ju, J.W., Sun, L.Z., 1999. A novel formulation for the exterior point Eshelby’s tensor of an ellipsoidal inclusion. Transactions of the ASME, 66, 570-574. Kang, H., Milton, G.W., 2008. Solution to the polya-sazego conjecture and the weak Eshelby
PT
conjecture. Archival Rational Mechanics and Analysis, 188, 93-116. 2008 Kinoshita, N., Mura, T., 1984. Eigenstrain problems in a finite elastic body. SIAM Journal on
RI
Applied Mathematics 44, 524–535.
SC
Kocks, U.F., 1998. Kinematics and kinetics of plasticity.in Kocks, U.F., Tomé, C.N., and Wenk, H.-R. ed. Texture and anisotropy: preferred orientations in polycrystals and their effect
NU
on materials properties. Cambridge University Press. Chapter 8, 326-389. Korn, G.A., and Korn, T.M. 1968. Mathematical handbook for scientists and engineers. 2nd
MA
edition, McGraw-Hill, New York.
Lebedev, V. I., 1977. Spherical quadrature formulas exact to orders 25–29. Siberian Mathematical Journal, 18(1), 99-107.
D
Lebensohn, R. A., Dawson, P. R., Kern, H. M., Wenk, H. R., 2003. Heterogeneous deformation
TE
and texture development in halite polycrystals: Comparison of different modeling approaches and experimental data. Tectonophysics, 370, 287-311.
AC CE P
Lebensohn, R.A., and Tomé, C.N., 1993, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys: Acta Metallurgica et Materialia, 41, 2611-2624. Lebensohn, R.A., Tomé, C.N. and Maudlin, P.J., 2004. A selfconsistent formulation for the prediction of the anisotropic behavior of viscoplastic polycrystals with voids. Journal of the Mechanics and Physics of Solids, 52(2): 249-278.Lebensohn, R.A., Turner, P.A., Signorelli, J.W., Canova, G.R. and Tomé, C.N., 1998. Calculation of intergranular stresses based on a large-strain viscoplastic self-consistent polycrystal model. Modelling and Simulation in Materials Science and Engineering, 6(4): 447-465. Lekhnitskii, S.G., 1950. Theory of elasticity of an anisotropic elastic body. Gostekhizdat, Moscow (in Russian). Holden-Day, San Francisco (in English, 1963). Li, C., Jiang, D., 2011. A critique of vorticity analysis using rigid clasts. Journal of Structural Geology 33, 203-219. Li, S., Sauer, R.A., Wang, G., 2007. The Eshelby tensors in a finite spherical domain – Part I: Theoretical formulations. Journal of Applied Mechanics 74, 770–783. 64
ACCEPTED MANUSCRIPT Li, S., and Wang, G., 2008. Introduction to micromechanics and nanomechanics. World Scientific. Lin, S., Jiang, D., and Williams, P.F. 1998. Transpression (or transtension) zones of triclinic
PT
symmetry: natural example and theoretical modelling. In: Holdsworth, R.E., Strachan, R.A. & Dewey, J.F. (eds) 1998. Continental transpressional and transtensional tectonics.
RI
Geological Society, London, Special Publications, No. 135, p. 41-57.
SC
Lister, G. and Williams, P.F., 1983. The partitioning of deformation in flowing rock masses. Tectonophysics, 92(1-3): 1-33.
NU
Liu, L., 2008. Solutions to the Eshelby conjectures. Proceedings of Royal Society of London. 464, 573-594.
MA
Lubarda, V.A., Markenscoff , X., 1998. On the stress field in sliding ellipsoidal inclusions with shear eigenstrain. Journal of Applied Mechanics, 65, 858–862. Mancktelow, N. S., 2011. Deformation of an elliptical inclusion in two-dimensional
D
incompressible power-law viscous flow. Journal of Structural Geology, 33(9): 1378-1393.
TE
Mancktelow, N. S., 2013. Behaviour of an isolated rimmed elliptical inclusion in 2D slow incompressible viscous flow. Journal of Structural Geology, 46: 235-254.
AC CE P
Mancktelow, N.S., Arbaret, L., Pennacchioni, G., 2002. Experimental observations on the effect of interface slip on rotation and stabilisation of rigid particles in simple shear and a comparison with natural mylonites. Journal of Structural Geology 24, 567-585. Mandal, N., Misra, S., Samanta, S.K., 2005a. Rotation of single rigid inclusions embedded in an anisotropic matrix: a theoretical study. Journal of Structural Geology 27, 731-743. Mandal, N., Samanta, S.K., Bhattacharyya, G., Chakraborty, C., 2003. Deformation of ductile inclusions in a multiple inclusion system in pure shear. Journal of Structural Geology 25, 1359-1370. Mandal, N., Samanta, S.K., Chakraborty, C., 2005b. Numerical models of flow patterns around a rigid inclusion in a viscous matrix undergoing simple shear: implications of model parameters and boundary conditions. Journal of Structural Geology 27, 1599-1609. Markenscoff, X., Gupta, A., 2006. Collected Works of J.D. Eshelby: the mechanics of defects and inhomogeneities. Springer, Dordrecht, Netherlands.
65
ACCEPTED MANUSCRIPT Marques, F.O., Bose, S., 2004. Influence of a permanent low-friction boundary on rotation and flow in rigid inclusion/viscous matrix systems from an analogue perspective. Tectonophysics 382, 229-245.
PT
Marques, F.O., Coelho, S., 2001. Rotation of rigid elliptical cylinders in viscous simple shear flow: analogue experiments. Journal of Structural Geology 23, 609-617.
RI
Marques, F.O., Taborda, R., Bose, S., Antunes, J., 2005a. Effects of confinement on matrix flow
SC
around a rigid inclusion in viscous simple shear: insights from analogue and numerical modeling. Journal of Structural Geology, 27, 379-396.
NU
Marques, F.O., Taborda, R., Antunes, J., 2005b. 2D rotation of rigid inclusions in confined bulk simple shear flow: a numerical study. Journal of Structural Geology, 27, 2171-2180.
MA
Marques, F.O., Taborda, R., Antunes, J., 2005c. Influence of a low-viscosity layer between rigid inclusion and viscous matrix on inclusion rotation and matrix flow: a numerical study. Tectonophysics 407, 101-115.
D
Means, W.D., Hobbs, B.E., Lister, G.S., Williams, P.F., 1980. Vorticity and non-coaxiality in
TE
progressive deformation. Journal of Structural Geology 2, 371-378. Mercier, S., Molinari, A., 2009. Homogenization of elastic-viscoplastic heterogeneous materials:
1048.
AC CE P
Self-consistent and Mori-Tanaka schemes. International Journal of Plasticity, 25, 1024-
Meng, C., Heltsley, W., Pollard, D.D., 2012. Evaluation of the Eshelby solution for the ellipsoidal inclusion and heterogeneity. Computers & Geosciences, 20, 20-48. Michelitsch, T.M., Gao, H., Levin, V.M., 2003. Dynamic Eshelby tensor and potentials for ellipsoidal inclusions. Proceedings of Royal Society of London A459, 863-890. Molinari, A. 2002. Averaging models for heterogeneous viscoplastic and elastic viscoplastic materials. Transactions of the ASME, Vol 124, 62-70. Molinari, A., Canova, G.R., and Ahzi, S., 1987, A self consistent approach of the large deformation polycrystal viscoplasticity. Acta Mettalurigica, v. 35, 2983-2994. Molinari, A., Ahzi, S., Kouddane, R., 1997. On the self-consistent modeling of elastic-plastic behavior of polycrystals. Mechanics of Materials, 26, 43-62. Molinari, A., Mercier, S. 2004. Homogenization of viscoplastic materials. In: Ahzi, S., et al. (eds.), Multiscale Modeling and Characterization of Elastic-Inelastic Behavior of Engineering Materials, 113-124. 66
ACCEPTED MANUSCRIPT Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 21, 571–574. Mulchrone, K.F., 2007. An analytical solution in 2D for the motion of rigid elliptical particles
PT
with a slipping interface under a general deformation. Journal of Structural Geology, 29(6), 950-960.
RI
Mulchrone, K.F., Walsh, K., 2006. The motion of a non-rigid ellipse in a general 2D deformation.
SC
Journal of Structural Geology 28, 392-407.
Mura, T., 1987. Micromechanics of defects in solids. Martinus Nijhoff Publishers,
NU
Dordrecht/Boston/Lancaster, 587 pp.
Mura, T., 1988. Inclusion problems. Applied Mechanics Reviews 41, 15–20.
MA
Mura, T., Shodja, H.M., Hirose, Y., 1996. Inclusion problems. Applied Mechanics Reviews 49, S118–S127.
Muskhelishvili, N.I. 1953. Some basic problems of the mathematic theory of elasticity. P.
D
Noordhoff Ltd., Groningen.
TE
Nemat-Nasser, S. and Hori, M., 1999. Micromechanics: overall properties of heterogeneous materials. Second Revised Edition, Elsevier.
AC CE P
Passchier, C.W., 1987. Stable positions of rigid objects in non-coaxial flow – a study in vorticity analysis. Journal of Structural Geology 9, 679-690. Passchier, C.W., 1998. Monoclinic model shear zones. Journal of Structural Geology, 20(8): 1121-1137.
Passchier, C.W., Simpson, C., 1986. Porphyroclast systems as kinematic indicators. Journal of Structural Geology 8, 831-843. Passchier, C.W., and Trouw, R.A.J., 1996, Microtectonics. Springer. Pelletier, D., Fortin, A., Camarero, R., 1989. Are FEM solutiobns of incompressible flows really incompressible – (or how simple flows can cause headaches). International Journal for Numerical Methods in Fluids, 9, 99-112. Pennacchioni, G., Di Toro, G., Mancktelow, N.S., 2001. Strain-insensitive preferred orientation of porphyroclasts in Mont Mary mylonites. Journal of Structural Geology 23, 1281-1298. Piazolo, S., Bons, P.D., Passchier, C.W., 2002. The influence of matrix rheology and vorticity on fabric development of populations of rigid objects during plane strain deformation. Tectonophysics 351, 315-329. 67
ACCEPTED MANUSCRIPT Piazolo, S., Passchier, C.W., 2002. Experimental modeling of viscous inclusions in a circular high-strain shear rig: implications for the interpretation of shape fabrics and deformed enclaves. Journal of Geophysical Research 107, 2242.
PT
Qu, J. 1993. The effect of slightly weakened interfaces on the overall elastic properties of composite materials. Mechanics of Materials, 14, 269–281
RI
Qu, J., and Cherkaoui, M., 2006. Fundamentals of micromechanics of solids. John Wiley & Sons,
SC
Inc.
Qu, M., Lu, X., Jiang, D. in review. An optimal method for calculating Eshelby tensors and a
NU
MATLAB package for simulation the motion of viscous inclusions. Submitted to Computers & Geosciences.
MA
Ramsay, J.G. 1967. Folding and fracturing of rocks. McGraw-Hill, New York. Ramsay, J.G., Graham, R.H. 1970. Strain variation in shear belts. Canadian Journal of Earth Sciences, 7, 786-813.
D
Rand, O., Rovenski, V., 2005. Analytical methods in anisotropic elasticity with symbolic
TE
computational tools. Birkäuser, Boston. Rudnicki, J.W., 1977. The inception of faulting in a rock mass with weakened zone. Journal of
AC CE P
Geophysical Research, 82, 844-854. Schmid, D.W., Podladchikov, Y.Y., 2003. Analytical solutions for deformable elliptical inclusions in general shear. Geophysical Journal International. 155, 269-288. Schmid, D.W., Podladchikov, Y.Y., 2004. Are isolated stable rigid clasts in shear zones equivalent to voids? Tectonophysics 384, 233-242. Schmid, D.W., Podladchikov, Y.Y., 2005. Mantled porphyroclast gauges. Journal of Structural Geology 27, 571-585. Spence, D.A., Wilmott, P., 1988. Rotation and deformation of a viscous inclusion in Stokes-flow. Proceedings of the Royal Society of London, Series A 418, 383-403. Spencer, A.J.M., 1980. Continuum mechanics. Longman, London. Stroh, A.N., 1958. Dislocations and cracks in anisotropic elasticity. Philosophic Magazine, 3, 625-646. Stroh, A.N., 1962. Steady state problems in anisotropic elasticity. Journal of Mathematical Physics. 41, 77-103. Synge, J.L. 1957. The hypercircle in mathematical physics. Cambridge University Press. 68
ACCEPTED MANUSCRIPT ten Brink, C.E., Passchier, C.W., 1995. Modelling of mantled porphyroclasts using nonNewtonian rock analogue materials. Journal of Structural Geology 17, 131-146. ten Grotenhuis, S.M., Passchier, C.W., Bons, P.D., 2002. The influence of strain localisation on
PT
the rotation behaviour of rigid objects in experimental shear zones. Journal of Structural Geology 24, 485-499.
RI
Ting, .C.T., 1996. Anisotropic elasticity: theory and applications. Oxford University Press.
SC
Ting, T. C. T. & Lee, V.-G. 1997. The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids. Quaterly Journal of Mechanics and Applied
NU
Mathematics. 50, 407–426.
Treagus, S.H., Treagus, J.E., 2002. Studies of strain and rheology of conglomerates. Journal of
MA
Structural Geology 24, 1541-1567.
Tullis, J. 2002. Deformation of granitic rocks: experimental studies and natural examples. Review of Mineralogy and Geochemistry. 51, 51-95.
D
Walke, K.P. 1993. Fourier integral representation of the Green function for an anisotropic elastic
TE
half-space. Proccdings of the Royal Society of London A. 443, 367-389. Walpole, L.J. 1978. A coated inclusion in an elastic medium. Math Proc Camb, 83:495–506
AC CE P
Weinberger, C., Cai, W., Barnett, D. 2005. Lecture notes – elasticity of microscopic structures. Stanford University (http://micro.stanford.edu/~caiwei/me340b/course_materials.html). White, S.H., Burrows, S.E., Carreras, J., Shaw, N.D., Humphreys, F.J., 1980. On mylonites in ductile shear zones. Journal of Structural Geology, 2, 175-187. Willis, D.G., 1977. A kinematic model of preferred orientation. Geological Society of America Bulletin 88, 883-894.
Willis, J.R., 1964. Anisotropic elastic inclusion problems. Quaterly Journal of Mechanics and Applied Mathematics 17, 157-174. doi:10.1093/qjmam/17.2.157. Xiang, B. and Jiang D., 2013. Small-scale ductile shear zones as transposed rheologically weak domains: A numerical modeling investigation and practical application. Journal of Structural Geology 54: 184-198. Xu, B.-X., Muller, R., Wang, M.-Z. 2011. The Eshelby property of sliding inclusions. Archive of Applied Mechanics., 81, 19-35. Yang, H.C., Chou, Y.T. 1976. Generalized plane problems if elastic inclusions in anisotropic solids. Journal of Applied Mechanics, 43, 424-430. 69
ACCEPTED MANUSCRIPT Yin, H.M., Lee, P.-H., Liu, Y. J., 2014. Equivalent inclusion method for the Stokes flow of drops moving in a viscous fluid. Journal of Applied Mechanics 81, 071010. doi:10.1115/1.4027312.
PT
Zheng, Q.-S., Zhao, Z.-H., Du, D.-X., 2006. Irreducible structure, symmetry and average of Eshelby’s tensor fields in isotropic elasticity. Journal of the Mechanics and Physics of
RI
Solids, 54, 368–383.
SC
Zhou, K., Hoh, H.J., Wang, X., Keer, L.M., Pang, J.H.L., Song, B., 2013. A review of recent works on inclusions. Mechanics of Materials, 60, 144-158.
NU
Zou, W.N., He, Q.C., Zheng, Q.S., 2012. Inclusions in a finite elastic body. International Journal
AC CE P
TE
D
MA
of Solids and Structures 49, 1627–1636.
70
ACCEPTED MANUSCRIPT Figure Captions Figure 1: Stereographic projection showing how, for any given
x , unit vector z in the plane
x (i.e., the x z 0 plane) can be expressed in terms of two mutually orthogonal unit vectors ( α and β ) on the x z 0 plane by z cos α sin β following Synge (1957) and
PT
normal to
SC
RI
cos sin Barnett (1972). When x is expressed in terms of spherical angles and by x sin sin , cos
NU
sin cos cos then α cos , and β sin cos . Green functions and their derivatives ( Gij , Gij ,l , and 0 sin
H i ) for anisotropic viscous materials are expressed as contour integrals along a unit circle on
MA
the x z 0 plane. See text for more details.
D
E Figure 2: A general anisotropic elliptic or ellipsoidal inclusion with viscous stiffnesses Cijkl
TE
embedded in an infinite matrix with viscous stiffnesses Cijkl subjected to far-field uniform
AC CE P
deformation. The mechanical fields inside the inclusion are uniform and denoted by lowercase symbols. The mechanical fields in the matrix (called exterior fields) vary with space (as functions of the position vector x ). The uniform loading or motion at infinity are denoted by uppercase symbols. See text for more details.
Figure 3: Schematic diagrams showing the geometries of the a1 a2 section of an inclusion that goes through a circular state at time t . The ellipse at time t t is reciprocal to the ellipse at time t t . Therefore, the two semiaxes swap identities prior and after the circular state but do not otherwise change directions. See text for more details.
Figure 4: An inclusion embedded in a matrix with planar anisotropy. The matrix rheology is defined by two distinct viscosities, n and s . m n defines the strength of anisotropy. The s far-field general shear flow (pure shearing and simple shearing ) is defined in the
coordinate system x1 x2 parallel to the anisotropy. The inclusion’s principal axes define another coordinate system ( x '1 x '2 ) at angle relative to the x1 x2 system.
71
ACCEPTED MANUSCRIPT Figure 5: Comparison of the interior fields of an inclusion ( R 5 , r 10 ) embedded in a matrix with m 25 ) calculated numerically from the 2D formulation of this paper with results calculated with Fletcher’s analytic solutions. (a) Strain rates differences (between inclusion and
NU
SC
RI
PT
far-field matrix) 11 and 12 . (b) Vorticity difference w11 . (c) Pressure differences p normalized against far-field stress invariant. The solid color curves are numerical results using the 2D formulation of this paper implemented in MATHCAD and the dashed black curves are calculated from analytical expressions (Eq.57). The far-field flow is plane-strain general shear with pure shearing 0.5 and simple shearing 1 both along the plane of anisotropy (Fig.4). Numerical results and those calculated directly using Eqs.57 (dashed lines). Numerical results agree with analytic results perfectly, suggesting excellent accuracy of the numerical approach.
AC CE P
TE
D
MA
Figure 6: Pressure perturbations in and around a strong isotropic inclusion ( R 3 , r 100 ). The pressure is normalized against the matrix deviatoric stress invariant. The far-field flow is pure shearing with horizontal stretching and vertical shortening. The inclusion long axis is horizontal, parallel to the maximum pure shearing axis. (a): Pressure perturbation in isotropic matrix ( m 1 ). (b)-(e): Matrix is anisotropic with m 10 . The anisotropy plane is horizontal in (b), 30°, 60°, and 90° clockwise with respect to the horizontal anisotropy in (c), (d), and (e) respectively. See text for more details.
Figure 7: Angular velocities (normalized against shearing) of rigid inclusions with a circular shape in the x1 x2 plane as functions of the matrix m for different inclusion a3 relative lengths. The far-field plane-straining flow is with pure shearing 0.5 and simple shearing 1 both along the plane of anisotropy (Fig.3). A circular cylinder ( 1:1: ) rotates at a constant angular speed of half the shearing, independent of m . For finite a3 inclusions, the angular velocities increase with m . As the relative length of a3 is reduced, the inclusion’s angular velocity increases. A disk with vanishing a3 has the highest angular velocity (the 1:1:0.001 case).
Figure 8: (a) Angular velocity variations of a rigid triaxial ellipsoid ( 6 : 1: 5 , with the a3 axis parallel to the x3 axis) as functions of its sectional orientation in the x1 x2 plane for different m . The matrix far-field flow is the same as for Fig.6. The horizontal axis is the angle between the inclusion’s sectional long axis and the anisotropy ( x1 axis in Fig.3). A greater m corresponds to a higher angular velocity at any given orientation, except near the stable orientation of the object. (b): Angular velocity variations of ellipsoids with the same sectional aspect ratio of 6:1 but 72
ACCEPTED MANUSCRIPT different relative a3 lengths in a matrix with fixed m =10. A relatively shorter a3 generally
PT
corresponds to a higher angular velocity.
Figure 9: Angular velocities of rigid ellipsoids with the same sectional aspect ratio (2:1) on the x1 x2 plane but varying a3 lengths, in a matrix with m =10 undergoing simple shearing parallel to
SC
RI
the anisotropy plane. The angular velocity of an inclusion depends on its sectional orientation ( ) as well as its relative a3 length. An elliptic cylinder rotates like in an isotropic matrix. As long as a3 is finite, its angular velocity increases with decreasing a3 . A disk with a vanishing a3
NU
rotates at a faster and finite rate.
AC CE P
TE
D
MA
Figure 10: Schematic diagrams showing that in the case of an inclusion with a slippery interface contact with the matrix the shear spin effect must be considered in the angular velocity equation. The system can be regarded as a matrix with a void and an inclusion in the void. (a) A state of the system of a matrix containing a slippery inclusion. A material point (p) on the major radius at the inclusion side of the interface coincides with the material point (p’) on the matrix side of the interface. The points p and p’ are at the same location but are drawn slightly apart across the interface for easy visualization. (b) The state of the system after an increment of deformation. Because of free slip along the interface, p’ has moved along the hoop with respect to p. That is the void’s principal axes have angular velocity swiping through material.
73
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 1
74
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 2
75
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 3
76
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 4
77
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 5
78
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 6
79
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 7
80
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 8
81
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 9
82
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 10
83
S0040-1951(16)30458-9 doi: 10.1016/j.tecto.2016.10.012 TECTO 127279
To appear in:
Tectonophysics
Received date: Revised date: Accepted date:
16 March 2016 7 October 2016 13 October 2016
Please cite this article as: Jiang, Dazhi, Viscous inclusions in anisotropic materials: Theoretical development and perspective applications, Tectonophysics (2016), doi: 10.1016/j.tecto.2016.10.012
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Viscous inclusions in anisotropic materials: theoretical development and perspective applications Dazhi Jiang1,2 Department of Earth Sciences, Western University, London, Ontario, N6A 5B7, Canada
PT
1
2
RI
([email protected])
State Key Laboratory of Continental Dynamics, Department of Geology, Northwest University,
SC
Xi’an, 710069, China
AC CE P
TE
D
MA
NU
Abstract: Theories and numerical solutions for a viscous ellipsoid in an infinite anisotropic viscous medium subjected to far-field homogeneous deformation lie at the heart of selfconsistent homogenization models and multiscale simulations of texture and fabric development in Earth’s lithosphere. There is considerable literature on ellipsoid inclusions, focused on anisotropic elastic materials, published in multi-disciplinary fields. To make this body of work more accessible as well as to advance viscous inclusion studies, an effort is made here to summarize recent advances and to further develop formally more explicit and, where possible, analytic solutions for incompressible viscous materials. The point-force concept and equivalent inclusion method of Eshelby are used together with the powerful Green function approach. This leads to generalized equations for ellipsoid inclusion behaviors in anisotropic materials. In the particular case of isotropic materials, the new mathematical development here enables the use of existing methods for elastic materials to get solutions for corresponding viscous inclusion problems efficiently and accurately. A 2D formulation is also presented for elliptic cylinders in plane-straining flows of anisotropic materials, using the same Green function method that is adopted for 3D inclusions. The 2D formulation is benchmarked with existing analytic solutions. A reconnaissance investigation to compare the behaviors of 2D elliptic inclusions and triaxial ellipsoids in a matrix of planar anisotropy undergoing far-field plane-straining flows suggests that conclusions drawn from 2D cannot be applied to 3D in anisotropic cases. The application of the viscous inclusion theory to the rheologically heterogeneous and non-linear lithosphere is discussed. By regarding a heterogeneous element as an ellipsoidal inclusion embedded in a hypothetical homogeneous equivalent matrix whose effective rheology is obtained through micromechanical homogenization and assuming that the conditions of scale separation are satisfied, the viscous inclusion theory forms the backbone of multiscale models for natural deformation and the accompanying texture and fabric development. Keywords: fabrics; viscous materials; multiscale modeling; inclusion problem; Eshelby solution; anisotropic viscosity; homogenization; rheology; flow partitioning
1
ACCEPTED MANUSCRIPT 1. Introduction
AC CE P
TE
D
MA
NU
SC
RI
PT
The instantaneous mechanical state and progressive evolution of a solitary inclusion in an infinite uniform material subjected to far-field homogeneous deformation, referred to as the “inclusion problems” below, have been studied extensively for a century and remain an active research subject today (see reviews by Mura, 1987, 1988; Ting, 1996; Zhou et al., 2013). Inclusion problems are significant in a wide variety of fields both theoretically and in practical applications. Micromechanics (e.g., Mura, 1987; Nemat-Nasser and Hori, 1993), the most recent development in applying continuum mechanics theories to real materials (Qu and Cherkaoui, 2006, p.2), stems from the seminal work of Eshelby (1957) on elastic inclusions. Earth scientists have applied and contributed to the advance of inclusion theories in the study of fabric development in the crust and strain and mixing of lithospheric heterogeneities in the mantle for half a century (e.g., Gay, 1968; Ramsay 1967, p.225; Hobbs et al., 1976, p.285; Passchier and Simpson, 1986; Passchier 1987; Spence et al., 1988; Spence and Wilmott, 1988; ten Brink and Passchier, 1995; Piazolo et al., 2002; Piazolo and Passchier, 2002; ten Grotenhuis et al., 2002; Treagus and Treagus, 2002; Ceriani et al., 2003; Mandal et al., 2003, 2005; Marques and Bose, 2004; Johnson et al., 2009; Li and Jiang, 2011; Mancktelow 2011, 2013; Xiang and Jiang, 2013; Chen et al., 2014). Among the considerable literature of inclusion studies, perhaps most influential to earth scientists are the work of Jeffery (1922) on the motion of a rigid ellipsoid in incompressible isotropic Newtonian fluids, its generalizations (Bretherton, 1962; Willis, 1977; Ferguson, 1979), and Eshelby’s (1957) work on isotropic elastic inclusions. The influence of the latter is largely through the works of Bilby et al. (1975) and Bilby and Kolbuszewski (1977) who applied Eshelby’s (1957) theory to isotropic incompressible viscous fluids. In the past two decades, Eshelby’s inclusion formalism has been extended to anisotropic viscous inclusions in nonlinear materials (Molinari et al., 1987; Lebensohn and Tomé, 1993). This has formed the backbone for the self-consistent viscoplastic model for texture development in poly-crystalline materials (Lebensohn and Tomé, 1993), the self-consistent MultiOrder Power Law Approach for multiscale fabric development in the heterogeneous lithosphere (Jiang and Bentley, 2012; Jiang 2014), and many self-consistent models of homogenization to obtain the bulk rheological properties from constituent properties (e.g., Molinari, 2002). Most works on inclusion problems in the past 5 decades are focused on elliptic inclusions in plane-strain (or plane-stress) and generalized plane deformation of anisotropic elastic materials (Jaswon and Bhargava 1961; Willis, 1964; Bhargava and Radhakrishna, 1964; Chen 1967; Yang and Chou, 1976; Hwa and Ting, 1989; Ting, 1996). These works adopt the complex variable formalism for plane problems of elasticity (Lekhnitskii 1950; Mushkelishvili 1953; Stroh 1958, 1962; Green and Zerna 1968). As will be shown in this paper, the work of Jaswon and Bhargave (1961) on elliptic inclusions in plane-strain (or plane-stress) deformation of isotropic elastic materials is directly applicable to corresponding 2D inclusions in isotropic viscous materials. However, it remains open how the works on 2D anisotropic elastic inclusions may be applied to corresponding viscous materials. The difficulty lies in the requirement of 2
ACCEPTED MANUSCRIPT
PT
considering incompressibility for viscous materials, which breaks the equivalence between linear elastic and incompressible Newtonian viscous problems. Molinari et al. (1987) first formally solve the incompressible 3D viscous inclusion problems in anisotropic viscous materials. Lebensohn et al. (1998) further develop this approach. Fletcher (2004, 2009) obtains analytic solutions for an elliptic isotropic viscous inclusion embedded in a viscous matrix with planar anisotropy, addressing incompressibility explicitly.
NU
SC
RI
Major progresses have also been made in developing numerical algorithms and programs to solve the equations of Jeffery (1922), Eshelby (1957, 1959), and more general ones for anisotropic materials (e.g., Freeman 1985, 1987; Ježek 1994; Jiang 2007a, 2007b, 2012; Meng, et al., 2012; Qu et al., 2016). Complete algorithms and software packages are now available for simulating single and multi-inclusion systems in general anisotropic viscous materials including non-linear power-law materials (Jiang 2012, 2013, 2014; Qu et al., 2016).
AC CE P
TE
D
MA
It is a bit challenging for one to follow the development on inclusion studies as the papers are published across a broad range of disciplinary fields and formulated in heavy and variable mathematical formalisms from one paper to another. Here I attempt to present a treatise on the subject. I review recent advances, develop more explicit formal and, where possible, analytical solutions for viscous inclusions in anisotropic materials in the narrative, and discuss the applications of the theoretical developments. The Green function method is used throughout the paper to render a more consistent mathematic notation. Many new equations governing the behavior of inclusions in anisotropic materials are established. New mathematical relations discovered in this paper enable efficient and accurate numerical solutions to viscous inclusion problems using existing equations for elastic inclusions. It is demonstrated that all existing results for ellipsoidal and elliptic inclusions can be derived readily from the general formulation as subsets or limiting cases. The more explicit formulation here greatly simplifies the efforts of numerical implementation of the inclusion theory in mathematic applications such as MATLAB, MATHCAD, and MAPLE.
2. Non-equivalence between compressible linear elasticity and incompressible linear viscosity It has been long held that linear elastic theories and linear viscous theories are formally equivalent (e.g., Spencer 1980, p.110-118). This implies that the former can be used, en masse, for the latter. However, elastic theories generally consider compressible solids whereas viscous materials are commonly assumed to be incompressible. This asymmetry does not cause undue difficulties for isotropic materials because their elastic stiffnesses are reduced to two scalar quantities (e.g., the bulk and shear moduli) which occur explicitly in elastic equations. One can address incompressibility easily by considering the limit case as the bulk modulus tends to infinity (or, equivalently, the Poisson’s ratio to 0.5). This was done by Bilby et al. (1975) and 3
ACCEPTED MANUSCRIPT
SC
RI
PT
Bilby and Kolbuszewski (1977). Later in this paper, it will be shown that Jaswon and Bhargava’s (1961) work on an elliptic inclusion in a 2D isotropic elastic medium can also be extended to viscous materials in a similar way. However, general anisotropic elastic materials have many non-vanishing stiffness terms and corresponding elastic equations are also more complex, generally expressed in convoluted and implicit forms. In 2D theories on anisotropic elastic solids, equations are also often expressed in heavy complex variable formalisms. This makes it non-trivial to consider incompressibility as a limiting case. It becomes more straightforward to regard incompressibility as an additional kinematic condition and solve the differential equation set for inclusion problems anew. The fundamental equations for linear anisotropic elasticity are listed below:
1 ui, j u j ,i 2
NU
ij
(1a)
el el ij Cijkl kl , ij M ijkl kl
constitutive law (generalized Hooke’s law)
(1b)
el Cijkl uk ,lj fi 0
equilibrium equations
(1c) (1d)
TE
el el el el s Cijpq M pqkl M ijpq C pqkl J ijkl
D
MA
compatible strain field
where ui , ij , ij , f i are respectively the displacement, elastic strain tensor, Cauchy stress
AC CE P
el el tensor, and body force, all referred to a fixed coordinate system xi i 1, 2, 3 ; Cijkl and M ijkl
are, respectively, 4th-order elastic stiffness and compliance tensors which are inverse of each s other in the sense of (1d), and J ijkl is the 4th-order symmetric identity tensor defined as
1 ik jl jkil ( ij being the Kronecker delta). The notation used in (1) and throughout 2 this paper follows Jiang (2014). A combination of Cartesian tensor notation and a standard tensor/matrix notation are used. Repeated indices in (1) and throughout the paper imply summation unless declared otherwise. A comma stands for differentiation with respect to the u Eulerian coordinate (e.g., ui , j i ). x j s J ijkl
To many readers, Eqs.1 are familiar. Eq.1a ensures that the elastic strain field is compatible. The generalized Hooke’s law (1b) relates the Cauchy stress tensor (symmetric, having 6 independent components) and the elastic strain tensor (also symmetric and having 6 el el independent components). Both Cijkl and M ijkl can be converted into 6 6 non-singular matrices (not the Voigt matrix notation) using a set of orthonormal second-order tensors so that their inversion (1d) calculations can be performed readily (see Jiang 2013 for more details).
4
ACCEPTED MANUSCRIPT In comparison, the corresponding fundamental equations for linear anisotropic viscosity are:
1 ui, j u j ,i 2
kk uk ,k 0
(2a)
PT
compatible strain rate field
RI
incompressibility
v v ijd Cijkl kl , ij M ijkl kld
rheology equations
(2c) (2d)
Cauchy stress field
v Cijkl uk ,lj p, j ij fi 0
equilibrium equations
NU
ij pij ijd
(2e) (2f)
MA
v v v v d Cijpq M pqkl M ijpq C pqkl J ijkl
(2b)
SC
ij
where, for easy comparison with (1), ui , ij , ij now stand for, respectively, the velocity,
D
viscous strain rate tensor, and the Cauchy stress tensor; ijd in Eq.2c stands for the deviatoric
AC CE P
TE
stress tensor; p is pressure. Eq.2b states the incompressibility condition. The rheology equations (2c) relate the deviatoric stresses (having 5 independent components) to the viscous strain rates v v (also having 5 independent components). The viscous stiffness ( Cijkl ) and compliance ( M ijkl ) are
1 d s J ijkl ij kl . inverse of each other in the deviatoric subspace (2f) where J ijkl 3 Note the viscous equilibrium equation (2e) differs from the elastic equilibrium equation (1c). Solving (2e) with the condition of (2b) for given boundary or initial conditions (a viscous problem) is different from solving (1c) for given boundary conditions (an elastic problem). This is why an elastic solution, in general, cannot be used directly for a corresponding viscous problem. In principle, one may still regard viscous incompressibility as a limiting case using a numerical approximation called ‘the penalty approach’ (e.g., Hutchinson 1976; Bilby et al. 1975; Lebensohn and Tomé 1993) which consists in using a modus operando ‘bulk viscosity’ ( b in the following expressions) to make an incompressible viscous material slightly compressible. Incompressible condition is approximated then by making b very large compared to other viscosity components. However, this method has serious drawbacks as explained below. A general penalty approach formulation is presented as follows. One can add a ‘slight compressibility term’ to express the viscous stiffness and compliance tensors as:
5
ACCEPTED MANUSCRIPT p p m v , M ijkl Cijkl 3b J ijkl Cijkl
1 m v J ijkl M ijkl 3b
PT
1 m v ij kl , and the superscript ‘p’ stands for the penalty approach. Unlike Cijkl where J ijkl and 3 v p p , the two tensors Cijkl and M ijkl are invertible in the Cauchy space (i.e., M ijkl
RI
p p p p s p p ). With Cijkl and M ijkl , equations for a viscous problem become Cijpq M pqkl M ijpq C pqkl J ijkl
SC
identical to those for a corresponding elastic problem and elastic solutions can be used. It would appear simple to get solutions for incompressible viscous problems by letting b tend to infinity.
NU
However, the difficulty lies in that fact that, except for isotropic materials, it is not possible to simply let b to get exact limit form expressions because of the implicit and complex nature of the solutions. This means that one has to use a large b numerically to obtain
MA
approximate solutions. Clearly a larger b would result in a better approximation of incompressibility. Unfortunately, because of the finite precision of the computer, too large a b v relative to the norm of Cijkl will deteriorate the numerical calculations due to ill conditioning of
TE
D
p and the inevitable accumulation of round-off errors (e.g., Pelletier et al., 1989). Cijkl
AC CE P
Because of this difficulty with the penalty approach, I shall consider incompressibility explicitly for general anisotropic viscous materials, advancing the work of Molinari et al. (1987) and Lebensohn et al (1998). This requires solving (2e) with condition (2b) afresh, which is the effort in the following section, through the powerful Green function method.
3. General solution for a single viscous inclusion in an infinite viscous matrix Two keystone ingredients in Eshelby’s (1957, 1959, 1961) [for a collection of Eshelby’s works, the reader is referred to Markenscoff and Gupta (2006)] treatment of the elastic inclusion problems are 1) the point force concept that allows an inclusion with eigenstrain to be treated as a distribution of body force and 2) the equivalent inclusion method that regards an inhomogeneous inclusion as a homogeneous one with proper eigenstrain. For a more detailed explanation of these two ingredients, the reader may refer to Jiang (2013, 2014). These ingredients are applied to solve the general three-dimensional viscous inclusion problems here, considering incompressibility explicitly. In the following, “viscous materials” shall always mean general “incompressible anisotropic viscous materials”. 3.1 Green functions for viscous materials The point force concept leads Eshelby (1957, 1959) to the use of elastic Green function for the solution of a general elastic inclusion problem. In linear elasticity, a single Green function 6
ACCEPTED MANUSCRIPT for displacement is used (e.g., Eshelby 1957). To address incompressibility explicitly, Molinari et al. (1987) and Lebensohn et al. (1998) use two separate Green functions: Gij x, x ' for velocity, and H i x, x ' for pressure. They relate, respectively, the velocity ui x and pressure
PT
p x at location x to the action of a unit point force at the source location x ' . With the Green
functions known, the velocity and pressure fields are then expressed by the convolution relations
RI
ui x Gij x, x ' f j x ' dV x ' and p x H i x, x ' fi x ' dV x ' , where f x ' is the body
kl* x ' x ' j
(e.g., Walker, 1993; Li and Wang, 2008, p.89;
NU
body force distribution fi x ' Cijkl
SC
force distribution in the deforming body and the volume containing all body force sources. The presence of an eigenstrain field kl* x ' in the medium is replaced by a layer of fictitious
where Cijkl is the elastic or viscous stiffness tensor of the material). With the equivalent
D
MA
inclusion method, a heterogeneous inclusion is then regarded as an equivalent homogeneous one with the appropriate eigenstrain (Eshelby 1957, see Jiang 2013, 2014 for more explanation). Therefore, the first important step toward solving the viscous inclusion problem is to find out expressions for the Green functions.
TE
Because in an infinite medium Gij x, x ' Gij x x ' and Hi x, x ' Hi x x ' , the Green functions are written as Gij x and H i x for simplicity below, setting the source point at the
AC CE P
origin. The Green functions for an infinite viscous material are solutions to the following auxiliary set of linear partial differential equations (Molinari et al. 1987, Eq.27 there): CijklGkm, jl x H m,i x im x 0
(3a)
Gkm,k x 0
(3b)
with boundary conditions
Gkm 0 , H m 0 at infinity
(3c)
where Cijkl is the viscous stiffness tensor of the material and x is the Dirac delta function. Applying Fourier transform to (3) converts them to the following algebraic system:
Aik K 2 gkm zi iKhm im
(4a)
zk K 2 gkm 0
(4b)
where 7
ACCEPTED MANUSCRIPT Aik Cijkl z j zl
(5a)
gkm (K )
G x exp iK x dx
(5b)
km
PT
hm (K )
H x exp iK x dx m
RI
(5c)
SC
The symmetric Aik is called the Christoffel stiffness tensor (e.g., Barnett, 1972; Walker, 1993),
gkm (K ) and hm (K ) are the Green functions in Fourier space, K Fourier wave vector, z unit
K ), and i 1 . The Green functions in Fourier space, gij and K hi respectively, can be solved from (4): i Aij ( z) , hi (K ) i 2 K K
MA
gij (K )
NU
vector parallel to K (i.e., z
(6)
ζ A T z
z 0
1
(7)
AC CE P
A T ζ
TE
D
where A , ζ , and , all functions of z , are submatrices of the following matrix:
We need Green functions in real space which are inverse Fourier transforms of (6). Analytical expressions for Gij x and H i x are known for isotropic materials (e.g., Bilby et al., 1975; Mura, 1987, p.22). Incompressible Green functions and their derivatives are unknown for anisotropic materials. In Appendix A, expressions of Gij , Gij ,l , and H i for viscous materials are obtained for the first time (A17) following the method of derivation of Barnett (1972) and Mura (1987) for elastic materials. The results are summarized below: Gij x
1
4 2 r 0
Aij d
(8a)
1 H i x 2 2 Fd i 4 r 0
(8b)
1 Gij ,l x 2 2 xl Aij zl M ij d 4 r 0
(8c)
8
ACCEPTED MANUSCRIPT where
r x , xi xi , Fi Ain mn m Aim i m xm , r
PT
M ij Aim A jn mn Aim j A jm i x m
RI
with C p q z p x q zq x p .
SC
The quantities Aij , Fi , and M ij in (8) are all functions of z which in the x z 0 plane
NU
can be expressed in terms of two mutually orthogonal unit vectors α and β (Synge, 1957; Barnett, 1972) shown in stereographic projection in Fig.1. The integrations of (8) are line integrals along a unit circle in the x z 0 plane.
Aij
1
ij
zi z j , Fi x i , and M ij
MA
For an isotropic viscous material with viscosity , it can be shown (Appendix A) that
1
z x i
j
z j x i . Explicit integration of (8) gives the
1 8 r
ij
xi x j
AC CE P
Gijiso x
TE
D
following well-known results (e.g., Bilby et al., 1975; Buryachenco, 2007, p.53, 56; Yin et al., 2014):
Gijiso,l x
H iiso x
1
8 r 2
x ij
l
jl
(9a)
x i il x j 3x i x j x l
xi 4 r 2
(9b)
(9c)
where the superscript “iso” stands for isotropic viscous materials. More explicit forms than (8) for Green functions and their derivatives are still not available for viscous materials. Some close-form expressions for elastic Gij and its derivatives have been found for transverse isotropic, orthotropic, and cubic materials, and on the symmetry plane of monoclinic materials (e.g., Ting and Lee, 1997 and references therein; Buryachenco, 2007). Unfortunately, because of the lack of full equivalence between compressible linear elasticity and incompressible linear viscous problems highlighted in Section 2, elastic Green function and their derivatives cannot be used readily for incompressible viscous materials. Nevertheless, the numerical evaluation of (8) is straightforward using simple quadrature techniques in standard mathematics applications such as MATHCAD or MATLAB. A MathCad worksheet to evaluate (8) for general anisotropic viscous materials is provided as online 9
ACCEPTED MANUSCRIPT Supplement Materials with this paper. The challenge remains in self-consistent simulations (below) where the integrations must be evaluated for many elements at every increment of deformation with a nontrivial amount of computation (8c, particularly) at each evaluation.
PT
3.2: Eshelby tensors and auxiliary quantities
RI
To present the solutions for the general viscous inclusion problems, two auxiliary tensors are defined: the second-order ij , related to the calculation of pressure (Lebensohn et al. 1998,
SC
their Eq.B14; Molinari et al., 1997, their Eq.31), and the fourth-order Tijkl (Molinari et al., 1987), sometimes referred to as Green’s interaction tensor (e.g., Lebensohn et al. 1998). They are respectively:
H x x 'n x ' dS x ' i
j
MA
NU
ij x H i , j x x ' dV x '
Tijkl x Gik , jl x x ' dV x '
G x x 'n x ' dS x ' ik ,l
j
(10a)
(10b)
TE
D
where the integrations are over the ellipsoid volume (first equal sign) or over the ellipsoid surface (the second equal sign) S x ' with n j being plane-normal unit vector component of the
AC CE P
area element dS on the ellipsoid surface. The second equal signs in (10a, b) are due to application of Gauss’s theorem (e.g., Korn and Korn, 1968, p.5.6.2). Two 4th order tensors, S and Π , referred to, respectively, as the symmetric and antisymmetric Eshelby tensors for viscous materials are defined as follows (Jiang 2013):
S J d : T : C and Π J a : T : C
(11)
d where C is the viscous stiffness (viscosity) tensor of the matrix, J ijkl is the same as in Eq.2f,
a J ijkl
1 ik jl jkil is the fourth-order anti-symmetric identity tensor as used in Jiang (2013, 2
2014). The operation “:” stands for double-index contraction between tensors. Eqs.11 are also equivalent to the following expressions used in Appendix B:
S Ts : C and Π Ta : C
10
ACCEPTED MANUSCRIPT s where Tijkl
1 1 a Tijkl T jikl Tijlk T jilk and Tijkl Tijkl T jikl Tijlk T jilk (Lebensohn et al., 4 4
PT
1998).
NU
SC
RI
An amazing result of Eshelby (1957) is that the mechanical fields inside the ellipsoid are uniform. This remarkable property, called the Eshelby uniformity (Zheng, et al., 2006), has played a key role in micromechanics. Note the Eshelby uniformity in ellipsoids (or ellipses in 2D) is a proven fact, not a conjecture, since Eshelby (1957). Eshelby (1961) did further conjecture ‘…Among closed surfaces, the ellipsoid alone has this convenient property…”. This conjecture of Eshelby’s has been proved for isotropic materials by Kang and Milton (2008) and Liu (2008). Eshelby’s conjecture remains open for anisotropic materials. As a consequence of the Eshelby uniformity, both ij and Tijkl are constant tensors for
MA
points within the ellipsoid (called interior points below). Accordingly, Sijkl and ijkl are also constant tensors for interior points. The mechanical fields in the matrix (exterior points) vary with space. To differentiate, the corresponding tensors for exterior points are denoted by ijE (x) ,
TE
D
E E E E , and ijkl (where superscript “E” stands for exterior field). Sijkl (x) , ijkl (x) , or simply ijE , Sijkl
AC CE P
3.3: Formal solutions for the interior and exterior fields of an incompressible linear viscous ellipsoid Finally, the mechanical fields inside the ellipsoid are related to the remote fields by the following set of equations: ε J d S1 : C1 : σ
(12a)
w Π : S1 : ε
(12b)
p Λ : C : S1 : ε
(12c)
1
where the mechanical field quantities inside the ellipsoid are represented by lowercase symbols ε , σ , w , and p , standing for, respectively, the strain rate tensor, the deviatoric stress tensor, the vorticity tensor, and the hydrostatic pressure; their uppercase counterparts Ε , Σ , W, and P shall stand for corresponding field quantities in the far-field matrix (Fig.2). The tilde over a quantity stands for the difference of the quantity between the field inside the ellipsoid and that in the far-field matrix, e.g., ε ε Ε . Eqs.12 collectively are called “solutions to an interior point Eshelby problem” or simply “the interior solutions”. 11
ACCEPTED MANUSCRIPT The exterior point fields are also expressed by a set of interaction equations as follows: (13a)
w(x) ΠE (x) : S1 : ε
(13b)
PT
ε(x) SE (x) : S1 : ε
(13c)
RI
p(x) Λ E (x) : C : S1 : ε
SC
This set is often referred to as “solutions to an exterior point Eshelby problem” or simply “the
NU
exterior solutions”.
MA
4. Viscous tensor identities
D
In this section we establish some previously unknown identities among viscous tensors, mostly for isotropic materials. These identities make the numerical evaluation of viscous inclusion problems more efficient and accurate.
TE
4.1: Interior field quantities
For the interior fields, ij and Tijkl , both constant tensors, can be expressed respectively
AC CE P
as (Lebensohn et al., 1998; Weinberger et al., 2005):
ij
a1a2a3 i z j 3 sin dd 2 0 0
(14a)
Tijkl
a1a2a3 Aik z j zl 3 sin dd 2 0 0
(14b)
where ai i 1, 2, 3 are the three semi-axes of the ellipsoid,
a1z1
2
a2 z2 a3 z3 , 2
cos sin and unit vector z is expressed in terms of spherical angles by z sin sin . These two cos integrals can be numerically evaluated efficiently using a product Gaussian quadrature (Jiang, 2013, 2014; Atkinson, 1982) or a Lebedev quadrature depending on the shape of the ellipsoid (see Qu et al., 2016).
12
2
ACCEPTED MANUSCRIPT From the definition of (7) it is easy to show that zi i 1 . Combining this with (14a), we a1a2a3 4
2
2
0 0
0 0
3 sin dd . Applying the identity
Michelitsch et al. (2003), we get the following identity:
kk 1
sin dd
4 due to a1a2a3
(15)
RI
(for general viscous materials)
3
PT
have kk
It is clear from (14a) that ij ji for anisotropic materials. For isotropic materials
SC
however, it can be shown (Appendix A) that i zi and the integrand in (14a) is odd with
ij 0 (if i j )
NU
respect to z if i j . This means: (for isotropic materials)
(16)
TE
D
MA
Since Bilby et al. (1975), it has been a common practice to use Eshelby tensor expressions for isotropic elastic materials to calculate the corresponding viscous Eshelby tensors by setting the Poisson’s ratio term ( in the following equations) to 0.5. In Appendix B, it is shown that the viscous and elastic symmetric Eshelby tensor components for isotropic materials are equal for off-diagonal components only. For diagonal components (i.e., components with indices like Siijj , including i j ), the following relations hold (see Appendix B for detail
ii
AC CE P
derivation):
1 el Sii11 Siiel22 Siiel33 0.5 3
el Siijj ii Siijj
0.5
(no sum)
(17a)
(no sum)
(17b)
where the superscript “el” stands for elastic. The antisymmetric Eshelby tensors for viscous and elastic isotropic materials are indeed equal: ( el ) ijkl ijkl
(18)
0.5
4.2: Exterior field quantities E For the exterior field solutions (13), the corresponding ijE , Tijkl , and Eshelby tensors E E , ijkl can be calculated numerically using integral expressions in (10) together with Green Sijkl
function expressions given in (8). The computation is however intensive even for isotropic materials (Eshelby 1959). Mura (1987) attempted to make the exterior solution for isotropic 13
ACCEPTED MANUSCRIPT
SC
RI
PT
elastic materials more explicit. He obtained rather long expressions (his Eqs.11.41 and 11.42) containing elliptic integrals (Mura 1987, p.85) and their derivatives that would still require considerable computational effort to evaluate. Thereafter, there have been two major advances for the exterior solution for isotropic elastic materials. First, Meng et al (2012) finally closed those unfinished derivatives in Mura’s expressions. They reduced the terms in Mura’s Eqs.11.41 and 11.42 that contain elliptic integrals and their derivatives to two simpler complete elliptic integrals of the first and second kind. This has greatly improved the efficiency and accuracy in numerical calculations of the exterior elastic field. The complete displacement field can be calculated using the far more explicit expression of Meng et al. (2012, Eq.21 there). Second, Ju and Sun (1999) developed a different and much more explicit approach to calculate S E ijkl [their
NU
G ijkl ] for isotropic elastic materials. Their approach allows the exterior strain and stress fields
(but not the complete displacement field) to be calculated even more efficiently. However, Ju and Sun’s (1999) approach cannot yield the complete displacement field because E ijkl is not
MA
calculated.
TE
D
The development below allows the complete viscous field quantities (deviatoric stress, strain rates, vorticity, and pressure) outside the inclusion to be computed efficiently and accurately. 4.2.1: Direct calculation of exterior field quantities
AC CE P
E To solve the exterior fields for viscous inclusions, we must obtain E ijkl , E ij , and Sijkl .
Let us investigate these quantities for incompressible isotropic viscous materials. First, because the Green functions for isotropic materials are available in analytic forms (9), it is straightforward, albeit a bit tedious, to get the following integral expressions: E Sijkl (x)
1 4
(19a)
xi x l jl x i x k ik x j x l il x j x k dV
(19b)
1 3 ij kl ik x j x l jk x i x l il x j x k jl x i x k 2 ij x k x l kl x i x j 15x i x j x k x l dV 3 2
r
E ijkl (x)
ijE (x)
3 8
1 4
1
r
1
r
3
where r x x ' , x i
3
jk
3x x i
j
ij dV
(19c)
xi x 'i , and the integration is over the volume of the ellipsoid x x'
x '12 x '2 2 x '32 2 2 1. a12 a2 a3
14
ACCEPTED MANUSCRIPT Second, using the two identities below:
xi x j x k 2 r
RI
xi x j x k xl 1 1 3 il x j x k jl x i x k kl x i x j 3 r 5 r x 'l
PT
x i x j 1 ij xi r3 3 r 3 x ' j r 2
SC
one can obtain the following relations with the Gauss theorem: xi x j 1 dV ij I 0 ij 3 r 3
(20a)
NU
xi x j x k x l 1 dV il jk jl ik kl ij ijkl 3 r 5
(20b)
MA
xi n dV x i x j x k nl , ij 2 j dS , and ijkl dS 3 r r r2
(20c)
TE
D
I0
AC CE P
Third, substituting (20) into (19) and then with simple variable substitutions to convert integrals over an ellipsoid surface to ones over a unit sphere, we get the following expressions:
1 4
E Sijkl (x)
E ijkl (x)
ijE (x)
1 ij kl 2 ik jl jk il il jk jl ik 3ijkl
1 jk il jlik ik jl il jk 8
(21a)
(21b)
1 ij 4
(21c)
where ij and ijkl are functions of x , which can be expressed explicitly as:
xi ai zi sin 2 2 2 0 x a z x a z x a z 2 2 2 3 3 3 1 11
2
ij (x) a1a2a3 0
2
ijkl (x) a1a2a3
aj
2
xi ai zi x j a j z j xk ak zk sin
d d
zl d d 2 2 2 2 a l x a z x a2 z2 x3 a3 z3 1 1 1 2
0 0
zj 3
5
15
(21d)
(21e)
ACCEPTED MANUSCRIPT and unit vector z is the same as in (14). With (21), solving the exterior fields eventually boils down to the evaluation of ij and
ijkl (21d and e), both of which can be performed efficiently and accurately with built-in
RI
PT
quadrature in Matlab, Maple, and Mathcad or with a combination of product Gaussian quadrature and Lebedev quadrature for extremely elongate or flattened ellipsoids (Qu, et al., 2016).
SC
4.2.2: Efficient calculation of exterior field quantities using Ju and Sun’s (1999) method The solution to exterior viscous inclusion problems in isotropic materials can be further
NU
E simplified. First, it is clear from Eq.19a that for incompressible isotropic viscous materials, Sijkl E E Sklij has the major symmetry, i.e., Sijkl . Whether this is also true for incompressible anisotropic
MA
materials remains open. Second, from (20a) it is clear that ij ji and kk 0 . Therefore ijE E 0 ) for isotropic materials. In summary, for must be symmetric ( ijE Eji ) and traceless ( kk
D
isotropic incompressible materials we have,:
(22)
TE
E E E Sijkl Sklij 0 , ij ji , kk 0 , ijE Eji , and kk
E From Eq.21a, it is easy to confirm that S kl 0 . This, together with the major symmetry
AC CE P
E of Sijkl , means that for incompressible materials, we have: E E E S kl Sij Sii 0 (sum applies to repeated Greek letter only)
(23)
Eqs.22 and 23 are now combined with Ju and Sun’s (1999) approach to greatly simplify the calculation of exterior field quantities. Ju and Sun (1999) develop an explicit approach to calculate the exterior-point symmetric Eshelby tensor for isotropic elastic materials [their G ijkl (x) ]. Comparing their Eq.4 (see also Li and Wang, 2008, p.105) with (19a) above, we get the following relation: E Sijkl (x) G ijkl (x)
0.5
ijE (x) kl
(24)
which can be regarded as a generalized relation of (B9) that we obtain for the interior field quantities in Appendix B.
16
ACCEPTED MANUSCRIPT E Applying the conditions of Eqs.23 and kk 0 to (24), we arrive at the following E E relations which allow all viscous tensors, ijE (x) , Sijkl (x) , and ijkl (x) , required for the solution
1 G ij (x ) , 3 0.5 E ij (x ) 1 G ij (x ) Gij ( x) , 0.5 3 0.5
ijE (x) kl ,
RI
(i j )
(25a)
except i j k l
i
j, i k , and j k
NU
E E E Siiii (x) Siijj (x) Siikk (x) ,
1 jk ilE jl ikE ik Ejl il Ejk 2
(25b) no sum
(25c) (25d)
MA
E ijkl (x)
(i j)
SC
E Sijkl (x) G ijkl (x)
PT
of exterior fields at a point, to be calculated simply from G ijkl (x) of Ju and Sun (1999):
TE
D
This greatly reduces the computation resources and improves the accuracy of the results. A MathCad worksheet to calculate exterior field quantities based on (25) is provided as online Supplement Materials with this paper.
AC CE P
5. Equations for rotational behaviors of ellipsoids in anisotropic viscous materials and principles of numerical implementation In this Section, we use the results of Section 4 to establish some new equations governing the motion of inclusions in anisotropic viscous materials and discuss numerical simulations of inclusion evolution. 5.1 Angular velocity of a viscous ellipsoid The rotation of an ellipsoidal inclusion in a viscous flow is determined by both the vorticity in the ellipsoid (12b), which determines the angular velocity of the material of the ellipsoid, and the migration of the principal axes through material due to ellipsoid’s shape change. The latter is described by a shear spin, denoted by w s , defined as (Goddard and Miller, 1967; Bilby and Kolbuszewski, 1977; Jiang 2012): ai2 a 2j 2 2 ij wijs ai a j w ij
ai a j
(no sum)
ai a j
17
(26)
ACCEPTED MANUSCRIPT The angular velocity tensor for an ellipsoid, written in ellipsoid’s coordinate system, ij , is then (e.g., Kocks, 1998; p. 395): ij wij wijs
PT
(27)
RI
Combining (12b) and (27), the angular velocity equation for the most general case of a deformable anisotropic inclusion in a dissimilar viscous anisotropic matrix is: (28)
SC
1 ij Wij ijmn Smnkl kl wijs
rigid in anisotropic
1 Wij ijmn Smnkl Ekl
MA
ij
NU
For a rigid inclusion, kl 0 (equivalent to kl kl ), (28) reduces to the following equation (29)
which is a generalization of Jeffery’s (1922) equation to an anisotropic viscous matrix. When the
D
1 matrix is isotropic, the term ijmn Smnkl Ekl becomes
TE
reduces to the well-known Jeffery’s result, ij
a I2 a J2 ij (Appendix C) and Eq.29 further a I2 a J2
rigid in isotropic
Wij
aI2 aJ2 ij (Jiang, 2012, Eq.15 aI2 aJ2
AC CE P
there). Here and below the notation adopted by Mura (1987) is used, whereby repeated lowercase indices are summed up from 1 to 3, while uppercase indices which take on the same numbers as the lowercase ones are not summed up. The behavior of an inviscid ellipsoid is another limiting case where the deviatoric stresses in the inclusion vanish (i.e., σ Σ ). Submitting this into (12a) and (12b), we get, after arranging: ε general J d S : Ε , 1
inviscid
ε general S1 J d : Ε 1
(30a)
inviscid
w general Π : J d S : Ε 1
(30b)
inviscid
ij
general inviscid
1 a2 a2 d d Wij ijkl I2 J2 J ijkl J S klmn Emn aI aJ
s
(30c)
5.2 Physical meaning of w , w , and justification for setting wijs wij when ai a j 18
ACCEPTED MANUSCRIPT
RI
PT
The vorticity in the inclusion w , which can be calculated from (12b), defines the rotation of the material element of the ellipsoid. w is always well defined whether the inclusion is a triaxial ellipsoid, spheroid, or sphere. Because the flow inside the ellipsoid is generally nonsteady, w can be decomposed to an internal vorticity and a spin in the sense of Astarita (1979) and Means et al (1980). The internal vorticity is the ‘objective’ (reference-frame invariant) part of vorticity that should be used to define the characteristics of flows (Jiang, 1999, 2010; Hobbs and Ord, 2015, p.77). Xiang and Jiang (2013) considered the decomposition of w in their investigation of the deformation of weak inclusions regarded as proxies for small-scale ductile s
NU
SC
shear zones. As pointed out earlier, the shear spin w is an entirely different quantity describing the phenomenon exhibited by a deforming ellipsoid whose principal axes swipe through the material as the ellipsoid changes shape (e.g., Goddard and Miller, 1967; Bilby and Kolbuszewski, s
MA
1977). It is called a ‘shear spin’ by Kocks (1998, p.385) in recognition that w arises from shearing between the instantaneous principal axes of the ellipsoid – wijs 0 if ij ( i j ) vanish. The shear spin must be considered for a deformable inclusion. However, the original definition of wijs is singular if the inclusion is a spheroid ( ai a j ). This occurs because a
TE
D
perfect circle does not have definable principal directions in the first place. Jiang (2007b) considers all possible singular cases that may be encountered in a numerical calculation of isotropic inclusion problems. Later, Jiang (2012) simplifies the singularity handling by setting wijs wij when ai a j which is the second equation in (26). Davies et al (2013) have
AC CE P
questioned the justification of this treatment. They instead suggested using a limit for wijs as ai a j based on the assumption that the orientation matrix Q is continuous at the ai a j
state. In the following, the second equation in (26) is further justified and it is shown that Q is intrinsically discontinuous at the ai a j state. The definition of wijs wij when ai a j first given in Jiang (2012) was based on the same line of reasoning as Jiang (2007b): Suppose an inclusion does go through a spheroid state at time t when two semiaxes are equal, say, a1 a2 (Fig.3). Call the partitioned flow field in the inclusion L= + w at time t (Fig.3). Immediately after t , at t t , the sectional shape of the inclusion must represent the infinitesimal strain ellipse in the inclusion. The two semiaxes ( a1 and a2 ) at t t are thus parallel to the 2 eigenvectors of ε , in the limit of t 0 . In an allisotropic system considered by Jiang (2007b), the eigenvectors of ε are the same as the eigenvectors of imposed far-field strain rate tensor. In general, the eigenvectors of ε do not coincide with that of the far-field strain rate tensor, but this does not affect the following argument. Immediately prior to the spheroid state, at t t , the inclusion’s sectional shape must be an ellipse which is reciprocal to the one at t t . This means that, at t t , the twosemiaxes are eigenvectors of ε (Fig.3). Since the eigenvalues of ε and ε are opposite in 19
ACCEPTED MANUSCRIPT sign but equal in magnitudes, it follows that the a1 - and a2 - axes at t t are parallel, respectively, to a2 - and a1 - at t t . In other words, the two semi-axes swap their identities but
PT
do not otherwise change orientations as the inclusion passes through the spheroid state (Fig.3). This behavior is fully captured by setting ij 0 which is equivalent to wijs wij . Numerically,
SC
RI
this treatment amounts to suspending the orientation update (achieved by ij 0 ), as the inclusion goes through a spheroid state, and updating the axial lengths only so that the inclusion’s new principal directions emerge. This fully justifies the second equation in (26) first introduced by Jiang (2012). The assumption that the orientation matrix Q is continuous at the ai a j stage is
NU
invalid. If is the angle between the major eigenvector of ε and one coordinate axis (Fig.3),
MA
cos sin sin cos then lim Qij t t = Qij t t and lim . As the left and t 0 x 0 sin cos cos sin right limits are not equal, Qij t is intrinsically discontinuous at the a1 a2 state.
D
One could make Q continuous if the ordering of a1 a2 is not imposed. In such a case,
TE
the two principal axes merely change relative lengths as the inclusion goes through the spheroid state with a1 a2 before t , a1 a2 at t , and a1 a2 after t . Even in this situation, we still have dQij
0 which is equivalent to setting ij 0 . A great deal of inconvenience arises if the a1 a2
AC CE P
dt
a1 a2 convention is removed. To use elliptic functions to calculate Eshelby tensors for
isotropic materials (Eshelby 1957; Mura 1987, p.78) for instance requires that the three semiaxes be ordered in terms of relative length. In a multi-inclusion system, imposing the a1 a2 a3 condition is the most straightforward way to keep track of the relative length of all inclusions. The intrinsic discontinuity of Qij if an inclusion does go through a spheroid state causes no trouble at all in numerical simulations if the complete definition (26) for shear spin is used. 5.3 Evolution of an inclusion through numerical calculations Given the formulation so far, the evolution of an ellipsoid inclusion in a general viscous material can be calculated numerically for any anisotropic viscous material in a given flow field. The procedure is as follows. The state of the ellipsoid at any time t is defined by its orientation matrix Q (Bunge, 1982; Jiang, 2007a, b) and the three semi-axes are expressed by a vector a1 a a2 . The evolution of the inclusion is found by solving the following two coupled a 3
differential equations (Jiang 2014): 20
ACCEPTED MANUSCRIPT dQ da ˆ , εa ΘQ dt dt
(31)
PT
where ˆij is a diagonal matrix made of the diagonal components of ij . Both ij and Θij are obtained from the above formulation. For a rigid ellipsoid, one only needs to solve the second
RI
equation.
NU
SC
Eqs.31 have been solved numerically by a fourth-order Runge-Kutta method (Jiang, 2007b, 2012), by using the Rodrigues rotation (e.g., Başar and Weichert, 2000, p.31) (Jiang 2013), and by a combination of the two (Jiang 2014). Clearly, it is the evaluation of associated Eshelby tensors in the process to obtain ˆij and Θij that takes the most computational time.
D
MA
Recently, Qu et al. (2016) have developed an optimal algorithm for calculating the Eshelby tensors of highly elongate and/or flat ellipsoids, combining the product Gaussian quadrature used in Jiang (2014) and the Lebedev quadrature.
TE
6. Isotropic systems
AC CE P
While for general anisotropic viscous materials, inclusion solutions must be expressed using integrals (12) and (13), it is shown in this section that, for all-isotropic (both inclusion and matrix are isotropic) systems, explicit analytic expressions can be derived directly from the general formulation. 6.1: Ellipsoid in 3D flows
From (12a and b), one can get the following explicit strain rate and vorticity partitioning equations (see Jiang (2013): ε J d S : C1 : Cinc J d : Ε
(32a)
w W Π : S1 : ε Ε
(32b)
1
inc
where C
stands for viscous stiffness of the inclusion. When the matrix and the inclusion are 1
inc
both isotropic, the C : C term in (32a) reduces to a scalar, r (not to be confused with r in Green functions), standing for the ratio of ellipsoid viscosity to the matrix viscosity and (32a) is simplified to (see Jiang (2013, Eq.11 there): 1
ε isotropic J d (r 1)S : Ε
(33)
21
ACCEPTED MANUSCRIPT All strain rate components in (33) can be solved explicitly (B7 and B8, Appendices B and C) to give:
1 2 r 1 Sijij
RI
, i j
(no sum)
(34b)
NU
where
ij
(34a)
SC
isotropic
ij isotropic
PT
1 1 (r 1) S2222 S2233 11 (r 1) S1122 S1133 22 1 22 isotropic (r 1) S2233 S2211 11 1 (r 1) S1111 S1133 22 33 11 22
11 isotropic
1 (r 1) S1111 S1133 1 (r 1) S2222 S2233 (r 1) 2 S1122 S1133 S2211 S2233 .
This
MA
reproduces Freeman’s (1987) equations with corrections by Jiang (2007b). All viscous Sijkl in (34) can be calculated using Eqs.17 and 18 from explicit expressions for elastic Eshelby tensors given by Eshelby (1957) and Mura (1987).
isotropic
aI2 aJ2 ij aI2 aJ2
(35)
TE
wij
D
In Appendix C, it is shown (Eq.C21) that:
wij
AC CE P
which combined with (34b) gives:
isotropic
ai2 a 2j 2 r 1 Sijij 2 ij (no sum) a a 2 1 2 r 1 S i j ijij
(36)
The shear spin is obtained by combing (26) with (34b) to get: s ij isotropic
w
ai2 a 2j ij , 2 2 a a 1 2 r 1 S j ijij i
(no sum)
(37)
Substituting (36) and (37) to (27) yields the following angular velocity equation for an isotropic ellipsoid embedded in an isotropic matrix: ai2 a 2j ai2 a 2j ij ij Wij 2 r 1 Sijij 2 2 2 2 a a a a 1 2 r 1 S isotropic j i j ijij i
(no sum)
(38)
Jeffery’s result is reproduced by letting r in (38). The angular velocity of an incompressible inviscid ellipsoid in an isotropic matrix is obtained by setting r 0 in (38) to get:
22
ACCEPTED MANUSCRIPT
ij
ai2 a 2j ai2 a 2j Wij 2Sijij 2 1 2Sijij ai2 a 2j ai a 2j ij
isotropic inviscid
(no sum)
(39)
PT
Mulchrone (2007) derives a 2D version of (37) for ellipses in isotropic plane flows.
ijd 2c ij , p 2m iI SiIjJ Jj
(40)
NU
1
SC
isotropic materials, the expressions are as follows:
RI
The deviatoric stresses inside the inclusion can be calculated from its strain rate tensor (obtained above) and the viscous stiffness. The pressure in the inclusion can be calculated using (12c). Because ij is diagonal (16) and Sijkl are zero except for terms with repeated indices for
MA
where the subscripts “c” and “m” stand for, respectively, the inclusion clast and the embedding matrix. As explained earlier, summation applies to repeated lowercase indices while uppercase indices take their lowercase values but are not summed.
TE
D
As already covered in Section 4.2, the mechanical fields outside the inclusion vary with space but can be easily obtained combining the formal solutions (13) with the relations (25) by which all relevant tensors can be calculated from the explicit approach of Ju and Sun (1999) efficiently and accurately.
AC CE P
6.2: Ellipses in 2D flows
An isotropic viscous ellipse in a 2D plane-straining flow can be regarded as an elliptic cylinder in 3D. Elastic Eshelby tensors for isotropic elliptic cylinders are already given in Mura (1987, p.80). With the Poisson’s ratio set to 0.5, the components are: el S1111
2R 1
R 1
2
el , S1122
1
R 1
2
el , S1133
el el el el S3311 0 , S3322 0 , S3333 0 , and S1212
where R
1 R R2 2R R2 el el el , S2211 , , S2233 , S 2222 2 2 R 1 R 1 R 1 R 1 R2 1
2 R 1
(41)
2
a ( a and b being the two semi-axes of the elliptic section with a b ). b
The corresponding non-vanishing viscous Eshelby tensor components can be calculated from (41) using (17) to get: S1111 S2222 S1122 S2211
R
R 1
2
23
(42a)
ACCEPTED MANUSCRIPT R2 1 2 R 1
11
(42b)
2
1 R , 22 R 1 R 1
(42c)
PT
S1212
RI
1 The non-vanishing components of Sijmn for elliptic cylinders under plane-straining flows 1 2D can be obtained using the identity Sijmn [where, as explained in more detail in Smnkl J ijkl
S
1 1212
R 1
1 1122
S
S
1 2211
R 1 4R
2
(43a)
2
2 R 2 1
(43b)
D
S
S
1 2222
MA
1 1111
1 1 ik jl jk il ij kl with i, j, k , l 1,2 ]. They are: 2 2
NU
2D deformations defined as: J ijkl
SC
2D Section 7, J ijkl is the fourth order deviatoric identity tensor for plane-straining incompressible
TE
With (42) and (43), we can express the interior strain rate and vorticity fields of the inclusion in terms of any given far-field plane-straining flow.
AC CE P
If the far-field plane-straining flow is a 2D general shear with vorticity parallel to x3 and the elliptic cylinder is parallel to the x3 axis as well, submitting (42) to (34) and (36) produces the following relations, equivalent to the results of Bilby and Kolbuszewski (1977):
2D 11
2D 22
R 1
2
R 2 2rR 1
11 ,
2D 12
R 1
2
rR 2 2 R r
2D 12
12 , w
1 r R 2 1 rR 2 2 R r
12
(44)
Submitting (42b) into (38) yields the angular velocity of the ellipse (Note for the elliptic cylinder 2D
2D
considered here the non-vanishing components are 12 21 ): 2D 12
isotropic
2 R 1 r R 1 2 R W12 12 2 R 1 rR 2 R r
(45)
which is equivalent to the result of Bilby and Kolbuszewski (1977, their Eq.15). The deviatoric stresses and pressure inside the inclusion can also be calculated similar to 3D cases (40). As ii and Sijkl are given in (42), the pressure in the ellipse, using (12c), is:
24
ACCEPTED MANUSCRIPT
p
2m11 1 r R 2 1
(46)
R 2 2rR 1
PT
which is equivalent to Schmid and Podladchikov (2003, their Eq.67).
NU
6.3: The work of Jaswon and Bhargava (1961)
SC
RI
The exterior fields can be obtained in a manner similar to the 3D case. Jin et al. (2011) give an explicit solution for elliptical cylinders as a special case of Ju and Sun (1999). Their solution can be combined with (25) to calculate the complete mechanical fields outside an inclusion. Another alternative to get the exterior fields is to apply the explicit solutions of Jaswon and Bhargava (1961) for plane-strain isotropic elastic inclusions to viscous conditions as described below.
MA
Analytic solutions for 2D elastic inclusion problems are found by Jaswon and Bhargava (1961). Although they use a complex variable formalism, most of their results are in explicit real-form expressions that can be readily extended to isotropic incompressible viscous materials.
D
The following expressions are from Eqs.19 and 26 of Jaswon and Bhargava (1961) recast in the current notation. They relate the constrained inclusion strain state ijc (Eshelby 1957) and
TE
the eigenstrains ij* :
1 R * * 11* 22 * 22 2 11 R 1 R 1
(47a)
c 22
1 R * * 11* 22 * 22 2 11 R 1 R 1
(47b)
12c
AC CE P
11c
R2 1
R 1
2
12*
(47c)
* Considering the incompressibility condition ( 11* 22 0 ) and using the relation ijc Sijkl kl* (c.f. 1 Eshelby, 1957), (47) immediately lead to all the expressions of Sijkl and Sijkl , given in (42) and
(43) respectively, which lead to the first two equations of (44) for the interior strain rates. Jaswon R 1 * and Bhargava (1961) also give an expression (their Eq.26) for inclusion rotation: w12c 12 R 1 2D R2 1 which readily leads to w12 2 12 and then the interior vorticity (the 3rd equation of 42) R 1 upon substitution of 12 12 E12 . This in turn leads to (45) upon consideration of the shear spin term. 25
ACCEPTED MANUSCRIPT Jaswon and Bhargava (1961) also discuss the limiting cases of a rigid inclusion and a void inclusion which lead to equations later rediscovered by Ghosh and Ramberg (1976) and Mulchrone (2007). The interior pressure can also be obtained from Eq.14a of Jaswon and Bhargava (1961)
PT
R 1
11* . This yields (46)
R 1 which follows from the 11 2
2R
incompressibility condition.
SC
1 1 upon substitution of 11* S1111 S1122 11
2m R 1
RI
which upon consideration of incompressibility becomes p
p ,
MA
NU
Solutions for the exterior fields are also included in Jaswon and Bhargava’s (1961). The expression for the exterior pressure field (their Eqs.13 and 21 combined, using the principle of superposition) is recast below: * 8m R sinh 2 sin 2 12* 1 11 2 R 1 cosh 2 cos 2 cosh 2 cos 2
D
where , are confocal elliptic coordinates which are related to the Cartesian coordinates by:
TE
x1 a 2 b2 cosh cos , and x2 a 2 b2 sinh sin . The eigenstrain terms are, for a 2
AC CE P
heterogeneous inclusion, replaced by an equivalent inclusion with
R 1 2 12 2
* 12
R 1
R 1 11 * 11
2R
and
to yield an explicit expression for the exterior pressure field:
p ,
8m R R 1 11 12 sinh 2 sin 2 1 R 1 cosh 2 cos2 2 R cosh 2 cos2 R 2 1
which can be further recast in the more familiar polar coordinates , , related to the Cartesian coordinates by x1 cos and x2 sin , as:
p , 8m
R R 1 11 12 B , 2 1 A , R 1 2R R 1
where
A ,
F 1 h cos 2 2F 26
(48)
ACCEPTED MANUSCRIPT F 1 h cos 2 2F
F 1 h 2h cos 2 , and h 2
a 2 b2 2
PT
B ,
B , takes positive sign when is in the first and 3 quadrants and negative sign when it is
RI
rd
NU
SC
in the second and fourth quadrants. 11 and 12 in (48) can also be expressed in terms of the far-field strain rates using (44). The expression of (48) is the most succinct real-form explicit expression available for the exterior pressure field.
AC CE P
TE
D
MA
Jaswon and Bhargava’s (1961) work has been overlooked by earth scientists who have, in their investigation of elliptic inclusions in 2D plane-straining flows of isotropic incompressible viscous materials, rediscovered equations available from Jaswon and Bhargava’s work. Ghosh and Ramberg (1976) obtained the angular velocity equation for a rigid ellipse in plane-straining flows of isotropic viscous materials, which is a limiting case explicitly discussed by Jaswon and Bhargava (1961, p. 677). Schmid and Podladchikov (2003) obtained close-form solutions (some expressed in complex variables) for the interior and exterior mechanical fields of a viscous ellipse. Their investigation of inclusions with a mantle layer of dissimilar viscosity is new contribution. Their other results overlap with Jaswon and Bhargava (1961). Mulchrone (2007) derives the angular velocity equation for an inviscid ellipse, which was also a limiting case explicitly considered by Jaswon and Bhargava (1961, p.677).
7. Two dimensional inclusion problems in anisotropic viscous materials 7.1 General statement
There is considerable literature on two-dimensional elastic inclusion problems since Jaswon and Bhargava (1961). Most of the work can be found in the books by Mura (1987) and Ting (1996). The work has been largely motivated by the attempt to avoid the ‘‘formidable” integrals (Yang and Chou, 1976) in the Green-function and point-force approach of Eshelby (Section 3). It has already been shown above that there is no longer any formidable integrals for isotropic viscous materials. The development in Section 4 together with the explicit approach of Ju and Sun (1999) and Jin et al (2011) has circumvented all formidable integrals. Only algebraic calculations and evaluation of two complete elliptic integrals of the first and second kind are necessary when one solves the interior and exterior fields of an ellipsoid of any shape. While it is straightforward to implement the general formulation in this paper numerically to solve any 3D inclusion problems in anisotropic materials, it is computationally intensive especially for the exterior field. Two-dimensional formulations are useful in at least three 27
ACCEPTED MANUSCRIPT
PT
respects. First, they lead to less intensive numerical calculations or even analytical solutions that allow inclusion behaviors to be investigated more thoroughly. Second, certain practical problems like deformation of thin elastic or viscous plates may be approximated by 2D plane-strain or plane-stress problems. Third, two-dimensional solutions, analytic or numerical, can serve as valuable benchmarks for more intensive and complex 3D inclusion systems.
AC CE P
TE
D
MA
NU
SC
RI
Theories for 2D anisotropic elasticity use a complex variable formalism pioneered by Lekhniskii (1950), Muskhelishvili (1953), and Green and Zerna (1968). Jaswon and Bhargava (1961) used the approach to study the 2D inclusion problem in isotropic elastic materials. Their work was extended to anisotropic elastic materials with cubic symmetry by Wallis (1964), with orthotropic symmetry by Bhargava and Radhakrishan (1964) and by Yang and Chou (1976) who gave real-variable expressions (see Mura 1987, p.142). Chen (1967) presented a solution procedure for monoclinic elastic materials based on the formulation of Green and Zerna (1960). Yang and Chou (1976), Hwa and Ting (1989), and Ting (1996) presented solution procedures for elastic inclusions in general anisotropic elastic materials under generalized plane deformation. As pointed out in the Introduction, it is not trivial to extend these works on 2D anisotropic elastic inclusions to incompressible viscous materials. The reason is that the “solutions” in these works are in fact menus, written in complex variable formalism, to be followed to get the solutions from the material’s anisotropic properties and boundary conditions, and it is often not trivial to execute the menu procedure. Take Chen’s (1967) “solution” for an elliptic inclusion in a monoclinic elastic material for example. According to the procedure, one must first solve 2 quartic equations (one for the matrix and one for the inclusion) of complex coefficients which are to be determined from the matrix and inclusion elastic compliances. Second, one solves two linear systems, each with four complex unknowns (one system is presented in Green and Zerna, 1968, p. 324 and the other in Eq.22-25 of Chen 1967). One must also consider possible degenerate cases if any of the quartic equations have repeated roots (Green and Zerna, 1968, p.205). The solutions in Hwa and Ting (1989) and Ting (1996, Chapter 10) require solving sextic equations of complex coefficients and manipulation of many matrices of complex components. Incompressibility may lead to new degenerate situations (Ting, 1996) that demand additional non-trivial treatment not available in the literature. In the following I use the same Green-function approach used in Section 3 for 3D inclusions to develop a 2D inclusion formulation for elliptic inclusions in anisotropic viscous materials. The formulation leads to integral expressions that can be evaluated efficiently in standard mathematic applications. 7.2 Generalized plane flows in anisotropic viscous materials Before proceeding, it is important to point out that a plane-straining flow in which the velocity field is confined to, say the x1 x2 plane [with velocity components being u1 u1 x1, x2 , u2 u2 x1 , x2 , and u3 0 ] is not possible in a general sense in anisotropic materials. This is
28
ACCEPTED MANUSCRIPT because the in-plane motion and anti-plane motion are not always separable in an anisotropic material (Yang and Chou, 1976; Ting, 1996). A flow field in which the anti-plane velocity component u3 does not vanish but nevertheless all velocity components depend on x1 and x2
PT
only [i.e., u1 u1 x1, x2 , u2 u2 x1 , x2 , and u3 u3 x1, x2 ] is called a generalized plane flow. In a generalized plane flow, the shear strain rates 13 and 23 do not vanish although 33 0 .
RI
Generalized plane flows are really special types of 3D flows.
SC
A plane-straining flow can exist in an anisotropic material only if the material has a symmetry plane parallel to the x1 x2 plane. Thus, triclinic materials can never have plane-straining
MA
NU
flows in them. If materials with monoclinic or higher symmetries are oriented such that no symmetry plane is parallel to the symmetry plane of the imposed velocity field, plane-straining flows will not occur in these materials. By inference, the occurrence of monoclinic flows (e.g., Jiang and Williams, 1998; Lin et al., 1998; Passchier, 1998) also requires that the flowing rock material have one symmetry plane parallel to the symmetry plane of the imposed velocity field.
D
The formulation below is valid for plane-straining flows only. Generalized plane flows can be treated using the formulation of Section 3.
TE
7.3 Formulation
The two Green functions, Gij2 D and H i2D respectively, (the superscript “2D” stands for
2D ij
G
AC CE P
the two-dimension case) are the inverse Fourier transforms of (6) in 2D space, respectively:
1
Aij ( z) 1 2 exp iK x dK and H i2 D 2 2 4 K 4
i i exp iK x dK K
Their derivatives, Gik2 D, jl and H i2, Dj , which are used in (10), can be obtained by simple differentiation and are expressed as: Gik2 D, jl
H
2D i, j
1
zzA
4 2
j l
1 4 2
ik
exp iK x dK
(49a)
z i
j
exp iK x dK
(49b)
Submitting (49) into (10), noting that the integrations are over the ellipse area for the 2D case, 2D and carrying out the explicit integration, we obtain the following expressions for ij2 D and Tijkl : 2D Tijkl x
ab
0
Aik z j zl B x, z d
(50a)
29
ACCEPTED MANUSCRIPT ab
0
i z j B x, z d
(50b)
i xz cos 1 1 where z , B r , z 2 2 sin 2 xz
stands for the real part of the complex function in the parentheses. Weinberger et al. (2005)
RI
with 2 a 2 cos2 b2 sin2 , and
PT
ij2 D x
SC
present a detailed procedure for getting (50a) from (49a). Using a similar procedure one obtains (50b) from (49b). For points inside the ellipse (the interior inclusion problem), 2 x z 0 and 1
2
NU
B x, z
2
2D , rendering constant ij2 D and Tijkl as expected. For exterior points (the exterior
TE
D
MA
inclusion problem), the two tensors vary in space. From (50) the corresponding Eshelby tensors for the interior and exterior points can be readily obtained which together with (12) and (13) solve the inclusion problem. The numerical evaluation of (50) can be done using built-in quadrature in any standard mathematics applications such as MAPLE, MATHCAD, or MATLAB, and the computation is efficient as only line-integrals are involved. 7.4 Application to materials with a planar anisotropy
AC CE P
Eqs.50 apply to any monoclinic material with symmetry plane parallel to the x1 x2 plane and subjected to plane-straining flows in the x1 x2 plane. As an illustration of using this 2D formulation, let us consider a viscous material of planar anisotropy specifically. Well-layered or foliated rocks (e.g., Johnson and Fletcher, 1994, p. 419) and minerals with a dominant slip system such as mica (Chen et al., 2014) have been treated as an idealized material with two distinct viscosities: a normal viscosity n that measures the resistance of the material to pure shearing along and perpendicular to the anisotropy plane and a shear viscosity s measuring the resistance to shearing parallel to the anisotropy. Fletcher (2004, 2009) obtains real-form analytic solutions for an isotropic viscous elliptic inclusion enclosed in such a planar anisotropic matrix undergoing plane-straining flows. Numerical calculations are compared with his analytic results. It should be pointed out that heterogeneous strain near the inclusion will lead to distortion of the planar anisotropy around the inclusion (Dabrowski and Schmid, 2011; Griera, et al., 2011, 2013). Such heterogeneity is not captured by the solution here or that of Fletcher (2004, 2009). As discussed above, in order for plane-straining flow to occur, the anisotropy plane (layering) must be parallel to the x3 axis. When the anisotropy plane is perpendicular to the x2 axis, the viscous stiffness tensor components can be written in the following block-matrix form:
30
ACCEPTED MANUSCRIPT 1 2D J 12 kl m 2D J 22 kl
(51)
PT
cijkl
2D J 11kl 2n 1 J 2 D 12 kl m
where the lowercase cijkl is used for 2D viscous stiffness tensor with all indices varying from 1 to
1 1 ik jl jk il ij kl ( i, j, k , l 1,2 ), m n is a measure of the strength of s 2 2
SC
2D J ijkl
RI
2D 2, J ijkl is the fourth order deviatoric identity tensor for plane-straining deformations defined as:
NU
anisotropy (Treagus, 2002; Fletcher 2004, 2009; Chen et al., 2014) which also corresponds to a similar concept in transverse isotropic elastic materials (Brace 1965). In (51) each component is itself a 2 2 matrix.
MA
In the ellipse’s coordinate system at angle relative to the anisotropy plane (Fig.4), the stiffness tensor components can be obtained by tensor transformation, c 'ijkl QiaQ jbQkcQld cabcd ,
TE
D
cos sin where Q . The non-vanishing components of c 'ijkl can be calculated and they sin cos are:
AC CE P
1 c '1111 = c '2222 = c '1122 =n cos2 2 sin 2 2 , m 1 c '1112 = c '1222 n 1 sin 2 cos 2 , m
(52)
1 c '1212 n sin 2 2 cos2 2 m 2D The corresponding compliances can also be obtained using the identity c 'ijmn s 'mnkl J ijkl . The
non-vanishing components are: s '1111 = s '2222 = s '1122 =
s '1112 = s '1222
s '1212
m 1 m 1 cos 4 , 8n
m 1 sin 4 ,
(53)
8n
m 1 m 1 cos 4 8n
31
ACCEPTED MANUSCRIPT which are equivalent to Fletcher (2009, Eq.5 there).
sin 2 2 m cos2 2 , A12 A21 n 1 m sin 4 2 , m 2m
n
m sin 2 2 cos2 2 m
RI
A22
n
(54)
SC
A11
PT
From (52), the corresponding Christoffel stiffness tensor components defined in (5a) can be calculated:
which in turn leads to the following expressions for Aik and i defined in (7) through matrix
n
ij
zi z j
m 1 cos m 1 cos(4 3 ) 2 m 1 sin m 1 sin(4 3 )
D
ζ
m
MA
Aij
NU
inversion:
TE
where m sin 2 2 cos2 2
1
(55a)
(55b)
(55c)
AC CE P
Substituting (55) into (50), the mechanical fields inside and outside the ellipse can be calculated. Below we consider the case in which the remote flow field is a plane-straining general shear with the shear plane parallel to the anisotropy and compare the numerical calculations based on the above with Fletcher’s (2009) analytic solutions.
7.5 Comparison with Fletcher (2004, 2009) Let the pure shearing strain rate and simple shear rate parallel to the anisotropy plane (the x1 x2 coordinate system in Fig.4) be respectively and . In the inclusion coordinate system (the x '1 x '2 coordinate system in Fig.4) at angle relative to anisotropy plane, the remote flow field strain rate and vorticity components are, expressed as functions of , , and :
11 , , 22 , , cos 2 sin 2 2
32
(56a)
ACCEPTED MANUSCRIPT
12 , , 21 , , sin 2 cos 2 2
(56c)
2
PT
W12 W21
(56b)
SC
m R 2 1
r m R 2 1 2 R
E12 , , , m, r
MA
12 , , , m, r
2 mR E11 , , , m, r R 2r mR 1 2
NU
11 , , , m, r
RI
The analytic expressions obtained by Fletcher (2009, his Eqs.10, 14, and 15) for the mechanical fields inside an isotropic elliptic inclusion can be recast into the following forms:
mR
11
2n R 2 1
11 , , r11 , , R 2r mR 1 m
TE
p
n R 2 1
D
m R 2 1 R2 1 w12 , , , m, r 2 12 E 12 , , , m, r R 1 r m R 2 1 2 R
2
(57a)
(57b)
(57c)
(57d)
as:
AC CE P
where r is the ratio of the inclusion viscosity to the matrix n and E ij , , , m, r are defined
E ij , , , m, r rEij , , Eij , , m One can easily confirm that (57) reduce to their isotropic versions in (44) and (46) when m 1 . Fig.5 presents the results for an inclusion ( R 5 , r 10 , and m 25 ) in plane-strain general shear ( 0.5 and 1 ) obtained from numerical calculations using the 2D formulation of this paper implemented in MATHCAD (solid lines) and those calculated directly using (57) (dashed lines). The results agree perfectly. The external mechanical fields can be calculated readily using the formulation here. Fig.6 shows the pressure fields normalized against the matrix deviatoric stress invariant around a strong inclusion ( R 3 , r 100 ). The inclusion long axis is horizontal and the far-field flow is pure shearing with horizontal stretching and vertical shortening. Fig.6a is for an isotropic matrix ( m 1 ). Figs.6b-e are when the matrix is anisotropic with m 10 . The anisotropy plane is horizontal in Fig.6b and increases by 30° increments clockwise until 90° in Fig.6e. It is seen clearly that the pressure perturbations are significantly subdued in an anisotropic matrix than in 33
ACCEPTED MANUSCRIPT an isotropic one (Fig.6a). This has implications for the study of tectonic pressure fluctuations during metamorphism.
PT
7.6: Comparison with 3D inclusion behaviors
NU
SC
RI
Although numerical calculations of 2D inclusion problems are far less intensive than 3D inclusion problems, most natural deformations involving heterogeneous inclusions cannot be treated as two-dimensional problems (see Li and Jiang, 2011 for more discussion). A question that naturally arises is: how well can 2D behaviors be applied to the understanding of 3D inclusion behaviors? For example, from Fletcher’s (2009) analytic solutions, one concludes that the rotational behavior of a rigid ellipse is unaffected by the planar anisotropy in the matrix. This is clear by considering the limit of (57c) as r . It reduces to the Ghosh and Ramberg (1976) relation for rigid inclusion in isotropic matrix. Is it also true for 3D ellipsoids? Is the rotation of a rigid ellipsoid unaffected by the presence of a planar anisotropy in the matrix?
MA
To answer this question, one can examine the general equation in (29) for a rigid ellipsoid 1 in anisotropic viscous matrix. As the term ijmn Smnkl in (29) must certainly depend on the
TE
D
anisotropy of the matrix, it cannot be true in general for rigid inclusion behaviors to be independent of matrix anisotropy. Is it true for the special case of planar anisotropy under the plane-straining flows? We can use the general formulation in Section 3 to answer this question.
Cijkl
AC CE P
In 3D, the viscous stiffnesses of a planar anisotropic material with the anisotropy plane normal to the x2 axis can be expressed in a block matrix form similar to (51):
d J 11kl 1 2n J 12d kl m J 13d kl
1 d J 12 kl m J 22d kl
1 d J 23kl m
J 13d kl 1 d J 23kl m J 33d kl
(58)
d where J ijkl are the 3D 4th-order deviatoric identity tensor as defined in (2).
Let us consider the rotational behavior of a rigid ellipsoid ( a1 : a2 : a3 ) with its a3 semi-axis parallel to the x3 axis embedded in a viscous matrix whose rheology is defined by (58). Suppose the far field flow is a plane-straining general shear like the 2D case considered above. Fig.7 presents the angular velocities (normalized against simple shearing) for spheroids ( a1 a2 ) in a plane-straining flow as functions of matrix anisotropy measured by m . The distinct
a3 axes of the spheroids are parallel to the vorticity of the flow (the x3 axis). The plane-straining flow in the x1 x2 plane is like the 2D case above with 0.5 and 1 along the anisotropy. A 34
ACCEPTED MANUSCRIPT circular cylinder ( 1: 1: ) rotates at a constant rate of half the shearing, independent of m , as (57c) predicts. The angular velocities depend on the relative lengths of a3 and m . A shorter a3
PT
and higher m imply a faster rotation rate. The angular velocity curves vary between the curve for a cylinder ( 1: 1: ) and the one for a disk with diminishing a3 ( 1: 1: 0 ). The same conclusions
RI
can be drawn for non-spherical objects. Fig.8a presents the angular velocities of a rigid triaxial ellipsoid ( 6 : 1: 5 , with the a3 axis parallel to the x3 axis) as a function of the orientation of the ellipsoid long axis with respect to the anisotropy plane ( x1 axis) for different matrix m . A
SC
greater m corresponds to a higher angular velocity at any given orientation, except near the stable orientation of the object. Fig.8b plots the angular velocity variations of the ellipsoid as its a3 axis varies from infinity to 0.01 for a fixed matrix m =10. A relatively shorter a3 generally
NU
corresponds to a higher angular velocity. Fig.9 shows the angular velocities of elliptic cylinder and ellipsoids, all with sectional aspect ratio of 2:1 on the x1 x2 plane but different a3 lengths, in a
MA
matrix with m =10 undergoing simple shearing. The angular velocity increases as a3 varies from infinity to nearly 0.
AC CE P
TE
D
The results above are based on a specific material with planar anisotropy and with the anisotropic plane parallel to the shearing. The behaviors of viscous inclusions in general anisotropic materials are terra incognita. There is indeed a great deal to be learnt using the formulation of this paper. The results from the simple example examined here already suggest that anisotropy affects the behavior of rigid and deformable inclusions significantly. In a 3D case, the ellipsoid can have any initial orientation with principal axes oblique to the far-field flow and the material anisotropy.
8. Applying the inclusion theory to lithosphere deformation So far, we have considered a single elliptic or ellipsoidal inclusion in an infinite homogeneous matrix. To apply the inclusion theory to the deformation of rocks in Earth’s lithosphere requires some assumptions. First, one must address the issue of non-linear rheology of the lithosphere as the inclusion theory is based on linear (viscous or elastic) constitutive behaviors. Second, rheologically distinct elements are irregularly shaped. How a theory based on ellipsoidal and elliptic inclusions may be applied to them also requires justification. Third, rock masses are heterogeneous poly-element aggregates. How and in what sense can they be treated as a “uniform matrix” needs to be clarified. Finally, assumptions like “an infinite matrix” and perfect bonding at inclusion interfaces used in the inclusion theory are not always justified in natural cases. We discuss these issues in more detail below. 8.1: Non-linear rheology and partitioning equations
35
ACCEPTED MANUSCRIPT
(59)
NU
ε(σ) M(sec) (σ) : σ , or σ(ε) C(sec) (ε) : ε
SC
RI
PT
There are no exact solutions like (12) and (13) for nonlinear materials such as the Earth’s ductile lithosphere which is often regarded as power law viscous (Kohlstedt et al., 1995; Tullis, 2002). However, it has been demonstrated that (13) can be used as approximate solutions for nonlinear materials using the linearization approach (Molinari et al., 1987; Lebensohn and Tomé, 1993; Ponte Castañeda 1996; Masson et al., 2000; Molinari 2004). The linearization approach consists in using a linear-form constitutive relation to approximate the actual nonlinear constitutive behavior in the vicinity of a given stress or strain (rate) state. For a power-law viscous material, the nonlinear constitutive equation is often written in a pseudo-linear form as (e.g., Hutchinson, 1976, p.104, his Eq.2.5):
where M(sec) (σ) and C(sec) (ε) are referred to, respectively, as the secant viscous compliance tensor and the secant viscous stiffness tensor (e.g., Lebensohn and Tomé, 1993). C(sec) is
D
MA
equivalent to the effective viscosity (e.g., Ranalli, 1995, p.78). M(sec) and C(sec) are no longer material constants but depend on the current state of deviatoric stress or strain rate. The linearization approach seeks the following linear-form relations: ε M : σ ε0 , or σ C : ε σ0
TE
(60)
to approximate (59). C and M in (60) are respectively the linearized stiffness and compliance
AC CE P
tensors, ε 0 and σ 0 are associated pre-strain rate and pre-stress terms. The two equations in (60) 1
are equivalent if C M , σ0 C : ε0 , and ε0 M : σ0 . A variety of linearization schemes have been proposed, differing in how M and ε 0 (or
C and σ 0 ) are constructed. For example, the tangent linearization uses the first-order Taylor expansion of (59) below:
ε(σ ) ε(σ 0 )
ε ε ε : σ σ0 : σ ε( σ 0 ) : σ0 σ σσ0 σ σσ0 σ σσ0
(61a)
σ( ε ) σ( ε 0 )
σ σ σ : ε ε0 : ε σ( ε 0 ) : ε0 ε ε ε0 ε ε ε0 ε ε ε0
(61b)
where M(tan)
ε σ and C(tan) are the tangent compliance and tangent stiffness σ σσ0 ε ε ε0
tensors respectively. For a power-law viscous material, (61) can be expressed as:
36
ACCEPTED MANUSCRIPT 1 ε(σ) nM(sec) (σ0 ) : σ ε0 , σ(ε) C(sec) (ε 0 ) : ε σ0 n
(62)
SC
RI
PT
1 where n is the power-law stress exponent, and ε0 1 n ε(σ0 ) , σ0 1 σ(ε 0 ) are the n back-extrapolated strain rate and stress terms (Hutchinson 1976; Molinari et al, 1987; Lebensohn and Tomé, 1993). If one uses a linearized stiffness or compliance, such as M(tan) , (12) and (13) are approximate solutions for the nonlinear viscous material. For simple presentation, (12a) is commonly recast in the following form (e.g., Jiang, 2014) with linearized stiffnesses and compliances used:
NU
ε H : σ , σ H : ε
(63)
where H S1 J s : M and H C : S1 J s .
MA
1
The following partitioning equations are obtainable from (63) (c.f. Jiang, 2014, Eq.10
D
there):
ε( k ) A( k ) : E ; σ( k ) B( k ) : Σ
TE
(64)
where A( k ) J d S( k ) : nmC(sec) 1 : C( k )(sec) J d : J d nm 1 S( k ) is the strain-rate 1
AC CE P
partitioning tensor and B( k ) nmC( k )(sec) : A( k ) : C(sec) 1 is the stress partitioning tensor, and nm is the power-law stress exponent of the matrix. The super- or sub- script “k” stands for the k-th element. Alternatively, the partitioning equations are expressed in terms of linearized stiffnesses or compliances (Jiang 2014): ε( k ) A
where A
(k )
(k )
(k )
H
(k )
( k ) 1
H( k ) C
partitioning tensors, A
and β
σ( k ) B
:Eα;
M
(k )
( k ) 1
(k )
(k )
: Σβ
: H( k ) C and α
(k )
(65)
( k ) 1
H( k ) C
being 4th-order and α second order; B
(k )
: Σ0( k ) σ0( k ) are strain rate
H
(k )
: E0( k ) ε0( k ) are stress partitioning tensors, B
37
M
(k )
: H
( k ) 1
(k )
M
being 4th-order and
ACCEPTED MANUSCRIPT β
(k )
second order. Note that the secant-based partitioning tensors ( A( k ) , B( k ) ) are distinct from (k )
(k )
,α ,B
(k )
,β
(k )
).
PT
the linearization-based partitioning tensors ( A
The partitioning equations relate the local deformation in a rheologically distinct element
RI
(RDE, regarded as an ellipsoidal inclusion) to the bulk deformation of the matrix.
SC
8.2: Irregular shape of rheological heterogeneities
AC CE P
TE
D
MA
NU
The simple reason for the common practice of regarding heterogeneous elements as Eshelby ellipsoids in micromechanics (e.g., Mura, 1987; Nemat-Nasser and Hori, 1999; Li and Wang, 2008) is because exact solutions for inclusion problems are limited to ellipsoids and cuboids in an infinite matrix and the Eshelby uniformity is limited only to ellipsoids (Eshelby, 1961; Mura, 1987). Where an inclusion is non-ellipsoidal, the Eshelby solutions (12) have been taken as an approximation for the average fields in the inclusion. It has been hard to justify this practice quantitatively (see Mancktelow, 2011; Jiang 2013 and references therein for some discussions). Budiansky and Mangasarian (1960) appear to be the first to treat individual grains in a polycrystalline material as ellipsoidal inhomogeneities embedded in the polycrystalline medium which is taken to be mechanically homogeneous (Qu and Cherkaoui, 2006). Their treatment has been followed in subsequent wide applications of the Eshelby formalism in micromechanics to obtain the macroscopic effective mechanical properties of a material from the properties of its constituent elements (Mura, 1987, p.421-433; Nemat-Nasser and Hori, 1999, Part I; Molinari, 2002). The practical success of this approach gives confidence on its use. Dislocations, stacking faults, cracks, weakened zones, and other discontinuities have also been treated as Eshelby “ellipsoids” (e.g., Mura, 1987, p.15-20, 240-379; Rudnicki, 1977; Healy et al., 2006; Exner and Dabrowski, 2010). In geology literature, rigid or deformable clasts have commonly been treated as ellipsoidal objects (e.g., Ramsay, 1967, p.209-221; Dunnet, 1969; Gay, 1968; Ghosh and Ramberg, 1976; Ježek et al., 2002). Many physical experiments (e.g., Ferguson, 1979; Arbaret, et al., 2001; Ghosh and Ramberg, 1976) and some theoretical considerations (Bretherton, 1962; Wallis, 1977) find that the rotational behavior of rigid non-ellipsoidal objects is very close to that of their best-fit ellipsoids. The assumption allows modern micromechanics to be applied to the deformation of heterogeneous materials like Earth’s lithosphere. It is also worth pointing out that the shape of an ellipsoid can be varied to fit a wide range of rheologicallydistinct elements. By changing the relative length of its three semi-axes, an ellipsoid varies from rod-like ( a1 , a2 a3 0 ), to spherical ( a1 a2 a3 ), to plane-like ( a1 a2 a3 0 ), and many other shapes. In a non-ellipsoid inclusion, the Eshelby tensor is no longer a constant but an Eshelby Tensor Field (ETF). Zou et al (2010) recently obtain explicit expressions of the ETF and its average for a wide variety of non-elliptical inclusions including arbitrary polygonal ones and 38
ACCEPTED MANUSCRIPT
8.3: Heterogeneous matrix and homogenization
SC
RI
PT
those characterized by the finite Laurent series. They conclude that the ETF inside a nonelliptical inclusion can be quite non-uniform. However, using the elliptical approximation as the average of the ETF is valid for convex non-elliptical inclusions. The approximation becomes less acceptable for non-convex non-elliptical inclusions. Where the shape effect is important, one must take the numerical challenge to address the so-called Eshelby’s problem of non-ellipsoidal inclusions for which there has been progress in isotropic elastic materials (e.g., Zou et al., 2010). In structural geology and tectonics studies, mechanically more competent elements tend to exhibit convex bodies. It is therefore more acceptable to regard Eshelby solutions as the average fields in these elements.
AC CE P
TE
D
MA
NU
In micromechanics, the heterogeneous rock mass composed of rheologically distinctive elements (RDEs hereafter) is replaced by a hypothetical homogeneous equivalent (or effective) medium (HEM). The mechanical interaction between a specific RDE (regarded as an Eshelby ellipsoid) and the surrounding matrix is then replaced by a heterogeneous ellipsoid with a rheologically uniform HEM. The rheological properties of the HEM at a point are taken to be the overall effective properties of all RDEs contained in a representative volume element (RVE) centered at that point. The procedure of getting HEM properties from the properties of constituent RDEs is called homogenization. Most homogenization methods are described in the books by Mura (1987), Nemat-Nasser and Hori (1999), and reviews by Zaoui (2002), Molinari (2002), and Mercier and Molinari (2004). Jiang (2014) summarizes the self-consistent homogenization method for poly-element continua based on Lebensohn et al (2003, 2004). A more complete summary of self-consistent homogenization is given below briefly. When the deforming rock mass is made of a more or less homogeneous matrix containing dilute heterogeneous inclusions, the contribution of the inclusions to the bulk rheology is negligible because the interactions among inclusions are insignificant. The HEM property is essentially the matrix property and the Eshelby solutions (63) and (64) can be used directly to find the mechanical fields inside individual inclusions (using the interior solutions) or to get the mechanical fields surrounding an inclusion (using the exterior solutions). As the inclusion volume fraction increases, interactions among inclusions become significant (e.g., Ildefonse et al., 1992, Mandal et al., 2003). A self-consistent homogenization is necessary to obtain the effective properties of the HEM, which is conducted by considering a representative volume element (RVE) of the material sufficiently large so that it statistically represents the material. It is also tacitly assumed that the constituent RDEs in the RVE are uniformly mixed and have a fractal-like structure in shape and orientations regardless of sizes. A salient reference medium is defined. Where the rock mass has a connected matrix containing heterogeneous inclusions, the Mori-Tanaka (1973) approach can be adopted. For a heterogeneous poly-element aggregate without a connected matrix, the reference medium is the HEM. In either case, each element is regarded as an ellipsoidal inclusion surrounded by the 39
ACCEPTED MANUSCRIPT reference medium and subjected to remote loading ( Em for strain rate and Σm for stress) which are, respectively, the average strain rate and associated deviatoric stress tensors in the reference material. The self-consistent homogenization works as follows. Consider a representative volume element (RVE) made of N different constituent elements and label them by k (k=1, 2, … N). The
where H
: σ( k ) Σm
(66)
RI
(k )
(k )
SC
ε( k ) Em H
PT
average strain rate ε ( k ) and stress σ ( k ) in a given element is given by (63):
is evaluated at Em .
NU
The strain rate partitioning equations follow from (64) and (65) as:
ε( k ) A
(k )
: Em α
(67a)
MA
ε( k ) A( k ) : Em (k )
(67b) means volume average over
D
Imposing the consistency condition ε ( k ) E on (65), where
1
:E
AC CE P
Em A( k )
TE
the RVE, gives:
Em A
(k )
1
:E A
(k )
1
: α
(68a)
(k )
(68b)
Combining the nonlinear response of each element (59), σ( k ) C( k )(sec) : ε( k ) , with (67a) and (68a) yields σ( k ) C( k )(sec) : A( k ) : A( k )
1
: E . Imposing the stress consistency condition
σ( k ) Σ on this equation gives the following expression for the homogenized secant viscosity: (sec)
C
C( k )(sec) : A( k ) : A( k )
1
(69) (k )
Combining the linearized response of each element (48), σ( k ) C
: ε( k ) σ0( k ) , with
(67b) and (68b) and then imposing the stress consistency condition σ( k ) Σ gives the following expressions for the linear-form of homogenized stiffness and back-extrapolated term: C C
(k )
Σ0 C
:A
(k )
(k )
:α
(k )
: A σ0
(k )
(k )
1
C
(k )
:A
(k )
: A
40
(k )
(70)
1
: α
(k )
ACCEPTED MANUSCRIPT A similar procedure can be followed with the stress partitioning equations from (64) and (65) to get the following homogenized compliances and the back-extrapolated term: M
M M
(k )
:B :β
(k )
(k )
: B( k ) : B( k )
: B
(k )
1
(71a)
1
ε 0( k ) M
(k )
:B
(k )
: B
(k )
1
: β
(71b)
(k )
SC
E0 M
(k )
( k )(sec)
PT
(sec)
RI
M
NU
Eqs.71b are the homogenization equations used in the self-consistent MOPLA of Jiang (2014). An initial set of guesses for the macroscale stiffness and compliance quantities are made to calculate them iteratively (Jiang 2014).
relations [like M nM
(sec)
MA
For a poly-element material, its secant properties, C
(sec)
and M
(sec)
, do not have simple
] with the linearized properties, C , M , Σ 0 , and
E0 that one may
AC CE P
TE
D
expect of mono-phase polycrystals (Hutchinson, 1976; Lebensohn and Tomé 1993). The rheological response of a general deforming heterogeneous material varies with time because fabric develops in the material as strain increases. Such a material is not expected to have a constant stress exponent. The common use of a steady-state, always-isotropic power-law rheology to represent Earth’s lithosphere in geodynamic models is likely unrealistic. Cook et al. (2015) arrived at a similar conclusion . 8.4: Inclusions in a finite space
To apply the inclusion theory to lithospheric deformation, the conditions of scale separation (Zaoui 2002; Jiang 2014) must be satisfied. To state simply, this requires that the size of an average RDE be small compared to the size of the RVE which itself must be small compared to the scale of the deformation. The scale separation conditions justify the use of the inclusion theory developed for an infinite body. Where the size of the inclusion phase becomes appreciable compared to its distance from the deformation boundary, boundary effects have to be considered. Both translational and angular velocities of the inclusion are affected (e.g., Goldsmith and Mason, 1967). This constitutes inclusion problems in a finite space that has been considered by Kinoshita and Mura (1984). Existing work is still limited to the spherical shape. Li et al. (2007) obtained closed-form solutions for the stress fields of a spherical inclusion concentrically placed in a finite spherical space. Zou et al. (2012) proposed a superposition framework to study inclusions in a finite space. Further progress to these finite space problems would enable better estimation of the mechanical properties for actual materials.
41
ACCEPTED MANUSCRIPT On the other hand, analogue experiments and numerical modeling have improved understanding of the rotational behaviors and surrounding flow patterns of rigid inclusions in confined simple shear (e.g., Marques et al., 2005a, b).
PT
8.5: Interface properties
AC CE P
TE
D
MA
NU
SC
RI
It is well known that the mechanical interactions between the inclusion and the embedding matrix and therefore the overall effective properties of heterogeneous materials are significantly influenced by the interface properties between the constituents. The interfaces may represent weak materials due to one or a combination of changes in mineral phase and texture, metamorphic reaction zones, and fluids. Consideration of interface properties constitutes problems of coated inclusions (Buryachenko, 2007, p.78-79). A weak interface may be due to the presence of a low-viscosity rim (Schmid and Podlachikov, 2003), creation of Somigliana’s dislocations (Mura, 1987, p.484-492) or other mechanisms of interface failure (Hashin, 1991; Huang et al., 1993; Lubarda and Markenscoff 1998). In any of these cases, the shear stress on the inclusion interface is relaxed. The simplest end-member situation is a perfectly-slippery interface which supports no shear stress. This implies that the inclusion, rigid or deformable, supports no deviatoric stresses and must behave like an inviscid inclusion (Xu et al., 2011). Eq.30c describes the instantaneous angular velocity of such an inclusion. For a perfectly-slippery inclusion, it is important to consider the shear spin effect, as in (30c), although the inclusion is ‘rigid’ in the sense that it does not change shape (because of vanishing deviatoric stresses). This is because a freely-slipping interface by definition requires that the hoop particles in contact with the inclusion slide freely (Fig.10) which is equivalent to allowing the principal axes of the inclusion to swipe away from their current positions. Failure to consider the shear spin contradicts the slippery interface assumption. There is good evidence from analogue experiments and natural observations that slippery interface may be significant in natural deformations (Ildeonse and Mancktelow, 1993; Arbaret et al., 2001; Marques and Coelho, 2001; Pennacchioni et al., 2001; Ceriani et al., 2003; Marques and Bose, 2004; Schmid and Podladchikov, 2004, 2005; Marques, et al., 2005). Eq.30c can be used to describe the angular velocity of a perfectly-slippery ellipsoid in a general anisotropic viscous material and Eq.39 is for the special situation where the matrix is isotropic. Mulchrone (2007) derived a 2D version of (39). Johnson et al. (2009) have validated the 2D equation numerically. The work of Johnson et al. (2009) confirms the importance of shear spin as discussed above.
9. Discussion 9.1 Mechanical interactions is key to natural deformation
42
ACCEPTED MANUSCRIPT
TE
D
MA
NU
SC
RI
PT
In essence, the inclusion theory is about the mechanical interactions between a heterogeneous element (an Eshelby inhomogeneity) and the embedding material (the hypothetical uniform HEM). In fact, the Eshelby-type interior solutions (12) have often been referred to as “interaction laws” (Molinari, 2002; Molinari and Mercier, 2004). To be complete, the interaction laws should also include the exterior solutions (13). Such mechanical interactions are indeed the key to understanding natural deformation and the accompanying fabric development. At the structural geology scale, the behaviors of deformable and rigid elements, which give rise to the development of shape fabrics, are manifestation of the mechanical interaction of the elements with the embedding material. The resulting interior fields in an element are responsible for the development of the fabrics inside the element. The development of tails and other geometric patterns around a heterogeneous element (e.g., Passchier and Trouw, 1996) results from the exterior flow fields and can be investigated using exterior solutions (13). The kinematics and mechanics in a Ramsay and Graham (1970) type ductile shear zone may be treated as interior problems of a flattened weak ellipsoidal inclusion (e.g., Xiang and Jiang, 2013) embedded in a country rock mass. The mechanical interactions governed by the interaction laws (12) and (13) determine the deformation history of the shear zone (interior) as well as the manner in which the active zone may increase or decrease in volume (White et al., 1980) during deformation. It is not appropriate to separate the deformation within the zone from that of the “country rock”.
AC CE P
The significance of mechanical interactions extends to all scales of natural deformation. A regional geological unit such as a plutonic body or a continental block may be regarded as a rheologically heterogeneous inclusion in a broader orogenic belt undergoing large-scale deformation modulated by plate tectonics. The deformation and related fabric development within the body constitute an interior problem. The deformation surrounding the body forms an exterior problem. In this context, one considers the orogenic system as a whole subjected to a far-field loading or displacement condition rather than treating the deformation of a geological body in isolation with ad hoc imposed boundary conditions. In such a whole system approach, the deformation of an individual element is the consequence of partitioned deformation due to rheological interactions between the element and the whole deforming system. The partitioning varies with time as the configuration of the deforming system varies continually. To capture this, partitioning (interaction laws) and homogenization equations must be considered simultaneously in a self-consistent way (Jiang 2014). Treating the deformation of an element in isolation, such as a ductile shear zone under a prescribed constant boundary condition, fails to consider the dynamic evolution of the whole deforming system and mechanical interactions among rheologically distinct geological elements. 9.2 A full mechanical and multiscale approach to natural deformation Studies of natural deformation involving large finite strains have been an important subject in structural geology and tectonics (e.g., Ramsay 1967; Hobbs et al., 1976; Passchier and Trouw, 1996). Traditionally the focus has been based on kinematic analysis. In short, the 43
ACCEPTED MANUSCRIPT
MA
NU
SC
RI
PT
kinematic approach investigates how lines, planes, and 3D elements rotate and change shape in a given flow field and, how all this relates to fabric development. The choice of the flow field is left to the intuition of the investigator. While on the plate scale, one may reasonably define the movement (the velocity field) across plate boundary regions because relative plate motions can be steady for many million years (e.g., DeMets and Wilson, 1997; DeMets and Traylen, 2000; Bird, 2003), how the encompassing rheologically heterogeneous lithosphere (made of multiscale elements) actually accommodates the imposed boundary motion, on progressively smaller scales, is governed by mechanical interactions among the constituent elements. It is impractical for one to assign a realistic flow field for the deformation within a rheological element as the latter is a partitioned field due to the interactions. The partitioning varies with time because of the continual change of the configurations of the elements and hence their mechanical interactions. The inclusion theory and associated micromechanics provide a powerful and rigorous mean to tackle such mechanical interactions and continual partitioning. It allows the kinematic aspects (strain rates, vorticity, finite strains) and dynamics of deformation (deviatoric stresses and pressure) to be derived completely from the solutions (12) and (13) rather than from simplistic kinematic models with ad hoc boundary conditions.
AC CE P
TE
D
Furthermore, the partitioning and homogenization equations are tools to bridge the scale gaps. The inclusion theory can be applied on multiple levels. Lower-order (usually larger) inclusions can contain higher-order (usually smaller in scale) inclusions. The partitioned mechanical fields in a lower-order inclusion can be further partitioned to give mechanical fields in higher-order inclusions contained in the lower-order one. This approach can thus address multi-scale deformations (Jiang and Bentley, 2012; Jiang 2014). In the self-consistent MultiOrder Power-Law Approach (MOPLA, Jiang and Bentley 2012; Jiang 2014), two levels of partitioning were considered. A macroscale deformation is first partitioned self-consistently into the constituent rheologically distinct elements making up the macroscale deforming rock masses. The partitioned flow fields in individual heterogeneous elements are used to investigate fabric development, on a smaller scale, in the elements. Through this approach, observations made on smaller scales are related to deformation boundary conditions at larger scales.
10. Summary and conclusions The long history, continued interest, and recent advances in inclusion studies are testament of the enormous significance of the problem. The theory of viscous inclusions is the backbone for the visco-plastic self-consistent model for texture development in polycrystalline materials (Lebensohn and Tomé, 1993), the self-consistent Multi-Order Power Law Approach (MOPLA) for ductile deformation of Earth’s heterogeneous lithosphere (Jiang, 2014), and many self-consistent homogenization models to estimate the overall rheological properties of heterogeneous materials (Molinari, 2002). The considerable literature over many disciplinary fields on the problem is, however, a challenge to earth scientists. As pointed out in this paper, 44
ACCEPTED MANUSCRIPT
RI
PT
duplicate investigations have occurred. In an effort to summarize and advance the current status on viscous inclusion studies, the powerful Green function method has been used in this paper to unify the mathematic notations and presentation. In order to address viscous incompressibility explicitly, two Green functions are used, following Molinari et al. (1987). This necessitates deriving expressions for 3D Green functions and their derivatives for incompressible materials using Barnett’s (1973) method. More explicit solutions than Molinari et al. (1987) are obtained for both the interior fields, expressed in (12), and exterior-point fields, expressed in (13). The general solutions are expressed in terms of Eshelby tensors and an auxiliary tensor
NU
SC
for the pressure field (the interior S , Π , Λ and exterior S E , Π E , Λ E ). It is straightforward to develop numerical algorithms to evaluate these quantities using standard mathematics applications. For isotropic viscous materials, evaluations of all these tensors are simpler, more
MA
efficient and more accurate with the advances made in this study. The evaluations of S E and Π E remain computationally intensive for general anisotropic materials. Fortunately, for the selfconsistent MOPLA and homogenization models, only interior solutions are needed.
TE
D
For isotropic viscous inclusion problems, many mathematic identities among viscous tensor quantities are discovered in this paper that greatly simplify the solutions and their numerical evaluations. These identities enable the use of existing analytic expressions for isotropic elastic interior-point Eshelby tensors to calculate the interior-point viscous Eshelby tensors ( S , Π ) and Λ . The identities also enable the use of the explicit method of Ju and Sun
AC CE P
(1999) for exterior-point elastic Eshelby tensor to calculate viscous exterior-point S E , Π E , and Λ E . As a result, numerical solutions for isotropic inclusion problems can be obtained much more efficiently and accurately than before. The most general equation is established for the rotation of a deformable anisotropic 1 kl wijs , Eq.28). ellipsoid in a rheologically dissimilar anisotropic matrix ( ij Wij ijmn Smnkl 1 The particularization of this equation for rigid ellipsoids leads to ( ij Wij ijmn Smnkl Ekl , Eq.29)
which is the generalization of Jeffery’s (1922) equation to an arbitrary anisotropic viscous matrix. The rotational behavior of a slippery ellipsoid is also obtained 1 a 2 aJ2 d d [ ij Wij ijkl I2 J ijkl J S mn , Eq.30c]. 2 klmn aI aJ
The shear spin of a deformable ellipsoid describes the swiping-through-material phenomenon of the ellipsoid’s principal axes. This phenomenon has nothing to do with material rotation of the ellipsoid which is described by the vorticity in the ellipsoid (the partitioned field). This helps to clarify and justify Jiang’s (2012) treatment for shear spin in the singular event where a deformable ellipsoid passes through a spheroid or sphere state.
45
ACCEPTED MANUSCRIPT
PT
In applying the general solutions (12) and (13) to isotropic inclusions, it is shown that all existing results that have been developed using different approaches in the literature can be derived directly from the general formulation without undue difficulty. This again simplifies inclusion research. Isotropic solutions can serve as benchmarks for numerical anisotropic solutions, and together they can be used to benchmarks more complex numerical models.
MA
NU
SC
RI
There has been significant progress on 2D inclusion studies in anisotropic elasticity. The general complete solutions are now available for generalized plane deformation (Ting, 1996). This body of work adopts the complex variable formalism. Because of the nature of the formalism and resulting generally implicit expressions, it remains open how or whether the results can be applied to incompressible viscous anisotropic materials. Considering this and the fact that numerical evaluation of 3D inclusion problems in anisotropic materials is computationally intensive, a 2D formulation is developed, using the same Green function method that was used in 3D, specifically for anisotropic viscous elliptic inclusions in a rheologically dissimilar anisotropic viscous matrix that can host plane-straining flows. Numerical calculations based on this 2D formulation agree perfectly with the analytic solutions of Fletcher (2009).
AC CE P
TE
D
A reconnaissance investigation on elliptic cylinders and ellipsoids in a matrix of planar anisotropy undergoing plane-straining flows shows that elliptic cylinders behave quite differently from 3D ellipsoids. This implies that conclusions drawn from 2D inclusion studies may not apply to 3D situations for anisotropic materials. As rocks are generally anisotropic, one should always consider the 3D effect of deformations where heterogeneous inclusions are investigated. More detailed investigation on this is under way. The application of the viscous inclusion theory based on linear rheology to the rheologically non-linear lithosphere is through the linearization approach. The heterogeneous lithosphere is approximated by a homogeneous equivalent matrix whose effective rheology is obtained through a micromechanical homogenization method. A heterogeneous rheological element is regarded as an ellipsoidal inclusion which interacts with the homogeneous equivalent matrix. The conditions of scale separation must also be satisfied in order for a finite volume deformation problem to be regarded as one in an infinite body. The interface between the inclusion and the matrix is assumed to be coherent. These conditions and assumptions may be justified from certain problems under given circumstances. Where they cannot be justified, modification of the theory is necessary to account for, such as irregular inclusion shapes, incoherent interface, and boundary effects of deformation. The viscous inclusion theory with its extension to rheologically heterogeneous and nonlinear materials forms the backbone theory to address the ubiquitous mechanical interactions on all scales in natural rock deformations. It enables a full mechanical and multiscale approach toward lithospheric deformation and the accompanying texture and fabric development.
46
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Acknowledgements: The author thanks the Journal editor and three reviewers for their critical comments and acknowledges financial support for the research from Natural Science and Engineering Research Council of Canada (NSERC) through a Discovery Grant, the China National Natural Science Foundation (grant NSFC41472184), a Distinguished Visiting Professorship from China Shaanxi Province and Northwest University as well as internal funding from the State Key Laboratory for Continental Dynamics at Northwest University.
47
ACCEPTED MANUSCRIPT Appendices A: Expressions of Gij , Gij ,l , and H i
PT
Integral expressions for Gij and Gij ,l for anisotropic elastic materials are given by
RI
Barnett (1972) and Mura (1987, p.33-34). Barnet’s derivations are repeated here in order to get expressions of Gij , Gij ,l , and H i for incompressible viscous materials.
SC
The two Greens functions in real space are inverse Fourier transforms of the expressions in (6): Aij (z) exp iK x dK 8 K 2
Hi x
i i exp iK x dK 8 K
3
3
MA
1
(A1)
(A2)
D
K and the integrations are over infinite 3D space. K
TE
where z
1
NU
Gij x
With Euler’s formula, the expressions become:
1
Aij (z) cos Kz x dK 8 K 2
AC CE P
Gij x
Hi x
3
1
8
3
(A3)
i
K sin Kz x dK
(A4)
Here only the real parts of the integrals are concerned as Gij and H i are real. The expression for Gij , l can be obtained by differentiating (A3) with respect to xl :
Gij ,l x
1
zl Aij (z) sin Kz x dK 8 K 3
(A5)
With the following set of variable substitution: x
x k dk , K , and therefore dK 3 r r r
(A6)
(A3-5) are reduced to:
48
ACCEPTED MANUSCRIPT
1 8 r
i
3 2
Gij ,l x
Aij (z) cos kz x dk k2
r
k
sin kz x dk
1
(A7a)
PT
8
3
zl Aij (z) sin kz x dk 8 r k
3 2
RI
Hi x
1
SC
Gij x
(A7b)
(A7c)
A spherical coordinate system aligned with x can be used to evaluate these integrals:
NU
d 3k k 2 sin dkd d , z x cos
(A8)
MA
where is a polar angle in the plane perpendicular to x (i.e., z x 0 ). With (A8), the three integrals in (A7) become:
Gij x 3 Aij (z)sin cos k cos dk d d 8 r 0 0 0
(A9a)
H i x 3 2 i sin k sin k cos dk d d 8 r 0 0 0
(A9b)
Gij ,l x 3 2 zl Aij (z)sin k sin k cos dk d d 8 r 0 0 0
(A9c)
2
TE
D
1
2
AC CE P
1
1
2
Making use of the following identities (e.g., Jones, 1966): cos k cos dk cos , 0
k sin k cos dk sin cos , and
f ( x ) ( n ) x dx ( 1)n f ( n ) (0) , we get from (A9):
0
Gij x
2
1 8 2 r
Hi x
A (z ) ij
0
1 8 r
2 2
2
i
0
d
(A10a)
2
d
(A10b)
2
49
ACCEPTED MANUSCRIPT
z A ij
l
0
More explicit expressions for
i
d
(A10c)
2
and
zl Aij
2
2
SC
First, z x cos implies (Barnett, 1972):
zi xi
2
Aij Aij z l Aij zl x l Aij zl
(A12)
MA
zl Aij
(A11)
NU
2
thus
are derived below.
PT
8 r
2 2
2
RI
Gij ,l x
1
2
D
i The expressions for Aij and are obtained by differentiating the identity
AC CE P
TE
A z A ζ I and reorganizing: T T z 0ζ
Aij A Aim A jn mn Aim j A jm i x m
(A13a)
i A Ain mn m Aim i m x m
(A13b)
Using the definition of the Christoffel stiffness matrix in (5a), we have A
2
z p zq C p q z p zq C p q z p x q zq x p
(A14)
2
Calling C p q z p x q zq x p Fi
A
and inserting it into (A13), we get:
i Ain mn m Aim i m x m
(A15a)
50
ACCEPTED MANUSCRIPT M ij
Aij Aim A jn mn Aim j A jm i x m
(A15b)
8 2 r 2
F d i
0
2
1 8 2r 2
x A
ij
l
zl M ij d
0
RI
Gij ,l (x)
2
1
SC
H i ( x)
PT
With (A15), (A10b) and (A10c) become respectively: (A16a)
(A16b)
NU
Finally, from (5a), it is clear that Aij z Aij z . From (7) it is easy to show that the following relations are true:
MA
1 1 zi i 1 , Aij Aim jm j zm Ajm im i zm , and Aij i z j
(A17)
D
Therefore, we have Aij (z) Aij (z) , (z) (z) , and i (z) i (z) . Considering this in
TE
(A15), it is clear that the integrands in (A10a) and (A16) are all even with respect to z. Thus Gij ,
AC CE P
Gij , l , and H i can be written as:
Gij (x)
1
4 2 r 0
Aij d
(A18a)
1 H i (x) 2 2 Fd i 4 r 0 Gij ,l (x)
1
4 2r 2 0
x A l
ij
(A18b)
zl M ij d
(A18c)
For isotropic materials, (A17) becomes:
i zi , Aij
1 z z , and ij i j
(A19)
B: Relations among tensor quantities for isotropic incompressible materials 51
ACCEPTED MANUSCRIPT Consider the integral of (14b) for isotropic incompressible viscous materials with
1
PT
s Cijkl 2 J ijkl , where is the viscosity. It can be shown that Aik ik zi zk , which turns
(14b) to:
ik
zi zk z j zl 3dS
RI
a1a2 a3 4
z 1
(B1)
SC
Tijkl
2 ij kl , ( and being respectively 1 2
NU
where is the same as defined in (14).
MA
el s For isotropic elastic materials, Cijkl 2 J ijkl +
the shear modulus and Poisson’s ratio, e.g., Mura 1987). We have Aij ij
el Tijkl
a1a2 a3 4
TE
D
1 1 zi z j . The interaction tensor is: ij 2 1 1 zi zk z j zl 3dS ik 2 1 z 1
AC CE P
Aij1
1 zi z j and 1 2
s The symmetric part of (B1) and (B2), Tijkl
(s) Tijkl
a1a2 a3 4
el(s) Tijkl
a1a2 a3 4
(B2)
1 Tijkl Tjikl Tijlk Tjilk , are respectively: 4
1 3 4 ik z j zl jk zi zl il z j zk jl zi zk zi z j zk zl dS z 1
(B3a)
1 1 zi z j zk zl 3dS (B3b) ik z j zl jk zi zl il z j zk jl zi zk 2 1 z 1 4
Therefore, we have from (B3):
Tijkl( s ) Tijklel(s)
(B4)
=0.5
The viscous Eshelby tensor and the elastic Eshelby tensor are respectively:
52
ACCEPTED MANUSCRIPT (s) Sijkl 2Tijkl
(B5) 2 el(s) Tij kl 1 2
(B6)
PT
el el(s) Sijkl 2Tijkl
Tijel(s)
a1a2 a3 1 2
zz 8 1 i
3
j
dS
1
el_s ij kl 2Tijkl
which yields: 0.5
el(s) ij kl 2 Tijkl
0.5
(B8)
D
el Sijkl
MA
NU
Submitting (B7) to (B6) reduces it to: el Sijkl
(B7)
SC
ξ 1
1 2 ij 2 1
RI
From (B3b) and (11), one gets:
TE
Considering (B4) and replacing the shear modulus by viscosity in (B8), the second term on the
AC CE P
right of (B8) becomes the viscous Sijkl because of (B5) and we have: el Sijkl ij kl Sijkl
(B9)
0.5
el Furthermore, because Sijkl are zero except for terms with repeated indices and because
Sii11 Sii 22 Sii 33 0 due to incompressibility (Jiang 2013), we get the following two relations
from (B9): ii
1 el Sii11 Siiel22 Siiel33 0.5 3
el Siijj ii Siijj
0.5
(no sum)
(no sum)
(B10) (B11)
which are respectively (17a) and (17b) in the text.
53
ACCEPTED MANUSCRIPT a Taking the anti-symmetric part of (B1) and (B2), Tijkl
1 Tijkl Tjikl Tijlk Tjilk , and 4
a1a2 a3 16
ik
z j zl jk zi zl il z j zk jl zi zk 3dS
z 1
RI
(a) ( a _ el ) Tijkl Tijkl
PT
a then contract Tijkl with the corresponding stiffness, it is easy to demonstrate that
ik
z j zl jk zi zl il z j zk jl zi zk 3dS
z 1
(B12)
NU
a1a2 a3 8
( el ) ijkl ijkl
SC
which is equivalent to:
MA
This is (18) in the text.
C: Derivation of equations for incompressible isotropic materials
ai2 a 2j ai2 a 2j
S ijij
(no sum)
(C1)
TE
ijij
D
First, we prove that the following relation holds for isotropic materials:
Siiii
AC CE P
Proof: Components of Eshelby tensors for isotropic incompressible materials are zero except for the following ones with repeated indices [combining Eshelby, 1957, (17) and (18) in the text]: 3 2 ai 4
ii
ii , Siijj
ijij ijji
j
8
i
3 2 aj 4
ij
ii , Sijij
3 ai 2 a j 2 8
, jiij jiji ijij
ij
,
(C2a)
(C2b)
where, as in the rest of this Appendix, repeated indices do not imply summation unless specified otherwise and the
terms are defined by the following integrals (Jeffery, 1922; Eshelby, 1957):
i
du 2 a1a2 a3 2 , 0 ai u
ii
2 a1a2 a3 0
a
2 i
du
u 2
, and
ij
2 du a1a2 a3 2 , ( i j) 2 3 0 ai u a j u
54
(C3)
ACCEPTED MANUSCRIPT a
where
2 1
u a22 u a32 u .
1 j i , ( i j , ai a j ) 3 ai2 a 2j
SC
Combining (C4) and the 3rd expression in (C2a) gives (C1).
(C4)
RI
ij
PT
The following relation was given by Eshelby (1957, Eq.3.13 there) but also follows directly from C3:
1 4Sijij
MA
1 Sijij
NU
Second, we prove the following relation which will be used in the derivation of many expressions:
d Proof: Because there are only 4 non-vanishing components of J ijkl if
(C5)
i j and they are
sum on lowercase indices only
TE
1 d SIJkl S klIJ J IJIJ 0.5
D
d d J ijij J djiij J ijji J djiji 0.5 . We therefore have by definition:
(C6)
AC CE P
Spelling out all nonzero terms in the contraction, the left side of (C6) is: 1 1 1 1 SIJkl SklIJ SIJIJ SIJIJ SIJJI S JIIJ 2SIJIJ SIJIJ
(no sum)
(C7)
Combining (C6) and (C7) gives (C5). Third, we derive (32) in the text. The tensor form equation (31) is rewritten below omitting the label: 1
ε J d (r 1)S : Ε
(C8)
For the three diagonal components of the strain rates ii (i 1, 2,3) , we expand (C8):
II (r 1)SIIkl kl II
(sum on repeated lowercase indices only)
Expanding (C9) and using the incompressible condition lead to the following linear set:
55
(C9)
ACCEPTED MANUSCRIPT 11 (r 1) S111111 S1122 22 S1133 33 11 22 (r 1) S221111 S2222 22 S2233 33 22
(C10)
PT
33 (r 1) S331111 S3322 22 S3333 33 33 11 22 33 11 22 33 0
RI
which can be reduced to two equations with two unknowns ( 11 and 22 ) written in the following
SC
matrix form:
(C11)
NU
1 (r 1) S1111 S1133 (r 1) S1122 S1133 11 11 (r 1) S2211 S2233 1 (r 1) S2222 S2233 22 22
MA
The solution to this system gives the following expressions (see Freeman, 1987 with errors corrected in Jiang 2007):
(C12)
AC CE P
TE
D
1 1 (r 1) S2222 S2233 11 (r 1) S1122 S1133 22 1 22 (r 1) S2233 S2211 11 1 (r 1) S1111 S1133 22 33 11 22
11
where 1 (r 1) S1111 S1133 1 (r 1) S2222 S2233 (r 1)2 S1122 S1133 S2211 S2233 . This is (32a) in the text.
To get (32b) for the off-diagonal terms of ij ( i j ), we proceed as follows. First by denoting Ξ J d (r 1)S , (C8) is rewritten as:
ε Ξ1 : Ε
(C13)
It follows readily that:
ijij J ijijd (r 1)Sijij 0.5 (r 1)Sijij , i j
(C14)
Applying the same argument and procedure that (C5) was proven to (C14), we get: 1 ijij
1 4 0.5 (r 1) Sijij
(C15)
Expanding (C13) and making use of (C15), we get: 56
ACCEPTED MANUSCRIPT 1 1 1 1 IJ IJkl kl IJIJ IJ IJJI JI 2 IJIJ IJ
IJ , IJ 1 2(r 1) S IJIJ
(C16)
PT
which is (32b) in the text. Fourth, we derive (33) below.
IJIJ
1 aI2 aJ2 4 aI2 aJ2
SC
S
1 IJIJ
RI
Combining (C1) and (C5) leads to the following relation:
(C18)
MA
1 wij ijkl Sklmn mn
NU
Now, write (12b) in index form:
(C17)
TE
to the kl indices, (C18) becomes:
D
Because there are only 2 nonzero terms for ijkl (i.e., ijij ijji ), after contraction with respect
AC CE P
1 wIJ 2 IJIJ SIJmn mn (sum on repeated lowercase indices only)
(C19)
Because there are only 2 nonzero terms of S IJmn for isotropic materials ( S IJIJ S IJJI ), after contraction with respect to the mn indices (C19) reduces to: 1 wIJ 4 IJIJ SIJIJ IJ
(no sum)
(C20)
Applying (C17) to (C20) yields: w IJ
aI2 aJ2 IJ aI2 aJ2
(C21)
which is (33) in the text. D: Orthonormal base tensors for plane-straining analysis To perform inversion of fourth-order symmetric tensors in 2D space, we can follow the approach of Lebensohn et al. (1998) and Jiang (2013) by introducing a set of 3 orthonormal 2nd order tensors: 57
ACCEPTED MANUSCRIPT b1
1 1 0 1 0 1 1 1 0 2 3 , b , b 2 0 1 2 1 0 2 0 1
(D1)
This set of tensors are orthonormal, satisfying b : b . A symmetric fourth-order tensor in
PT
2D, such as the symmetric Eshelby tensor S, can be reshaped into a 3 by 3 matrix by: ( , 1..3 )
(D2)
RI
Sλη = bij bkl Sijkl
λη
This inversion ensures that 1 ik jl jkil 2
MA
1 1 Sijmn Smnkl Sijmn Smnkl
(D3)
NU
-1
1 Sijkl bij bkl S
SC
The matrix can be inverted and then reshaped back to a fourth-order tensor through:
i, j, k , l 1,2
(D4)
AC CE P
TE
D
For incompressible materials in 2D, the deviatoric stress and strain rate tensors only have 2 independent components. The corresponding 4th order stiffness tensor is singular in the complete Cauchy stress space. The inversion of such fourth-order tensors can be constrained in the deviatoric space with , 1,2 (compare with 3D case illustrated in Jiang 2013). Such inversion ensures that: 1 1 2D Sijmn Smnkl Sijmn Smnkl J ijkl
1 1 ik jl jk il ij kl 2 2
58
i, j, k , l 1,2
(D5)
ACCEPTED MANUSCRIPT REFERENCES Arbaret, L., Mancktelow, N.S., Burg, J.P., 2001. Effect of shape and orientation on rigid particle rotation and matrix deformation in simple shear flow. Journal of Structural Geology 23,
PT
113-125.
RI
Astarita, G. 1979. Objective and generally applicable criteria for flow classification, J. NonNewtonian Fluid Mechanics 6, 69–76.
SC
Barnett, D.M., 1972. The precise evaluation of derivatives of the anisotropic elastic Green’s functions. Physica Status Solidi (b), 49, 741-748.
NU
Başar, Y., and Weichert, D., 2000. Nonlinear Continuum Mechanics of Solids. Springer. Bhargava, R.D., Radhakrishna, H..C. 1964. Elliptic inclusion in orthotropic medium. Journal of
MA
the physical society of Japan, 19(3), 396-405.
Bilby, B.A., Eshelby, J.D., and Kundu, A.K., 1975, The change of shape of a viscous ellipsoidal region embedded in a slowly deforming matrix having a different viscosity:
D
Tectonophysics, v. 28, p. 265-274.
TE
Bilby, B.A., Kolbuszewski, M.L., 1977. The finite deformation of an inhomogeneity in twodimensional slow viscous incompressible flow. Proceedings of the Royal Society of
AC CE P
London, A355, 335-353.
Bird, P. 2003. An updated digital model of plate boundaries. Geochemistry Geophysics Geosystems, 4(3), 1027, doi:10.1029/2001GC000252., Brace, W.F., 1965. Relation of elastic properties of rocks to fabric. Journal of Geophysical Research, 70, 22-35.
Bretherton, F.P., 1962. The motion of rigid particles in a shear flow at low Reynolds number. Journal of Fluid Mechanics, 14, 284-304. Budiansky, B., Mangasarian, O., 1960. Plastic stress concentration at circular hole in infinite sheet subjected to equal biaxial tension. Journal of Applied Mechanics, 27, 59-64. Bunge, H.J. 1982. Texture analysis in materials science, mathematical methods. (English edition translated by Morris, P.A.), Butterworths & Co. Buryachenko, V.A., 2007. Micromechanics of heterogeneous materials. Springer. Ceriani, S., Mancktelow, N.S., Pennacchioni, G., 2003. Analogue modelling of the influence of shape and particle/matrix interface lubrication on the rotational behaviour of rigid particles in simple shear. Journal of Structural Geology 25, 2005-2021. 59
ACCEPTED MANUSCRIPT Chen, T., 1967. On an elliptic elastic inclusion in an anisotropic medium. Quarterly Journal of Mechanics and Applied Mathematics. 20, 307-313. Chen, Y., Jiang, D., Zhu, G., Xiang, B. 2014. The formation of micafish: a modeling
PT
investigation based on micromechanics. Journal of Structural Geology 68, 300-315. Cook, A.C., Vel, S.S., Gerbi, C., Johnson, S.E., 2014. Computational analysis of nonlinear creep
RI
of polyphaser aggregates: Influence of phase morphology. Journal of Geophysical
SC
Research: Solid Earth, 119, 6877-6906. doi:10.1002/2014JB011197. Davis, J.R., Titus, S.J., Horsman, E., 2013. Non-steady homogeneous deformations:
NU
Computational techniques using Lie theory, and application to ellipsoidal markers in naturally deformed rocks. Journal of Structural Geology, 56, 142-155.
MA
Dabrowski, M., Schmid, D.W. 2011. A rigid circular inclusion in an anisotropic host subject to simple shear. Journal of Structural Geology 33 (7), 1169-1177.
D
Dabrowski, M., Schmid, D.W., Podladchikov, Y.Y. 2012. A two-phase composite in simple
TE
shear: effective mechanical anisotropy development and localization potential. Journal of Geophysical Research, 117, B08406, doi:10.1029/2012JB009183.
AC CE P
DeMets, C., Traylen, S., 2000. Motion of the Rivera plate since 10 Ma relative to the Pacific and North American plates and the mantle. Tectonophyics, 318, 119–159. DeMets, C., Wilson, D.S., 1997. Relative motions of the Pacific, Rivera, North America, and Cocos plates since 0.78 Ma. Journal of Geophysical Research, 102, 2789–2807. Eshelby, J.D., 1957, The determination of the elastic field of an ellipsoidal inclusion, and related problems: Proceedings of the Royal Society of London, v. A241, p. 376-396. Eshelby, J.D., 1959. The elastic field outside an ellipsoidal inclusion. Proceedings of the Royal society of London, A 252, 561-569. Eshelby, J.D., 1961. Elastic inclusions and inhomogeneities. In: Snedden, N.I., Hill, R. (eds.) Progress in solid mechanics, Vol. 2, 89-140. Eshelby, J. D., Read, W.T., Shockley, W. 1953. Anisotropic elasticity with applications to dislocation theory. Acta Metallurgica. 1, 251-259. Exner, U. and Dabrowski, M., 2010. Monoclinic and triclinic 3D flanking structures around elliptical cracks. Journal of Structural Geology, 32(12): 2009-2021.
60
ACCEPTED MANUSCRIPT Ferguson, C.C., 1979. Rotations of elongate rod of particles in slow non-Newtonian flows. Tectonophysics 60, 247-262. Fletcher, R.C. 2004. Anisotropic viscosity of a dispersion of aligned elliptical inclusions in a
PT
homogeneous incompressible anisotropic viscous fluid. Journal of Structural Geology, 26, 1977-1987.
RI
Fletcher, R.C. 2009. Deformable, rigid, and inviscid elliptical inclusions in a homogeneous
SC
incompressible anisotropic viscous fluid. Journal of Structural Geology, 31, 382-387. Fletcher, R.C., Pollard, D.D. 1999. Can we understand structural and tectonic processes and their
NU
products without a complete mechanics? Journal of Structural Geology, 21, 1071-1088. Freeman, B. 1985. The motion of rigid ellipsoidal particles in three-dimensional slow flows:
MA
implications for geological strain analysis. Tectonophysics 113, 163-183. Freeman, B. 1987. The behavior of deformable ellipsoidal particles in three-dimensional slow flows: implications for geological strain analysis. Tectonophysics, 132, 297-309.
D
Gao, X.-L., Ma, H.M., 2010. Strain gradient solution for Eshelby’s ellipsoidal inclusion problem.
TE
Proceedoings of the Royal Society A, 466, 2425–2446. doi:10.1098/rspa.2009.0631 Gay, N.C., 1968. The motion of rigid particles embedded in a viscous fluid during pure shear
AC CE P
deformation of the fluid. Tectonophysics 5, 81-88. Ghosh, S.K., Ramberg, H., 1976. Reorientation of inclusions by combination of pure and simple shear. Tectonophysics, 34, 1-70. Goddard, J.D., Miller, C., 1967. Nonlinear effects in the rheology of dilute suspensions. Journal of Fluid Mechanics 28, 657-673. Goldsmith, H.L., Mason, S.G., 1967. The microrheology of dispersions. Rheology: theory and applications. Vol 4, 86-250. Griera A., Bons P.D., Jessell M.W., Lebensohn R., Evans L. and Gomez-Rivas E. 2011. Strain localization and porphyroclast rotation. Geology, 39, 275-278. Griera, A., Llorens, M-.G., Gomez-Rivas, E., Bons, P.D., Jessell, M.W., Evans, L.A. and Lebensohn, R. 2013. Numerical modelling of porphyroclast and porphyroblast rotation in anisotropic rocks. Tectonophysics, 587, 4–29.
61
ACCEPTED MANUSCRIPT Hashin, Z., 1991. Thermoelastic properties of particular composites with imperfect interface. Journal of Mechanics and Physics of Solids, 39, 745–762. Hill, R. 1961. Discontinuity relations in mechanics of solids. In: Snedden, N.I., Hill, R. (eds.)
PT
Progress in solid mechanics, Vol. 2, 245-276.
Hill, R. 1983. Interfacial operators in the mechanics of composite media. Journal of Mechanics
RI
and Physics of Solids, 31, 347–357
SC
Hobbs, B.E., Means, W.D., and Williams, P.F., 1976. An Outline of Structural Geology. John Wiley & Sons.
NU
Hobbs, B.E., Ord, A., 2015. Structural geology: the mechanics of deforming metamorphic rocks. Volume 1: principles. Elsevier.
MA
Huang, J.H., Furuhashi, R., Mura, T., 1993. Frictional sliding inclusions. Journal of Mechanics and Physics of Solids, 41, 247–265.
Hutchinson, J.W., 1976. Bounds and Self-Consistent Estimates for Creep of Polycrystalline
D
Materials. Proceedings of the Royal Society of London. A. Mathematical and Physical
TE
Sciences, 348(1652): 101-127.
Hwa, C., Ting, T.C.T., 1989. Two-dimensional problems of the anisotropic elastic solid with an
AC CE P
elliptic inclusion. Quarterly Journal of Mechanics and Applied Mathematics, 42, 553-572. Ildefonse, B., Mancktelow, N.S., 1993. Deformation around rigid particles: the influence of slip at the particle/matrix interface. Tectonophysics 221, 345-359. Ildefonse, B., Launeau, P., Fernandez, A., Bouchez, J.L., 1992. Effects of mechanical interactions on development of shape preferred orientations: a two-dimensional experimental approach. Journal of Structural Geology 14, 73-83. Jaswon, M.A., Bhargava, R.D., 1961. Two-dimensional elastic inclusion problems. Proc. Camb. Phil. Soc., 57, 669-680. Jeffery, G.B., 1922, The motion of ellipsoidal particles immersed in a viscous fluid: Proceedings of the Royal Society of London, v. A102, p. 161-179. Ježek, J., Melka, R., Schulmann, K., Venera, Z., 1994. The behaviour of rigid triaxial ellipsoidal particles in viscous flows -- modeling of fabric evolution in a multi-particle system. Tectonophysics 229, 165-180. Ježek, J., Schulmann, K., Segeth, K., 1996. Fabric evolution of rigid inclusions during mixed coaxial and simple shear flows. Tectonophysics 257, 203-221. 62
ACCEPTED MANUSCRIPT Jiang, D. 1999. Vorticity decomposition and its application to sectional flow characterization. Tectonophysics 301, 243-259. Jiang, D., 2007a. Numerical modeling of the motion of rigid ellipsoidal objects in slow viscous
PT
flows: A new approach. Journal of Structural Geology, 29(2): 189-200. Jiang, D., 2007b. Numerical modeling of the motion of deformable ellipsoidal objects in slow
RI
viscous flows. Journal of Structural Geology, 29(3): 435-452.
SC
Jiang, D., 2010. Flow and finite deformation of surface elements in three dimensional homogeneous progressive deformations. Tectonophysics, 487(1-4): 85-99.
NU
Jiang, D., 2012. A general approach for modeling the motion of rigid and deformable objects in ductile flows. Computers & Geosciences 38, 52-61.
MA
Jiang, D., 2013. The motion of deformable ellipsoids in power-law viscous materials: Formulation and numerical implementation of a micromechanical approach applicable to flow partitioning and heterogeneous deformation in Earth’s lithosphere. Journal of
D
Structural Geology 50, 22-34.
TE
Jiang, D. 2014. Structural geology meets micromechanics: a self-consistent model for the multiscale deformation and fabric development in Earth’s ductile lithosphere. Journal of
AC CE P
Structural Geology 68, 247-272.
Jiang, D. and Bentley, C., 2012. A micromechanical approach for simulating multiscale fabrics in large-scale high-strain zones: Theory and application. Journal of Geophysical Research 117: B12201, doi:12210.11029/12012JB009327. Jiang, D., Williams, P.F., 1998, High-strain zones: a unified model. Journal of Structural Geology, v. 20, p. 1105-1120. Jin, X., Keer, L.M., Wang, Q.J., 2011. A closed-form solution for the Eshelby tensor and the elastic field outside an elliptic cylindrical Inclusion. Journal of Applied Mechanics, Transactions ASME 78, 031009. Johnson, A.M., Fletcher, R.C., 1994. Folding of viscous layers. Columbia University Press, New York. Johnson, S.E., Lenferink, H.J., Price, N.A., Marsh, J.H., Koons, P.O., West Jr., D.P., Beane, R., 2009. Clast-based kinematic vorticity gauges: the effects of slip at matrix/clast interfaces. Journal of Structural Geology 31, 1322-1339. Jones, D.S., 1966. Generalized functions. McGraw-Hill. 63
ACCEPTED MANUSCRIPT Ju, J.W., Sun, L.Z., 1999. A novel formulation for the exterior point Eshelby’s tensor of an ellipsoidal inclusion. Transactions of the ASME, 66, 570-574. Kang, H., Milton, G.W., 2008. Solution to the polya-sazego conjecture and the weak Eshelby
PT
conjecture. Archival Rational Mechanics and Analysis, 188, 93-116. 2008 Kinoshita, N., Mura, T., 1984. Eigenstrain problems in a finite elastic body. SIAM Journal on
RI
Applied Mathematics 44, 524–535.
SC
Kocks, U.F., 1998. Kinematics and kinetics of plasticity.in Kocks, U.F., Tomé, C.N., and Wenk, H.-R. ed. Texture and anisotropy: preferred orientations in polycrystals and their effect
NU
on materials properties. Cambridge University Press. Chapter 8, 326-389. Korn, G.A., and Korn, T.M. 1968. Mathematical handbook for scientists and engineers. 2nd
MA
edition, McGraw-Hill, New York.
Lebedev, V. I., 1977. Spherical quadrature formulas exact to orders 25–29. Siberian Mathematical Journal, 18(1), 99-107.
D
Lebensohn, R. A., Dawson, P. R., Kern, H. M., Wenk, H. R., 2003. Heterogeneous deformation
TE
and texture development in halite polycrystals: Comparison of different modeling approaches and experimental data. Tectonophysics, 370, 287-311.
AC CE P
Lebensohn, R.A., and Tomé, C.N., 1993, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys: Acta Metallurgica et Materialia, 41, 2611-2624. Lebensohn, R.A., Tomé, C.N. and Maudlin, P.J., 2004. A selfconsistent formulation for the prediction of the anisotropic behavior of viscoplastic polycrystals with voids. Journal of the Mechanics and Physics of Solids, 52(2): 249-278.Lebensohn, R.A., Turner, P.A., Signorelli, J.W., Canova, G.R. and Tomé, C.N., 1998. Calculation of intergranular stresses based on a large-strain viscoplastic self-consistent polycrystal model. Modelling and Simulation in Materials Science and Engineering, 6(4): 447-465. Lekhnitskii, S.G., 1950. Theory of elasticity of an anisotropic elastic body. Gostekhizdat, Moscow (in Russian). Holden-Day, San Francisco (in English, 1963). Li, C., Jiang, D., 2011. A critique of vorticity analysis using rigid clasts. Journal of Structural Geology 33, 203-219. Li, S., Sauer, R.A., Wang, G., 2007. The Eshelby tensors in a finite spherical domain – Part I: Theoretical formulations. Journal of Applied Mechanics 74, 770–783. 64
ACCEPTED MANUSCRIPT Li, S., and Wang, G., 2008. Introduction to micromechanics and nanomechanics. World Scientific. Lin, S., Jiang, D., and Williams, P.F. 1998. Transpression (or transtension) zones of triclinic
PT
symmetry: natural example and theoretical modelling. In: Holdsworth, R.E., Strachan, R.A. & Dewey, J.F. (eds) 1998. Continental transpressional and transtensional tectonics.
RI
Geological Society, London, Special Publications, No. 135, p. 41-57.
SC
Lister, G. and Williams, P.F., 1983. The partitioning of deformation in flowing rock masses. Tectonophysics, 92(1-3): 1-33.
NU
Liu, L., 2008. Solutions to the Eshelby conjectures. Proceedings of Royal Society of London. 464, 573-594.
MA
Lubarda, V.A., Markenscoff , X., 1998. On the stress field in sliding ellipsoidal inclusions with shear eigenstrain. Journal of Applied Mechanics, 65, 858–862. Mancktelow, N. S., 2011. Deformation of an elliptical inclusion in two-dimensional
D
incompressible power-law viscous flow. Journal of Structural Geology, 33(9): 1378-1393.
TE
Mancktelow, N. S., 2013. Behaviour of an isolated rimmed elliptical inclusion in 2D slow incompressible viscous flow. Journal of Structural Geology, 46: 235-254.
AC CE P
Mancktelow, N.S., Arbaret, L., Pennacchioni, G., 2002. Experimental observations on the effect of interface slip on rotation and stabilisation of rigid particles in simple shear and a comparison with natural mylonites. Journal of Structural Geology 24, 567-585. Mandal, N., Misra, S., Samanta, S.K., 2005a. Rotation of single rigid inclusions embedded in an anisotropic matrix: a theoretical study. Journal of Structural Geology 27, 731-743. Mandal, N., Samanta, S.K., Bhattacharyya, G., Chakraborty, C., 2003. Deformation of ductile inclusions in a multiple inclusion system in pure shear. Journal of Structural Geology 25, 1359-1370. Mandal, N., Samanta, S.K., Chakraborty, C., 2005b. Numerical models of flow patterns around a rigid inclusion in a viscous matrix undergoing simple shear: implications of model parameters and boundary conditions. Journal of Structural Geology 27, 1599-1609. Markenscoff, X., Gupta, A., 2006. Collected Works of J.D. Eshelby: the mechanics of defects and inhomogeneities. Springer, Dordrecht, Netherlands.
65
ACCEPTED MANUSCRIPT Marques, F.O., Bose, S., 2004. Influence of a permanent low-friction boundary on rotation and flow in rigid inclusion/viscous matrix systems from an analogue perspective. Tectonophysics 382, 229-245.
PT
Marques, F.O., Coelho, S., 2001. Rotation of rigid elliptical cylinders in viscous simple shear flow: analogue experiments. Journal of Structural Geology 23, 609-617.
RI
Marques, F.O., Taborda, R., Bose, S., Antunes, J., 2005a. Effects of confinement on matrix flow
SC
around a rigid inclusion in viscous simple shear: insights from analogue and numerical modeling. Journal of Structural Geology, 27, 379-396.
NU
Marques, F.O., Taborda, R., Antunes, J., 2005b. 2D rotation of rigid inclusions in confined bulk simple shear flow: a numerical study. Journal of Structural Geology, 27, 2171-2180.
MA
Marques, F.O., Taborda, R., Antunes, J., 2005c. Influence of a low-viscosity layer between rigid inclusion and viscous matrix on inclusion rotation and matrix flow: a numerical study. Tectonophysics 407, 101-115.
D
Means, W.D., Hobbs, B.E., Lister, G.S., Williams, P.F., 1980. Vorticity and non-coaxiality in
TE
progressive deformation. Journal of Structural Geology 2, 371-378. Mercier, S., Molinari, A., 2009. Homogenization of elastic-viscoplastic heterogeneous materials:
1048.
AC CE P
Self-consistent and Mori-Tanaka schemes. International Journal of Plasticity, 25, 1024-
Meng, C., Heltsley, W., Pollard, D.D., 2012. Evaluation of the Eshelby solution for the ellipsoidal inclusion and heterogeneity. Computers & Geosciences, 20, 20-48. Michelitsch, T.M., Gao, H., Levin, V.M., 2003. Dynamic Eshelby tensor and potentials for ellipsoidal inclusions. Proceedings of Royal Society of London A459, 863-890. Molinari, A. 2002. Averaging models for heterogeneous viscoplastic and elastic viscoplastic materials. Transactions of the ASME, Vol 124, 62-70. Molinari, A., Canova, G.R., and Ahzi, S., 1987, A self consistent approach of the large deformation polycrystal viscoplasticity. Acta Mettalurigica, v. 35, 2983-2994. Molinari, A., Ahzi, S., Kouddane, R., 1997. On the self-consistent modeling of elastic-plastic behavior of polycrystals. Mechanics of Materials, 26, 43-62. Molinari, A., Mercier, S. 2004. Homogenization of viscoplastic materials. In: Ahzi, S., et al. (eds.), Multiscale Modeling and Characterization of Elastic-Inelastic Behavior of Engineering Materials, 113-124. 66
ACCEPTED MANUSCRIPT Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 21, 571–574. Mulchrone, K.F., 2007. An analytical solution in 2D for the motion of rigid elliptical particles
PT
with a slipping interface under a general deformation. Journal of Structural Geology, 29(6), 950-960.
RI
Mulchrone, K.F., Walsh, K., 2006. The motion of a non-rigid ellipse in a general 2D deformation.
SC
Journal of Structural Geology 28, 392-407.
Mura, T., 1987. Micromechanics of defects in solids. Martinus Nijhoff Publishers,
NU
Dordrecht/Boston/Lancaster, 587 pp.
Mura, T., 1988. Inclusion problems. Applied Mechanics Reviews 41, 15–20.
MA
Mura, T., Shodja, H.M., Hirose, Y., 1996. Inclusion problems. Applied Mechanics Reviews 49, S118–S127.
Muskhelishvili, N.I. 1953. Some basic problems of the mathematic theory of elasticity. P.
D
Noordhoff Ltd., Groningen.
TE
Nemat-Nasser, S. and Hori, M., 1999. Micromechanics: overall properties of heterogeneous materials. Second Revised Edition, Elsevier.
AC CE P
Passchier, C.W., 1987. Stable positions of rigid objects in non-coaxial flow – a study in vorticity analysis. Journal of Structural Geology 9, 679-690. Passchier, C.W., 1998. Monoclinic model shear zones. Journal of Structural Geology, 20(8): 1121-1137.
Passchier, C.W., Simpson, C., 1986. Porphyroclast systems as kinematic indicators. Journal of Structural Geology 8, 831-843. Passchier, C.W., and Trouw, R.A.J., 1996, Microtectonics. Springer. Pelletier, D., Fortin, A., Camarero, R., 1989. Are FEM solutiobns of incompressible flows really incompressible – (or how simple flows can cause headaches). International Journal for Numerical Methods in Fluids, 9, 99-112. Pennacchioni, G., Di Toro, G., Mancktelow, N.S., 2001. Strain-insensitive preferred orientation of porphyroclasts in Mont Mary mylonites. Journal of Structural Geology 23, 1281-1298. Piazolo, S., Bons, P.D., Passchier, C.W., 2002. The influence of matrix rheology and vorticity on fabric development of populations of rigid objects during plane strain deformation. Tectonophysics 351, 315-329. 67
ACCEPTED MANUSCRIPT Piazolo, S., Passchier, C.W., 2002. Experimental modeling of viscous inclusions in a circular high-strain shear rig: implications for the interpretation of shape fabrics and deformed enclaves. Journal of Geophysical Research 107, 2242.
PT
Qu, J. 1993. The effect of slightly weakened interfaces on the overall elastic properties of composite materials. Mechanics of Materials, 14, 269–281
RI
Qu, J., and Cherkaoui, M., 2006. Fundamentals of micromechanics of solids. John Wiley & Sons,
SC
Inc.
Qu, M., Lu, X., Jiang, D. in review. An optimal method for calculating Eshelby tensors and a
NU
MATLAB package for simulation the motion of viscous inclusions. Submitted to Computers & Geosciences.
MA
Ramsay, J.G. 1967. Folding and fracturing of rocks. McGraw-Hill, New York. Ramsay, J.G., Graham, R.H. 1970. Strain variation in shear belts. Canadian Journal of Earth Sciences, 7, 786-813.
D
Rand, O., Rovenski, V., 2005. Analytical methods in anisotropic elasticity with symbolic
TE
computational tools. Birkäuser, Boston. Rudnicki, J.W., 1977. The inception of faulting in a rock mass with weakened zone. Journal of
AC CE P
Geophysical Research, 82, 844-854. Schmid, D.W., Podladchikov, Y.Y., 2003. Analytical solutions for deformable elliptical inclusions in general shear. Geophysical Journal International. 155, 269-288. Schmid, D.W., Podladchikov, Y.Y., 2004. Are isolated stable rigid clasts in shear zones equivalent to voids? Tectonophysics 384, 233-242. Schmid, D.W., Podladchikov, Y.Y., 2005. Mantled porphyroclast gauges. Journal of Structural Geology 27, 571-585. Spence, D.A., Wilmott, P., 1988. Rotation and deformation of a viscous inclusion in Stokes-flow. Proceedings of the Royal Society of London, Series A 418, 383-403. Spencer, A.J.M., 1980. Continuum mechanics. Longman, London. Stroh, A.N., 1958. Dislocations and cracks in anisotropic elasticity. Philosophic Magazine, 3, 625-646. Stroh, A.N., 1962. Steady state problems in anisotropic elasticity. Journal of Mathematical Physics. 41, 77-103. Synge, J.L. 1957. The hypercircle in mathematical physics. Cambridge University Press. 68
ACCEPTED MANUSCRIPT ten Brink, C.E., Passchier, C.W., 1995. Modelling of mantled porphyroclasts using nonNewtonian rock analogue materials. Journal of Structural Geology 17, 131-146. ten Grotenhuis, S.M., Passchier, C.W., Bons, P.D., 2002. The influence of strain localisation on
PT
the rotation behaviour of rigid objects in experimental shear zones. Journal of Structural Geology 24, 485-499.
RI
Ting, .C.T., 1996. Anisotropic elasticity: theory and applications. Oxford University Press.
SC
Ting, T. C. T. & Lee, V.-G. 1997. The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids. Quaterly Journal of Mechanics and Applied
NU
Mathematics. 50, 407–426.
Treagus, S.H., Treagus, J.E., 2002. Studies of strain and rheology of conglomerates. Journal of
MA
Structural Geology 24, 1541-1567.
Tullis, J. 2002. Deformation of granitic rocks: experimental studies and natural examples. Review of Mineralogy and Geochemistry. 51, 51-95.
D
Walke, K.P. 1993. Fourier integral representation of the Green function for an anisotropic elastic
TE
half-space. Proccdings of the Royal Society of London A. 443, 367-389. Walpole, L.J. 1978. A coated inclusion in an elastic medium. Math Proc Camb, 83:495–506
AC CE P
Weinberger, C., Cai, W., Barnett, D. 2005. Lecture notes – elasticity of microscopic structures. Stanford University (http://micro.stanford.edu/~caiwei/me340b/course_materials.html). White, S.H., Burrows, S.E., Carreras, J., Shaw, N.D., Humphreys, F.J., 1980. On mylonites in ductile shear zones. Journal of Structural Geology, 2, 175-187. Willis, D.G., 1977. A kinematic model of preferred orientation. Geological Society of America Bulletin 88, 883-894.
Willis, J.R., 1964. Anisotropic elastic inclusion problems. Quaterly Journal of Mechanics and Applied Mathematics 17, 157-174. doi:10.1093/qjmam/17.2.157. Xiang, B. and Jiang D., 2013. Small-scale ductile shear zones as transposed rheologically weak domains: A numerical modeling investigation and practical application. Journal of Structural Geology 54: 184-198. Xu, B.-X., Muller, R., Wang, M.-Z. 2011. The Eshelby property of sliding inclusions. Archive of Applied Mechanics., 81, 19-35. Yang, H.C., Chou, Y.T. 1976. Generalized plane problems if elastic inclusions in anisotropic solids. Journal of Applied Mechanics, 43, 424-430. 69
ACCEPTED MANUSCRIPT Yin, H.M., Lee, P.-H., Liu, Y. J., 2014. Equivalent inclusion method for the Stokes flow of drops moving in a viscous fluid. Journal of Applied Mechanics 81, 071010. doi:10.1115/1.4027312.
PT
Zheng, Q.-S., Zhao, Z.-H., Du, D.-X., 2006. Irreducible structure, symmetry and average of Eshelby’s tensor fields in isotropic elasticity. Journal of the Mechanics and Physics of
RI
Solids, 54, 368–383.
SC
Zhou, K., Hoh, H.J., Wang, X., Keer, L.M., Pang, J.H.L., Song, B., 2013. A review of recent works on inclusions. Mechanics of Materials, 60, 144-158.
NU
Zou, W.N., He, Q.C., Zheng, Q.S., 2012. Inclusions in a finite elastic body. International Journal
AC CE P
TE
D
MA
of Solids and Structures 49, 1627–1636.
70
ACCEPTED MANUSCRIPT Figure Captions Figure 1: Stereographic projection showing how, for any given
x , unit vector z in the plane
x (i.e., the x z 0 plane) can be expressed in terms of two mutually orthogonal unit vectors ( α and β ) on the x z 0 plane by z cos α sin β following Synge (1957) and
PT
normal to
SC
RI
cos sin Barnett (1972). When x is expressed in terms of spherical angles and by x sin sin , cos
NU
sin cos cos then α cos , and β sin cos . Green functions and their derivatives ( Gij , Gij ,l , and 0 sin
H i ) for anisotropic viscous materials are expressed as contour integrals along a unit circle on
MA
the x z 0 plane. See text for more details.
D
E Figure 2: A general anisotropic elliptic or ellipsoidal inclusion with viscous stiffnesses Cijkl
TE
embedded in an infinite matrix with viscous stiffnesses Cijkl subjected to far-field uniform
AC CE P
deformation. The mechanical fields inside the inclusion are uniform and denoted by lowercase symbols. The mechanical fields in the matrix (called exterior fields) vary with space (as functions of the position vector x ). The uniform loading or motion at infinity are denoted by uppercase symbols. See text for more details.
Figure 3: Schematic diagrams showing the geometries of the a1 a2 section of an inclusion that goes through a circular state at time t . The ellipse at time t t is reciprocal to the ellipse at time t t . Therefore, the two semiaxes swap identities prior and after the circular state but do not otherwise change directions. See text for more details.
Figure 4: An inclusion embedded in a matrix with planar anisotropy. The matrix rheology is defined by two distinct viscosities, n and s . m n defines the strength of anisotropy. The s far-field general shear flow (pure shearing and simple shearing ) is defined in the
coordinate system x1 x2 parallel to the anisotropy. The inclusion’s principal axes define another coordinate system ( x '1 x '2 ) at angle relative to the x1 x2 system.
71
ACCEPTED MANUSCRIPT Figure 5: Comparison of the interior fields of an inclusion ( R 5 , r 10 ) embedded in a matrix with m 25 ) calculated numerically from the 2D formulation of this paper with results calculated with Fletcher’s analytic solutions. (a) Strain rates differences (between inclusion and
NU
SC
RI
PT
far-field matrix) 11 and 12 . (b) Vorticity difference w11 . (c) Pressure differences p normalized against far-field stress invariant. The solid color curves are numerical results using the 2D formulation of this paper implemented in MATHCAD and the dashed black curves are calculated from analytical expressions (Eq.57). The far-field flow is plane-strain general shear with pure shearing 0.5 and simple shearing 1 both along the plane of anisotropy (Fig.4). Numerical results and those calculated directly using Eqs.57 (dashed lines). Numerical results agree with analytic results perfectly, suggesting excellent accuracy of the numerical approach.
AC CE P
TE
D
MA
Figure 6: Pressure perturbations in and around a strong isotropic inclusion ( R 3 , r 100 ). The pressure is normalized against the matrix deviatoric stress invariant. The far-field flow is pure shearing with horizontal stretching and vertical shortening. The inclusion long axis is horizontal, parallel to the maximum pure shearing axis. (a): Pressure perturbation in isotropic matrix ( m 1 ). (b)-(e): Matrix is anisotropic with m 10 . The anisotropy plane is horizontal in (b), 30°, 60°, and 90° clockwise with respect to the horizontal anisotropy in (c), (d), and (e) respectively. See text for more details.
Figure 7: Angular velocities (normalized against shearing) of rigid inclusions with a circular shape in the x1 x2 plane as functions of the matrix m for different inclusion a3 relative lengths. The far-field plane-straining flow is with pure shearing 0.5 and simple shearing 1 both along the plane of anisotropy (Fig.3). A circular cylinder ( 1:1: ) rotates at a constant angular speed of half the shearing, independent of m . For finite a3 inclusions, the angular velocities increase with m . As the relative length of a3 is reduced, the inclusion’s angular velocity increases. A disk with vanishing a3 has the highest angular velocity (the 1:1:0.001 case).
Figure 8: (a) Angular velocity variations of a rigid triaxial ellipsoid ( 6 : 1: 5 , with the a3 axis parallel to the x3 axis) as functions of its sectional orientation in the x1 x2 plane for different m . The matrix far-field flow is the same as for Fig.6. The horizontal axis is the angle between the inclusion’s sectional long axis and the anisotropy ( x1 axis in Fig.3). A greater m corresponds to a higher angular velocity at any given orientation, except near the stable orientation of the object. (b): Angular velocity variations of ellipsoids with the same sectional aspect ratio of 6:1 but 72
ACCEPTED MANUSCRIPT different relative a3 lengths in a matrix with fixed m =10. A relatively shorter a3 generally
PT
corresponds to a higher angular velocity.
Figure 9: Angular velocities of rigid ellipsoids with the same sectional aspect ratio (2:1) on the x1 x2 plane but varying a3 lengths, in a matrix with m =10 undergoing simple shearing parallel to
SC
RI
the anisotropy plane. The angular velocity of an inclusion depends on its sectional orientation ( ) as well as its relative a3 length. An elliptic cylinder rotates like in an isotropic matrix. As long as a3 is finite, its angular velocity increases with decreasing a3 . A disk with a vanishing a3
NU
rotates at a faster and finite rate.
AC CE P
TE
D
MA
Figure 10: Schematic diagrams showing that in the case of an inclusion with a slippery interface contact with the matrix the shear spin effect must be considered in the angular velocity equation. The system can be regarded as a matrix with a void and an inclusion in the void. (a) A state of the system of a matrix containing a slippery inclusion. A material point (p) on the major radius at the inclusion side of the interface coincides with the material point (p’) on the matrix side of the interface. The points p and p’ are at the same location but are drawn slightly apart across the interface for easy visualization. (b) The state of the system after an increment of deformation. Because of free slip along the interface, p’ has moved along the hoop with respect to p. That is the void’s principal axes have angular velocity swiping through material.
73
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 1
74
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 2
75
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 3
76
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 4
77
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 5
78
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 6
79
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 7
80
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 8
81
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 9
82
ACCEPTED MANUSCRIPT
AC CE P
TE
D
MA
NU
SC
RI
PT
Fig. 10
83