# Contrast Expansion Method for Elastic Incompressible Fibrous Composites

Authors: Piotr Drygaś; Vladimir Mityushev
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

Contrast parameter expansion of the elastic fields for 2D composites is developed by Schwarz’s method and by the method of functional equations for the case of circular inclusions. A computationally efficient algorithm is described and implemented in symbolic form to compute the local fields in 2D elastic composites and the effective shear modulus for macroscopically isotropic composites. The obtained new analytical formula contains high-order terms in the contrast parameter and explicitly demonstrates dependence on the location of inclusions. As a numerical example, the hexagonal array is considered.


---

## Body

## 1. Introduction

The general potential theory of mathematical physics yields methods of integral equation to numerically solve various boundary value problems. Integral equations for plane elastic problems were constructed by Muskhelishvili [1], first, extended to doubly periodic problems in [2], and developed in [3–5] and papers cited therein. The obtained results were applied to computations of the effective properties of the elastic media. Integral equations are efficient for the numerical investigation of nondilute composites when interactions of inclusions have to be taken into account. Another type of integral equation based on the generalized alternating method of Schwarz was proposed by Mikhlin [6] and developed in [7].Consider a problem on the planeR2 isomorphic to the complex plane C. Let D≔C∪∞∖⋃k=1nDk∪∂Dk where the boundary ∂Dk of Dk is a simple closed Lyapunov curve oriented in clockwise sense. Schwarz’s method for a multiply-connected domain D=C^∖⋃k=1n(Dk∪∂Dk) is based on the separate solutions of the simple boundary value problems for simply connected domains C^∖(Dk∪∂Dk)(k=1,2,…,n). It was demonstrated in [7] that Schwarz’s method can be realized as the iterative schemes constructed on contrast and on concentration parameters considered as small perturbation parameters and precisely described below. Convergence was proved for Laplace’s equation for the contrast expansion.The main advantage of Schwarz’s method consists in analytical solution to the problems discussed when the physical parameters are presented in symbolic form in the final exact or approximate formulae. Such formulae were recently obtained in [8] for biharmonic functions which describes elastic materials for a circular multiply-connected domain. The results of [8] are based on the concentration expansion.In the present paper, we apply Schwarz’s method based on the contrast parameter expansion for a circular multiply-connected domain. Schwarz’s method is used in the form of the functional equations method [7–9]. Elastic isotropic materials are described through two independent moduli. In order to make the presentation clear, we simplify the problem by consideration of incompressible materials when Poisson’s ratio ν is equal to 0.5. Then, we introduce only one contrast parameter:(1)ρ=μ1/μ-1μ1/μ+1,where μ1 and μ denote the shear modulus of inclusions Dk and matrix D, respectively.In the present paper, we deduce a computationally efficient algorithm implemented in symbolic form to compute the local fields in 2D elastic composites and the effective shear modulusμe for macroscopically isotropic composites. The obtained new analytical formula is valid up to O(ρ3) and explicitly demonstrates dependence on the location of inclusions. Such an approach has advantages over pure numerical methods when dependencies of the effective constants on the mechanical properties of constitutes and on geometrical structure are required.The numerical examples from Section6 give sufficiently accurate values of μe for all admissible ρ, that is, for |ρ|≤1 and for f not exceeding 0.4.


## 2. Contrast Expansion Method

Consider a finite numbern of inclusions as mutually disjoint simply connected domains Dk(k=1,2,…,n) in the complex plane C≡R2. It is worth noting that the number n is given in a symbolic form with an implicit purpose to pass to the limit n→∞ later.The components of the stress tensor can be determined by the Kolosov-Muskhelishvili formulae [1]: (2)σxx+σyy=4Reφk′z,z∈Dk,4Reφ0′z,z∈D,σxx-σyy+2iσxy=-2zφk′′z¯+ψk′z¯,z∈Dk,-2zφ0′′z¯+ψ0′z¯,z∈D.The strain tensor components of incompressible materials ϵxx=-ϵyy,  ϵxy are determined by the formulae [1]: (3)ϵxx+iϵxy=-12μ1zφk′′z¯+ψk′z¯,z∈Dk,-12μzφ0′′z¯+ψ0′z¯,z∈D.Let the stress tensor be applied at infinity:(4)σ∞=σxx∞σxy∞σyx∞σyy∞.The uniform shear stress will be considered in the paper when (5)σxx∞=σyy∞=0,σxy∞=σyx∞=1.Then, (6)ψ0z=iz+ψz,where the functions φ0(z) and ψ(z) are analytical in D and bounded at infinity. The functions φk(z) and ψk(z) are analytical in Dk and twice differentiable in the closures of the considered domains.It follows from (1) that (7)μ1μ=1+ρ1-ρ.The perfect bonding at the matrix-inclusion interface can be expressed by two equations [1]: (8)φkt+tφk′t¯+ψkt¯=φ0t+tφ0′t¯+ψ0t¯,t∈∂Dk,k=1,2,…,n,(9)φkt-tφk′t¯-ψkt¯=1+ρ1-ρφ0t-tφ0′t¯-ψ0t¯,t∈∂Dk,k=1,2,…,n.The problem (8)-(9) is the classic boundary value problem of the plane elasticity having the unique solution up to additive constants corresponding to rigid shifts of medium. It follows from [1] that solutions of (8)-(9) analytically depend on ρ for sufficiently small |ρ|. One can see also Chapter 2, Section 3 in [3], where the problem (8)-(9) is reduced to the Fredholm integral equation shortly written as (10)Φ=μ1μ-1KΦ+g,where K is a compact integral operator in the space of the Hölder continuous functions.We are looking for the complex potentials in the contrast expansion form. For instance,φk(z) for sufficiently small |ρ| has the following form: (11)φkz=∑l=0∞φklzρl.Then, the problem (8)-(9) is reduced to the following cascade of boundary value problems. Equation (8) becomes(12)φklt+tφkl′t¯+ψklt¯=φ0lt+tφ0l′t¯+ψ0lt¯,t∈∂Dk,k=1,2,…,n,l=0,1,….Equation (9) yields the cascade (13)φk0t-tφk0′t¯-ψk0t¯=φ00t-tφ00′t¯-ψ00t¯,φklt-tφkl′t¯-ψklt¯=Fl-1t,t∈∂Dk,k=1,2,…,n,l=1,2,…,where (14)Fl-1t=2φ0l-1t-tφ0l-1′t¯-ψ0l-1t¯.Addition and subtraction of the boundary condition (12) for l=0 and the first condition (13) yield (15)φk0t=φ00t,ψk0t=ψ00t,t∈∂Dk,k=1,2,…,n.Application of principle of analytic continuation to (15) implies analytic continuation of the functions φk(0)(z) and φ0(0)(z) into all the domains Dk. Then, (6) gives the exact formulae for the zero approximation: (16)φ00z=0,ψ00z=izand F0(t)=2it¯.Ultimately, we arrive at the following cascade. First, we solven boundary value problems separately for every domain Dk: (17)φk1t-tφk1′t¯-ψk1t¯=2it¯,t∈∂Dk.Next, we solve the problem following from (13): (18)φ01t+tφ01′t¯+ψ01t¯=φk1t+tφk1′t¯+ψk1t¯,t∈∂Dk,k=1,2,…,n.This is a boundary value problem for the multiply-connected domain D on the functions φ0(1)(z), ψ0(1)(z) analytic in D and twice continuously differentiable in the closure of D including infinity. The described above step is the first step of the iterative scheme when we pass from φ0(0)(z) and ψ0(0)(z) to φ0(1)(z) and ψ0(1)(z). The step (l+1) consists in the solution to the problems for every domain Dk: (19)φkl+1t-tφkl+1′t¯-ψkl+1t¯=Flt,t∈∂Dk,k=1,2,…,n,and, further, the problem for the domain D(20)φ0l+1t+tφ0l+1′t¯+ψ0l+1t¯=φkl+1t+tφkl+1′t¯+ψkl+1t¯,t∈∂Dk,k=1,2,…,n.Therefore, the conjugation problem (8)-(9) is reduced to the sequence of the problems separately for the domains Dk (k=1,2,…,n) and D. The described iterative scheme is computationally effective if the inclusions Dk have simple shape. The next section is devoted to its explicit realization for circular inclusions.


## 3. Method of Functional Equations

Consider the circular inclusionsDk≔z∈C:z-ak<r(k=1,2,…,n) and D≔C∪∞∖⋃k=1nDk∪∂Dk, where ∂Dk≔t∈C:t-ak=r. Introduce the new unknown functions (21)Φkz=r2z-ak+ak¯φk′z+ψkz,z-ak≤r,analytic in Dk except at the point ak, where its principal part has the form r2z-ak-1φk′(ak).Letzk∗=r2z-ak¯-1+ak denote the inversion with respect to the circle ∂Dk. If a function f(z) is analytic in z-ak<r, then f(zk∗)¯ is analytic in z-ak>r. The problem (8), (9) was reduced in [9] (see equations (5.6.11) and (5.6.16) in Chapter 5), [10] to the system of functional equations: (22)φkz=ρ∑m≠kΦmzm∗¯-z-amφm′am¯-ρφk′ak¯z-ak+p0,(23)Φkz=ρ∑m≠kφmzm∗¯+ρ∑m≠kzk∗¯-zm∗¯Φmzm∗¯′-φm′am¯+1+ρiz-ρ∑k=1nr2φk′ak¯z-ak+q0,z-ak≤r,k=1,2,…,n,where p0 and q0 are constants. The unknown functions φk(z) and Φk(z) (k=1,2,…,n) are related by 2n equations (22)-(23). One can see that the functional equations do not contain integral operators but contain compositions of φk(z) and Φk(z) with inversions.Following [11], we first introduce the Hardy-Sobolev space H(2,2)(Dk) separately for each k=1,…,N as the space of functions analytic in Dk satisfying the conditions: (24)sup0<r<rk∫02πfqreiθ+ak2dθ<∞,q=0,1,2,where f(q) denotes the derivative dqf/dzq. The norm is introduced as follows:(25)fH2,2Dk2=fH2Dk2+f′H2Dk2+f′′H2Dk2,where the classic Hardy norm (26)fH2Dk2=sup0<r<rk∫02πfreiθ+ak2dθis used.Using the above designations (24)–(26) we introduce the space H(2,2)(⋃k=1NDk) shortly denoted by H(2,2) of functions f analytic in ⋃j=1NDk endowed with the norm: (27)fH2,22≔∑k=1NfkH2,2Dk2,where f(z)=fk(z) for z∈Dk(k=1,…,N).The functional-differential equations (22)-(23) include the meromorphic functions Φk(z) not belonging to H(2,2). They were written as(28)ϕ=Aϕ+ηon the vector-function ϕ(z)=(φk(z),ψk(z))T introduced in all Dk by substitution (21) in the space H(2,2)×H(2,2). The operator A and the given vector-function η are determined by (22)-(23). Equation (28) was explicitly written in [11] as the system of functional-differential equations (29)φkz=ρ∑m≠kzφm′zm∗¯+ψmzm∗¯-z-amφm′am¯-ρφk′ak¯z-ak+p0,z-ak≤rk,k=1,2,…,N,(30)ψkz=ρ∑m≠kφmzm∗¯-ρ∑m≠kzm∗¯φm′zm∗¯+zφm′zm∗¯′+ψmzm∗¯′-φm′am¯-ρ∑m≠krm2z-amψm′am¯+1+ρiz+ρak¯φk′ak¯+q0,z-ak≤rk,k=1,2,…,N.It was proved in [11] that the operator A is compact in H(2,2)×H(2,2). One can see that the contrast parameter ρ plays the role of the spectral parameter; hence, ϕ can be written in the form of power series in ρ. This implies that the method of successive approximations applied to (29) converges in H(2,2)×H(2,2) for sufficiently small ρ. Let φk∗(z), ψk∗(z) (k=1,2,…,n) be a solution of (29). This unique solution belongs to H(2,2)×H(2,2). The given vector-function η is twice differentiable in the closures of Dk. Hence, the pumping principle [9, page 22] can be applied as follows. The shifts in composition operators of the right part of (29) are the shift strictly into domains. Hence, if we substitute φk∗(z), ψk∗(z) into (29), we obtain that φk∗(z), ψk∗(z) are twice differentiable in the closures of Dk. The equivalent method of successive approximations can be applied to (22)-(23) more conveniently in computations.The functionsφ(z) and ψ(z) are determined through φk(z) and Φk(z) up to additive constants. For instance, (31)φz=ρ1+ρ∑m=1nΦmzm∗¯-z-amφm′am¯,z∈D.


## 4. Method of Successive Approximations

Application of successive approximations to functional equations is equivalent to Schwarz’s method described at the end of Section2. It follows, for instance, from the uniqueness of the analytic expansion in ρ near zero. Moreover, each iteration for the functional equations corresponds to an iteration step in Schwarz’s method since the coefficients in the series in ρ are also uniquely determined. It can be also established directly form formulae written below.Using the series (11) for φk(z) and analogous series for other functions and applying successive approximations to the functional equations we obtain the following iteration scheme. The zeroth approximation is (32)φk0z=0,Φk0z=iz.The next approximations for s=1,2,… are (33)φksz=∑m≠kΦms-1zm∗¯-z-amφms-1′am¯-z-akφks-1′ak¯,(34)Φksz=∑m≠kφms-1zm∗¯-∑m=1nr2z-amφms-1′am¯+∑m≠kzk∗¯-zm∗¯Φms-1zm∗¯′-φms-1′am¯+δ0s+δ1siz,where δj,k denotes the Kronecker symbol.Introduce the functions(35)Fsz=∑m=1nΦmszm∗¯-z-amφms′am¯.Then, (36)∑m=1nΦmzm∗¯-z-amφm′am¯=∑s=0∞Fszρs.It follows from (31) that φ(z) can be written in the following form: (37)φz=∑s=0∞∑j=0s-1-1s-j-1Fjzρs,where the following expansion is used: (38)ρ1+ρ=∑l=0∞-1lρl+1.The sth approximation for φ(z) becomes (39)φsz=∑j=0s-1-1s-j-1Fjz,s=1,2,…and φ(0)(z)=0.


## 5. Shear Modulus

### 5.1. General Iterative Scheme

Introduce the average value over a sufficiently large rectangleQn containing all the inclusions Dk: (40)wn=1Qn∬Qnwdx1dx2.The averaged shear modulus μe(n) of the considered finite composite is introduced as the ratio (41)μen=σ12n2ϵ12n.It is related to the effective shear modulus μe for macroscopically isotropic composites by the limit μe=limn→∞⁡μe(n). It was demonstrated in [8] that (41) can be transformed into (42)μeμ=1+ReA1-ReA,where(43)A=limn→∞12Qn∑k=1n∫∂Dkφtdt.Here, φ(z) is a solution to the problem with n inclusions.In order to determineA, we first calculate the integral ∫∂Dkφ(z)dz. Equation (37) implies that (44)∫∂Dkφtdt=∑s=0∞∑j=0s-1-1s-j-1∫∂DkFjtdtρs.Application of (35) and Cauchy’s integral theorem yield (45)∫∂DkFszdz=∑m=1n∫∂DkΦmszm∗¯dz.The presented iterative scheme can be easily realized numerically. But we are interested in analytical formulae which can be obtained by symbolic computations. In the next sections, we perform symbolic computations to determine φ(z) and the effective shear modulus in the third-order approximation.


### 5.2. Second-Order Approximation

We are looking forφ(z) in the form (46)φz=φ1zρ+φ2zρ2+φ3zρ3+Oρ4.It follows from (37) that for the third-order approximation we need the integrals (45) for j=0,1,2. Using the second equation (32), we have (47)∫∂DkF0zdz=∑m=1n∫∂DkΦm0zm∗¯dz=-i∑m=1n∫∂Dkzm∗¯dz.The functions zm∗¯=r2/z-am+am¯ are analytic in Dk for m≠k. Therefore, the integrals with m≠k vanish in (47) by Cauchy’s integral theorem and (47) can be calculated by residues (48)∫∂DkF0zdz=-i∫∂Dkzk∗¯dz=2πr2.In order to use (37) for s=1 we find from (34) (49)Φk1z=iz+∑m≠kzk∗¯-ak¯-zm∗¯-am¯+ak¯-am¯r2z-am2.Using the expansion for m≠k, (50)z-am-l=∑j=0∞-1jl+j-1jak-am-l+jz-akj,we have(51)zk∗¯-ak¯r2z-am2=r2z-akr2z-am2=r4ak-am-2z-ak-1-r4∑j=0∞-1jj+2ak-am-j+3z-akj,(52)zm∗¯-am¯r2z-am2=r42∑j=0∞-1jj+1j+2ak-am-j+3z-akj.Along similar lines (53)ak¯-am¯r2z-am2=r2∑j=0∞-1jj+1jak¯-am¯ak-am-j+2z-akj.Subtracting (51) and (52) we obtain (54)zk∗¯-ak¯r2z-am2-zm∗¯-am¯r2z-am2=r4ak-am-2z-ak-1-r42∑j=0∞-1jj+2j+3ak-am-j+3z-akj.Using the above formulae we obtain from (49) the exact formula:(55)Φm1z=iz+∑m1≠mr4am-am1-2z-am-1-r42∑j=0∞-1jj+2j+3am-am1-j+3z-amj+r2∑j=0∞-1jj+1am¯-am1¯am-am1-j+2z-amj.We are now ready to calculate the integral: (56)∫∂DkF1zdz=∫∂Dk∑m=1nΦm1zm∗¯dz=-i∫∂Dk∑m=1nzm∗¯+∑m1≠mr4am-am1¯-2zm∗-am¯-1-r42∑j=0∞-1jj+2j+3am-am1¯-j+3zm∗-am¯j+r2∑j=0∞-1jj+1am-am1am-am1¯-j+2zm∗-am¯jdz.Using Cauchy’s integral theorem and residues we obtain(57)∫∂Dk∑m=1nzm∗-am¯-1dz=∫∂Dk∑m=1nz-amr2dz=0,∫∂Dk∑m=1nzm∗-am¯jdz=∫∂Dkr2z-akjdz=2πir2forj=10otherwise.This yields (58)∫∂DkF1zdz=2πr2+12πr6∑m≠kak-am¯-4-4πr4∑m≠kak-amak-am¯-3.


### 5.3. Third-Order Approximation

The third-order approximation requires advanced and long computations presented below. In order to calculateφ(3), we need the following function:(59)Φm2z=-i∑m1≠m∑m2≠m1zm1∗m2∗¯-∑m1=1nr2z-am1∑m2≠m1r2am1¯-am2¯2+∑m1≠mzm∗¯-am¯-zm1∗¯-am1¯+am¯-am1¯r2z-am12-∑m2≠m1r2r2/z-am1+am1¯-am2¯2+2r4r2/z-am1+am1¯-am2¯3z-am1-3r6r2/z-am1+am1¯-am2¯4z-am12+2r4am1-am2r2/z-am1+am1¯-am2¯3z-am12-∑m2≠m1r2am1¯-am2¯2.It is obtained from (34) for s=2 by substitution Φk(1) given by (49) and φk(1) calculated by (33) with s=1.Below, we describe general recurrent formulae for an arbitrarysth approximation and explicitly write the third-order terms when s=3. The following double series in ρ and (z-ak) are used: (60)φkz=∑s=0∞∑j=0∞φk,jsz-akjρs,(61)Φkz=r2z-akφm,1+∑j=0∞γm,jz-akj.The coefficients in the power series (61) are presented as series in ρ: (62)φm,l=∑s=0∞φm,lsρs,γm,l=∑s=0∞γm,lsρs.Substitution of (60)–(62) into (22) and selection of coefficients in the same powers of ρ and (z-ak) yield(63)φk,js=-δj,1φk,1s-1¯+-1j∑m≠k∑l=1∞γm,ls-1¯r2ll+j-1jak-am-l+jwith the zeroth approximation φk,j(0)=0. Along similar lines, (23) yields (64)γk,js=iδj,1δ0,s+δ1,s+∑m≠k--1jr2φm,1s-1¯ak-am-j+1+-1j∑l=1∞φm,ls-1¯r2ll+j-1jak-am-l+j+-1j∑l=1∞lγm,ls-1¯r2l+1l+j+2j+1ak-am-l+j+2--1jak-am¯∑l=1∞lγm,ls-1¯r2ll+jjak-am-l+j+1with the zeroth approximation γk,j(0)=iδ1,j.The expanded form ofφ(z) follows from (31):(65)φz=ρ1+ρ∑k=1n∑s=0∞∑j=0∞γk,js¯r2jz-ak-jρs.Equation (65) is in accordance with (37) where(66)Fjz=∑m=1n∑l=0∞γm,lj¯r2z-am-l.Using (66) we calculate the integral (67)∫∂DkFjtdt=∫∂Dk∑m=1n∑l=0∞γm,lj¯r2t-am-ldt=∫∂Dk∑l=0∞γk,lj¯r2t-ak-ldt=2πiγk,1j¯r2.The coefficientsγm,l(j) for j=0,1,2 can be explicitly written by the iterations (64): (68)γk,10=i,γk,11=i-ir22∑m≠kak-am¯ak-am3-6r2∑m≠k1ak-am4,(69)γk,12=i-r4∑k1≠k∑k2≠k11ak-ak121ak1¯-ak2¯2+6r4∑k1≠k1ak-ak14-2r2∑k1≠kak¯-ak1¯ak-ak13+∑l=1∞-1l-r2+2ll∑k1≠k∑k2≠k11ak-ak1l+11ak1¯-ak2¯l+1-14r2l+3ll+2l+32∑k1≠k1ak¯-ak1¯l+3∑k1≠k1ak-ak1l+3+12r2l+2ll+1l+2l+3∑k1≠kak-ak1ak¯-ak1¯l+2∑k1≠k1ak-ak1l+3+12r2l+2ll+1l+2l+3∑k1≠k1ak¯-ak1¯l+3∑k1≠kak¯-ak1¯ak-ak1l+2-r2l+1ll+12∑k1≠kak-ak1ak¯-ak1¯l+2∑k1≠kak¯-ak1¯ak-ak1l+2.The limit (43) has to be found. Using (44) and (67) we express A through γk,1(j): (70)A=∑s=0∞∑j=0s-1-1s-j-1πir2limn→∞1Qn∑k=1nγk,1j¯ρs.Taking into account terms up to O(ρ4) we rest γk,1(j) for j=0,1,2; hence, (71)A=πir2limn→∞1Qn∑k=1nγk,10¯ρ+πir2limn→∞1Qn∑k=1nγk,11¯-γk,10¯ρ2+πir2limn→∞1Qn∑k=1nγk,12¯-γk,11¯+γk,10¯ρ3+Oρ4.The above limits are calculated by use of the formalism developed in [8, 12]. It was supposed that N equal nonoverlapping disks belong to a parallelogram (fundamental) cell periodically extended to the complex plane by two linearly independent translation vectors. As an example, we shortly present the transformation of the term of γk,1(1) given by (68) multiplied by πir2 and complexly conjugated in accordance with (71): (72)Tn=12πr6∑m≠k1ak-am¯4.We are interested in the limit associated with the Eisenstein summation [8, 12]:(73)T=limn→∞Tn2Qn.The normalized e-sums were introduced in [8, 12] by means of the Eisenstein summation: (74)ep=limn→∞1N1+p/2Qnp/2∑m≠k1ak′-am′p,(75)epp=limn→∞1N1+pQnp∑m≠k∑k≠k11ak′-am′pak′-ak1′¯p,(76)ep1=limn→∞1N1+p-1/2Qnp-1/2∑m≠kak′-am′¯ak′-am′p,where ak′=Qn-1/2ak. This means that the rectangle Qn is normalized by the C-linear transformation z→z′=Qn-1/2z to Qn′ having the unit area. Therefore, T=6f3/π2e4¯. Other terms are transformed by the same method. As a result we arrive at the following formula:(77)A=fρ+2f3f2π2e4¯-fπe31¯ρ2+f-f2π2e2,2+∑l=1∞-1l-fπl+1lel+1,l+1-14fπl+3ll+2l+32el+32+12fπl+2ll+1l+2l+3el+21¯el+3+12fπl+2ll+1l+2l+3el+3¯el+21-fπl+1ll+12el+212ρ3+Oρ4,where (78)f=limn→∞Nπr2Qndenotes the concentration. The absolute convergence of series in R3 follows from the standard root test (Cauchy’s criterion).The effective shear modulus is calculated by (42):(79)μeμ=1+fρ+R2ρ2+R3ρ31-fρ-R2ρ2-R3ρ3,where(80)R2=23e4¯π2f3-e31¯πf2,(81)R3=ff2π2e2,2+∑l=1∞-1l-fπl+1lel+1,l+1-14fπl+3ll+2l+32el+32+12fπl+2ll+1l+2l+3el+21¯el+3+12fπl+2ll+1l+2l+3el+3¯el+21-fπl+1ll+12el+212.Expansion of (79) in ρ yields (82)μeμ=1+2fρ+2f2+R2ρ2+2f3+2fR2+R3ρ3+Oρ4.


## 5.1. General Iterative Scheme

Introduce the average value over a sufficiently large rectangleQn containing all the inclusions Dk: (40)wn=1Qn∬Qnwdx1dx2.The averaged shear modulus μe(n) of the considered finite composite is introduced as the ratio (41)μen=σ12n2ϵ12n.It is related to the effective shear modulus μe for macroscopically isotropic composites by the limit μe=limn→∞⁡μe(n). It was demonstrated in [8] that (41) can be transformed into (42)μeμ=1+ReA1-ReA,where(43)A=limn→∞12Qn∑k=1n∫∂Dkφtdt.Here, φ(z) is a solution to the problem with n inclusions.In order to determineA, we first calculate the integral ∫∂Dkφ(z)dz. Equation (37) implies that (44)∫∂Dkφtdt=∑s=0∞∑j=0s-1-1s-j-1∫∂DkFjtdtρs.Application of (35) and Cauchy’s integral theorem yield (45)∫∂DkFszdz=∑m=1n∫∂DkΦmszm∗¯dz.The presented iterative scheme can be easily realized numerically. But we are interested in analytical formulae which can be obtained by symbolic computations. In the next sections, we perform symbolic computations to determine φ(z) and the effective shear modulus in the third-order approximation.


## 5.2. Second-Order Approximation

We are looking forφ(z) in the form (46)φz=φ1zρ+φ2zρ2+φ3zρ3+Oρ4.It follows from (37) that for the third-order approximation we need the integrals (45) for j=0,1,2. Using the second equation (32), we have (47)∫∂DkF0zdz=∑m=1n∫∂DkΦm0zm∗¯dz=-i∑m=1n∫∂Dkzm∗¯dz.The functions zm∗¯=r2/z-am+am¯ are analytic in Dk for m≠k. Therefore, the integrals with m≠k vanish in (47) by Cauchy’s integral theorem and (47) can be calculated by residues (48)∫∂DkF0zdz=-i∫∂Dkzk∗¯dz=2πr2.In order to use (37) for s=1 we find from (34) (49)Φk1z=iz+∑m≠kzk∗¯-ak¯-zm∗¯-am¯+ak¯-am¯r2z-am2.Using the expansion for m≠k, (50)z-am-l=∑j=0∞-1jl+j-1jak-am-l+jz-akj,we have(51)zk∗¯-ak¯r2z-am2=r2z-akr2z-am2=r4ak-am-2z-ak-1-r4∑j=0∞-1jj+2ak-am-j+3z-akj,(52)zm∗¯-am¯r2z-am2=r42∑j=0∞-1jj+1j+2ak-am-j+3z-akj.Along similar lines (53)ak¯-am¯r2z-am2=r2∑j=0∞-1jj+1jak¯-am¯ak-am-j+2z-akj.Subtracting (51) and (52) we obtain (54)zk∗¯-ak¯r2z-am2-zm∗¯-am¯r2z-am2=r4ak-am-2z-ak-1-r42∑j=0∞-1jj+2j+3ak-am-j+3z-akj.Using the above formulae we obtain from (49) the exact formula:(55)Φm1z=iz+∑m1≠mr4am-am1-2z-am-1-r42∑j=0∞-1jj+2j+3am-am1-j+3z-amj+r2∑j=0∞-1jj+1am¯-am1¯am-am1-j+2z-amj.We are now ready to calculate the integral: (56)∫∂DkF1zdz=∫∂Dk∑m=1nΦm1zm∗¯dz=-i∫∂Dk∑m=1nzm∗¯+∑m1≠mr4am-am1¯-2zm∗-am¯-1-r42∑j=0∞-1jj+2j+3am-am1¯-j+3zm∗-am¯j+r2∑j=0∞-1jj+1am-am1am-am1¯-j+2zm∗-am¯jdz.Using Cauchy’s integral theorem and residues we obtain(57)∫∂Dk∑m=1nzm∗-am¯-1dz=∫∂Dk∑m=1nz-amr2dz=0,∫∂Dk∑m=1nzm∗-am¯jdz=∫∂Dkr2z-akjdz=2πir2forj=10otherwise.This yields (58)∫∂DkF1zdz=2πr2+12πr6∑m≠kak-am¯-4-4πr4∑m≠kak-amak-am¯-3.


## 5.3. Third-Order Approximation

The third-order approximation requires advanced and long computations presented below. In order to calculateφ(3), we need the following function:(59)Φm2z=-i∑m1≠m∑m2≠m1zm1∗m2∗¯-∑m1=1nr2z-am1∑m2≠m1r2am1¯-am2¯2+∑m1≠mzm∗¯-am¯-zm1∗¯-am1¯+am¯-am1¯r2z-am12-∑m2≠m1r2r2/z-am1+am1¯-am2¯2+2r4r2/z-am1+am1¯-am2¯3z-am1-3r6r2/z-am1+am1¯-am2¯4z-am12+2r4am1-am2r2/z-am1+am1¯-am2¯3z-am12-∑m2≠m1r2am1¯-am2¯2.It is obtained from (34) for s=2 by substitution Φk(1) given by (49) and φk(1) calculated by (33) with s=1.Below, we describe general recurrent formulae for an arbitrarysth approximation and explicitly write the third-order terms when s=3. The following double series in ρ and (z-ak) are used: (60)φkz=∑s=0∞∑j=0∞φk,jsz-akjρs,(61)Φkz=r2z-akφm,1+∑j=0∞γm,jz-akj.The coefficients in the power series (61) are presented as series in ρ: (62)φm,l=∑s=0∞φm,lsρs,γm,l=∑s=0∞γm,lsρs.Substitution of (60)–(62) into (22) and selection of coefficients in the same powers of ρ and (z-ak) yield(63)φk,js=-δj,1φk,1s-1¯+-1j∑m≠k∑l=1∞γm,ls-1¯r2ll+j-1jak-am-l+jwith the zeroth approximation φk,j(0)=0. Along similar lines, (23) yields (64)γk,js=iδj,1δ0,s+δ1,s+∑m≠k--1jr2φm,1s-1¯ak-am-j+1+-1j∑l=1∞φm,ls-1¯r2ll+j-1jak-am-l+j+-1j∑l=1∞lγm,ls-1¯r2l+1l+j+2j+1ak-am-l+j+2--1jak-am¯∑l=1∞lγm,ls-1¯r2ll+jjak-am-l+j+1with the zeroth approximation γk,j(0)=iδ1,j.The expanded form ofφ(z) follows from (31):(65)φz=ρ1+ρ∑k=1n∑s=0∞∑j=0∞γk,js¯r2jz-ak-jρs.Equation (65) is in accordance with (37) where(66)Fjz=∑m=1n∑l=0∞γm,lj¯r2z-am-l.Using (66) we calculate the integral (67)∫∂DkFjtdt=∫∂Dk∑m=1n∑l=0∞γm,lj¯r2t-am-ldt=∫∂Dk∑l=0∞γk,lj¯r2t-ak-ldt=2πiγk,1j¯r2.The coefficientsγm,l(j) for j=0,1,2 can be explicitly written by the iterations (64): (68)γk,10=i,γk,11=i-ir22∑m≠kak-am¯ak-am3-6r2∑m≠k1ak-am4,(69)γk,12=i-r4∑k1≠k∑k2≠k11ak-ak121ak1¯-ak2¯2+6r4∑k1≠k1ak-ak14-2r2∑k1≠kak¯-ak1¯ak-ak13+∑l=1∞-1l-r2+2ll∑k1≠k∑k2≠k11ak-ak1l+11ak1¯-ak2¯l+1-14r2l+3ll+2l+32∑k1≠k1ak¯-ak1¯l+3∑k1≠k1ak-ak1l+3+12r2l+2ll+1l+2l+3∑k1≠kak-ak1ak¯-ak1¯l+2∑k1≠k1ak-ak1l+3+12r2l+2ll+1l+2l+3∑k1≠k1ak¯-ak1¯l+3∑k1≠kak¯-ak1¯ak-ak1l+2-r2l+1ll+12∑k1≠kak-ak1ak¯-ak1¯l+2∑k1≠kak¯-ak1¯ak-ak1l+2.The limit (43) has to be found. Using (44) and (67) we express A through γk,1(j): (70)A=∑s=0∞∑j=0s-1-1s-j-1πir2limn→∞1Qn∑k=1nγk,1j¯ρs.Taking into account terms up to O(ρ4) we rest γk,1(j) for j=0,1,2; hence, (71)A=πir2limn→∞1Qn∑k=1nγk,10¯ρ+πir2limn→∞1Qn∑k=1nγk,11¯-γk,10¯ρ2+πir2limn→∞1Qn∑k=1nγk,12¯-γk,11¯+γk,10¯ρ3+Oρ4.The above limits are calculated by use of the formalism developed in [8, 12]. It was supposed that N equal nonoverlapping disks belong to a parallelogram (fundamental) cell periodically extended to the complex plane by two linearly independent translation vectors. As an example, we shortly present the transformation of the term of γk,1(1) given by (68) multiplied by πir2 and complexly conjugated in accordance with (71): (72)Tn=12πr6∑m≠k1ak-am¯4.We are interested in the limit associated with the Eisenstein summation [8, 12]:(73)T=limn→∞Tn2Qn.The normalized e-sums were introduced in [8, 12] by means of the Eisenstein summation: (74)ep=limn→∞1N1+p/2Qnp/2∑m≠k1ak′-am′p,(75)epp=limn→∞1N1+pQnp∑m≠k∑k≠k11ak′-am′pak′-ak1′¯p,(76)ep1=limn→∞1N1+p-1/2Qnp-1/2∑m≠kak′-am′¯ak′-am′p,where ak′=Qn-1/2ak. This means that the rectangle Qn is normalized by the C-linear transformation z→z′=Qn-1/2z to Qn′ having the unit area. Therefore, T=6f3/π2e4¯. Other terms are transformed by the same method. As a result we arrive at the following formula:(77)A=fρ+2f3f2π2e4¯-fπe31¯ρ2+f-f2π2e2,2+∑l=1∞-1l-fπl+1lel+1,l+1-14fπl+3ll+2l+32el+32+12fπl+2ll+1l+2l+3el+21¯el+3+12fπl+2ll+1l+2l+3el+3¯el+21-fπl+1ll+12el+212ρ3+Oρ4,where (78)f=limn→∞Nπr2Qndenotes the concentration. The absolute convergence of series in R3 follows from the standard root test (Cauchy’s criterion).The effective shear modulus is calculated by (42):(79)μeμ=1+fρ+R2ρ2+R3ρ31-fρ-R2ρ2-R3ρ3,where(80)R2=23e4¯π2f3-e31¯πf2,(81)R3=ff2π2e2,2+∑l=1∞-1l-fπl+1lel+1,l+1-14fπl+3ll+2l+32el+32+12fπl+2ll+1l+2l+3el+21¯el+3+12fπl+2ll+1l+2l+3el+3¯el+21-fπl+1ll+12el+212.Expansion of (79) in ρ yields (82)μeμ=1+2fρ+2f2+R2ρ2+2f3+2fR2+R3ρ3+Oρ4.


## 6. Numerical Simulations

The asymptotic formula (82) is a new theoretical formula obtained in the present paper. The computationally effective formulae and algorithms for the absolutely convergent sums ep and epp (p>2) were developed in [13, 14]. However, the numerical implementation of the conditionally convergent series (74)-(75) for p=2 and (76) for p=3 requires further investigations.Formula (42) is similar to the famous Clausius-Mossotti approximation (Maxwell’s formula) [15, Section  10.4] applied for calculation of the effective conductivity of dilute composites. It is surprising that the Eisenstein summation and Maxwell’s self-consistent formalism are based on the different summation definitions of the conditionally convergent series [16, 17]. Using an analogy with the conductivity problem we now justify the proper limit value of e31, that is, (76) for p=3. The Hashin-Shtrikman bounds [15, Ch. 23] for 2D incompressible elastic media become (83)h+=1+fρ1-fρ,h-=ρ+1f-1ρ+1ρ-1f-1ρ-1forρ<0.In the case ρ>0, the bounds h+ and h- are replaced with each other. One can check that the upper and lower bounds (83) coincide up to O(ρ3): (84)h=1+2fρ+2f2ρ2+Oρ3.This implies that the coefficients in the term f2ρ2 of the expressions (82) and (84) must be equal. It follows from (80) that e31 must vanish in the framework of the considered Maxwell’s formalism. The definition of e2,2 is not essential in the final formulae since the terms with e2,2 are cancelled in (81). The above demonstration is based on the formal pure mathematical arguments. Its physical interpretation has been not clear yet.Consider a numerical example, the regular hexagonal lattice, when thee-sums become the Eisenstein-Rayleigh lattice sums ek(1)=Sk(1) and ek=Sk calculated by algorithm presented in [8, 18]. Only the nonzero values of Sk and Sk(1) are presented in Table 1.Table 1
The nonzero values of the lattice sums for the hexagonal array.
k
S
k
+
1
S
k
(
1
)
5
3.80815
4.2426
11
2.53532
2.93754
17
1.64402
1.89818
23
1.06787
1.23308
29
0.693603
0.800903
35
0.450508
0.520202
41
0.292614
0.337881
47
0.190058
0.21946
53
0.123446
0.142544
59
0.0801808
0.0925848The effective shear modulus for the regular hexagonal lattice becomes(85)μeμ=1+fρ+H3ρ31-fρ-H3ρ3+Oρ4,where (86)H3=ff2+∑l=1∞-1l-fπl+1lSl+12-14fπl+3ll+2l+32Sl+32+fπl+2ll+1l+2l+3Sl+21Sl+3-fπl+1ll+12Sl+212.The typical dependence of the effective modulus (85) on f is displayed in Figure 1 for ρ and in Figure 2 for negative ρ. A dependence on ρ is given in Figure 3. The presented and other graphs demonstrate sufficiently good precision for f<0.4 and for arbitrary |ρ|.Figure 1
Dependence of the effective shear modulus for the hexagonal array on the concentrationf for ρ=0.9. The data are computed by (85) (solid line) and its polynomial expansion up to O(ρ4) (dotted line). The Hashin-Shtrikman bounds (83) are shown by dashed lines.Figure 2
The same dependencies as in Figure1 but for ρ=-0.9.Figure 3
Dependence of the effective shear modulus for the hexagonal array onρ for f=0.3. The designations of lines are the same as in Figure 1. Equation (85) and its polynomial expansion coincide (solid line).


## 7. Conclusion

Schwarz’s alternating method is applied to 2D elastic problem for dispersed composites. It is realized for circular inclusions in symbolic form. Exact and approximate formulae for the local fields and for the effective shear modulus are established. In general, Schwarz’s method can be realized as expansion on the concentrationf and on the contrast parameter ρ. The concentration expansion was recently realized in [8]. In the present paper, we use the contrast parameter expansion. These two expansions in f and ρ yield two different computational schemes [7]. The effective modulus formulae are the same in the second-order approximations and begin to differ in the third-order terms in f and ρ, respectively. Schwarz’s method is used in the form of the functional equations method for circular inclusions [7–9].We develop a computationally efficient algorithm implemented in symbolic form to compute the local fields in 2D elastic incompressible composites and the effective shear modulus for macroscopically isotropic composites. The new analytical formula (79) contains the third-order term in ρ and explicitly demonstrates dependence on the location of inclusions. The theory is supplemented by a numerical example on the hexagonal array of inclusions. Figures 1–3 illustrate the dependence of the effective shear modulus on f and ρ. These numerical examples give sufficiently accurate values of μe for all admissible ρ, that is, for |ρ|≤1 and for f<0.4.One can expect that the precision ofμe will increase by using of the next approximation terms ρm (m>3). As it is noted in Introduction, Schwarz’s method can be based on two different expansion, in ρ used here and in f used in [12]. The both expansions contain the locations of inclusions in symbolic form. The expansions in f were held for smaller f but for an arbitrary ν [12]. At the present time, we are inclined to use the contrast parameter expansions. However, we suppose that the choice between different variants of Schwarz’s method will depend on the further implementation of high-order symbolic-numerical codes.



---

*Source: 4780928-2017-10-17.xml*