Effective elastic constants of hexagonal array of soft fibers

Computational Materials Science 139 (2017) 395–405

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Effective elastic constants of hexagonal array of soft fibers Piotr Drygas´ a, Simon Gluzman b, Vladimir Mityushev c,⇑, Wojciech Nawalaniec c a

Department of Mathematical Analysis, Faculty of Mathematics and Natural Sciences, University of Rzeszow, Pigonia 1, 35-959 Rzeszow, Poland 3000 Bathurst St, Apt. 606, Toronto, ON M6B 3B4, Canada1 c Pedagogical University, ul. Podchorazych 2, Krakow 30-084, Poland b

a r t i c l e

i n f o

Article history: Received 10 July 2017 Received in revised form 7 August 2017 Accepted 8 August 2017

Keywords: Effective elastic constants Hexagonal array Percolation regime Symbolic computation

a b s t r a c t Analytical formulae for the effective elastic constants of 2D composites with soft unidirectional fibers arranged in the hexagonal array are obtained for arbitrary concentration of fibers, from dilute case to percolation. It is supposed that every section of composites perpendicular to fibers is the hexagonal array with circular holes or soft inclusions. First, a polynomial approximation in concentration is obtained by application of functional equations. Further, an asymptotic analysis is applied to the obtained polynomial with using of the known asymptotic formulae in percolation regime. Finally, the formula for the effective shear modulus is suggested. It can be applied to the typical matrices made from ceramics, metals and polymers, and to other materials as well, e.g., for resin and thallium. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Accurate and consistent evaluation of the effective properties of composites and porous media is one of the fundamental tasks of applied mathematics and engineering. Experimental methods require advanced technological methods in order to obtain accurate results [1]. First of all this concerns the transverse shear modulus of fibrous composites. Therefore, theoretical investigations by means of analytical and numerical techniques are paramount. In particular, they are important for the regime with highconcentration of fibers. Solution to boundary value problems for a multiply connected domain in a class of periodic functions leads to theoretical evaluation of the effective properties. In the present paper, we restrict ourselves to plane strain elastic problems with circular holes that corresponds to fibrous composites with very soft fibers. Such elastic problems have application in poroelasticity, in particular, in biomechanics [2] in order to estimate a stressstrain state of reconstruction system bone-implant. Approximate and analytical formulae for porous media are applied in geophysics [3,4]. Applications of fiber composites in industry is outlined in [5]. Aluminium bricks with cylindrical holes were widespread in constructions of engines and it is still used in certain applications where it remains advantageous.

1

Present address.

⇑ Corresponding author. E-mail addresses: [email protected] (P. Drygas´), [email protected] (S. Gluzman), [email protected] (V. Mityushev), [email protected] (W. Nawalaniec). http://dx.doi.org/10.1016/j.commatsci.2017.08.009 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.

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 [6], first, extended to doubly periodic problems in [7] and developed in [8–10]. The obtained results were applied to computations of the effective properties of the elastic media. Integral equations are efficient for the numerical investigation of a non-dilute composites when interactions of inclusions have to be taken into account. Approximate analytical formulae were recently obtained for biLaplace’s equation which describes elastic materials [11] for a circular multiply connected domain. It is worth noting that many efforts were applied to get analytical formulae. The results were obtained for two extremal regimes when the concentration of inclusions f is low [12–15] and near the percolation threshold [12,15]. In the present paper, we apply the method of functional equations [16,17,11] to elastic plane problems for the regular hexagonal (triangular) array of holes. In this particular, but important in application case, we show how to deduce high order concentration formulae for general random 2D composites. The second step of investigations is asymptotic analysis to establish analytic formula for the effective shear modulus near the percolation threshold and to obtain a universal formula valid for all f. For an array of circular holes on the cites the hexagonal lattice effective shear modulus is expected to decay in the vicinity of p  0:9069 as a power-law, f c ¼ 12

le ’ Aðf c  f ÞT ;

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

396

with positive critical index T . The phase-interchange theorem [15], does not hold in the case of shear modulus for such array of holes, i.e., it should be two different values for the critical index for holes and inclusions, but does hold for the bulk-modulus [15]. For the antiplane shear problem the elastic displacements are zero in the plane but are non-zero in the direction perpendicular to the plane, and in the isotropic case only scalar elastic shear modulus l participates. The 2D dielectric constant (or electrical conductivity) problem is rigorously analogous to the antiplane shear elasticity problem [15,18], so that all results concerning effective properties discussed above can be applied to the effective shear l modulus l as a function of contrast parameter ql ¼ ll1  , and conþl

Dk

D

1

centration of the inclusions f. Here, l1 and l stand for the shear modulus of the inclusions and matrix, respectively. De Gennes [19] conjectured that in the vicinity of percolation threshold, such properties as conductance of a resistor network and the effective elastic modulus for one-component elastic displacement of point-like monomers, behave analogously. I.e., in antiplane shear an effective elastic modulus for perfectly rigid inclusions should diverge as

le  ðf c  f ÞS in the vicinity of the critical concentration f c . According to [19], the critical index S for the elastic modulus is equal to the superconductivity index s. By analogy we expect S ’ 1:3 in random case, and S ¼ 12 for regular lattice arrangements of inclusions. In case of plane strain elasticity such quantity as bulk modulus behaves similarly and with the same critical index at least in the regular case [15], and will be considered elsewhere. Another interesting exact result is independence of the effective Young modulus of a 2D sheet containing circular holes on the Poisson coefficient of the matrix [20,21]. Although it does not hold for rigid or other inclusions, actual dependence of the Poisson ratio m is rather weak. The holes are punched in the matrix with different elastic properties. Plane strain elastic problem is considered for such composite and the effective elastic modulus is obtained in the form of power series in the inclusions concentration and elastic constants for holes and matrix. For holes the shear modulus l1 ¼ 0 and the Poisson ratio m1 ¼ 0. We construct an expansion for the effective shear modulus avail14

able up to Oðf Þ inclusively.

Fig. 1. Section of fibrous composite.

rxx  ryy þ 2irxy

xx þ yy ¼

where the curves @Dk are oriented in clockwise sense. Let D be b (see Fig. 1). the complement of [nk¼1 ðDk [ @Dk Þ to C Let the uniform shear stress are applied at infinity

ð1Þ

The component of the stress tensor can be determined by the Kolosov-Muskhelishvili formulae [6]



4Re u0k ðzÞ; z 2 Dk ; 4Re u0 ðzÞ;

z 2 D;

( j1 1

l1 Re uk ðzÞ; j1 Re u0 ðzÞ; l 0

(

xx  yy þ 2ixy ¼

z 2 Dk ;

ð4Þ

z 2 D;

   l1 Re zu00k ðzÞ þ w0k ðzÞ ; z 2 Dk ; 1    l1 Re zu00 ðzÞ þ w00 ðzÞ ; z 2 D;

ð5Þ

m ; m is the 2D Poisson ratio [22]. The same notation is where j ¼ 3 1þm used for j1 . We have w0 ðzÞ ¼ iz þ wðzÞ. The functions uðzÞ; wðzÞ are analytic in D, twice differentiable in D [ @D and bounded at infinity. The functions uk ðzÞ and wk ðzÞ are analytic in Dk and twice differentiable in the closures of the considered domains. The perfect bonding at the matrix-inclusion interface can be expressed by two conditions [6]

uk ðtÞ þ tu0k ðtÞ þ wk ðtÞ ¼ uðtÞ þ tu0 ðtÞ þ w0 ðtÞ;  l  j1 uk ðtÞ  tu0k ðtÞ  wk ðtÞ ¼ 1 juðtÞ  tu0 ðtÞ  w0 ðtÞ : l

ð6Þ ð7Þ

Introduce the new unknown functions

We begin our study with a finite number n of inclusions on the infinite plane. This number n is given in a symbolic form with an implicit purpose to pass to the limit n ! 1 later. The shear modulus of inclusions is also arbitrary taken as l1 to pass to the limit l1 ! 0 in the final formulae. Introduce the complex variable z ¼ x þ iy where i denotes the imaginary unit. Let inclusions be disks Dk ¼ fz 2 C : jz  ak j < rg ðk ¼ 1; 2; . . . ; nÞ in the extended b ¼ C [ f1g. Denote @Dk :¼ ft 2 C : jt  ak j ¼ rg complex plane C

rxx þ ryy ¼

ð3Þ

where Re denotes the real part and the bar the complex conjugation. The component of the deformation tensor have the form

2. Method of functional equations for local fields

1 1 r1 r1 xx ¼ ryy ¼ 0; xy ¼ ryx ¼ 1:

8 h i > < 2 zu00k ðzÞ þ w0k ðzÞ ; z 2 Dk ; h i ¼ > : 2 zu00 ðzÞ þ w00 ðzÞ ; z 2 D;

ð2Þ



Uk ðzÞ ¼

 r2 þ ak u0k ðzÞ þ wk ðzÞ; z  ak

jz  ak j 6 r;

analytic in Dk except the point ak where its principal part has the form r 2 ðz  ak Þ1 u0k ðak Þ. Introduce the inversion with respect to the circle @Dk

zðkÞ ¼ r2 ðz  ak Þ1 þ ak : The problem (6), (7) was reduced in [16] (see Eqs. (5.6.11) and (5.6.16) in Chapter 5), [23] to the system of functional equations









i Xh l1 l1 þ j1 uk ðzÞ ¼ 1 Um ðzðmÞ Þ  ðz  am Þu0m ðam Þ l l m–k   l1   1 u0k ðak Þðz  ak Þ þ p0 ; jz  ak j 6 r; l ð8Þ

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405



  X l l j 1 þ 1 Uk ðzÞ ¼ j 1  j1 um ðzðmÞ Þ

l

l m–k   X 2  l1 r r2 1 þ ak   am þ l z  ak z  am m–k h 0 i l  0  Um ðzðmÞ Þ  um ðam Þ þ 1 ð1 þ jÞiz l þ

n X r 2 qk þ q0 ; z  ak k¼1

jz  ak j 6 r;

k ¼ 1; 2; . . . ; n; ð9Þ

where p0 and q0 are constant and

    l l qk ¼ u0k ðak Þ ðj  1Þ 1  ðj1  1Þ  u0k ðak Þ 1  1 ;

l

l

k ¼ 1; 2; . . . ; n:

ð10Þ

The unknown functions uk ðzÞ and Uk ðzÞ ðk ¼ 1; 2; . . . ; nÞ are related by 2n Eqs. (8) and (9). The functional equations contain compositions of uk ðzÞ and Uk ðzÞ with inversions which determine a compact operator in a functional space [24]. The functions uðzÞ and wðzÞ are expressed by means of uk ðzÞ and wk ðzÞ

 X n h i l1 l1 ð1 þ jÞuðzÞ ¼ 1 Um ðzðmÞ Þ  ðz  am Þu0m ðam Þ þ p0 ; l l m¼1 ð11Þ n X l1 r 2 qk ð1 þ jÞwðzÞ ¼ þ q0 l z  ak k¼1  X h n  0 i l1 r2  1 þ am Um ðzðmÞ Þ  u0m ðam Þ l z  am m¼1  X n l1 þ j  j1 um ðzðmÞ Þ; z 2 D: ð12Þ l m¼1

The additive constants p0 ; q0 do not impact on the final formulae, hence they will be omitted in the forthcoming computations.

397

We use the computational scheme to obtain an approximate analytical expression for le shortly outlined as follows. First, the complex potential uðzÞ for a finite number of inclusions is calculated analytically by (11) where um ðzÞ and Um ðzÞ are approximately found from the system of functional Eqs. (8) and (9) as series in r2 . The passage to the limit in (16) for statistically homogeneous composites can be accomplished for doubly periodic composites. The justification of the formal scheme presented below was given in [11] by means of the Eisenstein summation method. The integral (16) can be calculated explicitly up to practically arbitrarily fixed precision in r2 by using approximations of the function uðzÞ. As an example, we consider now the precision Oðr2 Þ. Introduce the constant

c1 ¼

l1 l 1 : j ll1 þ 1

ð17Þ

Then, the functional Eqs. (8) and (9) yield up to Oðr 4 Þ

uk ðzÞ ¼ 0; Uk ðzÞ ¼ ic1 z:

ð18Þ

Substitution (18) into (11) yields uðzÞ ¼ ic1

Pn

r2 k¼1 zak

up to Oðr4 Þ.

The next terms Oðr4 Þ are calculated by the same method. Substituting (18) into the right hand part of (8) and (9) we get simple but long formulae for uk ðzÞ and Uk ðzÞ within Oðr4 Þ. After substitution of the obtained expressions into (11) we arrive at the approximate analytical formula

uðzÞ ¼ ic1

n X

n X r2 r4 X am  ak þ 2ic21 þ Oðr 6 Þ: z  ak z  ak m–k ðam  ak Þ3 k¼1 k¼1

ð19Þ

Applying the residue theorem to the integral (16) with uðzÞ given by (19) we obtain

1 2

Z @Dk

6 uðtÞdt ¼ pr2 c1  2pr4 c21 eð1Þ 3 ðnÞ þ Oðr Þ;

ð20Þ

where 3. Effective constants ð1Þ

After approximate solution to the functional equations following [11] we can find the averaged shear modulus le of the considered finite composite with n inclusions on the plane. The average over a sufficiently large rectangle Q n containing all the inclusions Dk is introduced by means of the integral

e3 ðnÞ ¼

X am  ak m–k ðam

 ak Þ3

:

ð21Þ

ðnÞ

hwin ¼

1 jQ n j

Z

w dx1 dx2 :

ð13Þ

Qn

Then, the macroscopic shear modulus can be computed by formula

lðnÞ e ¼

hrxy in : 2hxy in

ð14Þ

Further, Q n is extended to the infinite point (n ! 1) and we arrive at the macroscopic shear moduli

le 1 þ Re A ¼ ; l 1  jRe A

le ¼ limn!1 leðnÞ and ð15Þ

n Z 1 X uðtÞdt: n!1 2jQ n j k¼1 @Dk

Introduce A ¼ fpx1 þ qx2 ; p; q 2 Zg. One can consider the set A as the reordered set of points fa1 ; a2 ; . . .g. It was proved in [25] that for the hexagonal array

1 ð1Þ lim e ðnÞ n!1 jQ n j 3 Substitution

¼

p 2

of

:

ð22Þ (20)

and

(22)

into

(16)

yields

pr denotes the concenRe A ¼ f c1  f c21 þ Oðf Þ, where f ¼ limn!1 jQ nj 2

3

2

tration of inclusions. Then, (15) implies that

where

A ¼ lim

It is problematic to calculate without computer the next approximations for uðzÞ because of huge manipulations. High order formulae obtained by use of the advanced symbolic computations will be described in the next section. Let the set of centers of disks generate the hexagonal double periodic structure formed by two fundamental translation vectors qffiffi qffiffi pffiffi expressed by complex numbers x1 ¼ 4 43 and x2 ¼ 4 43 12 þ i 23 .

ð16Þ

Eqs. (15) and (16) are obtained by substitution of the representations (2)–(5) into (14) and by use of the Cauchy integral theorem for analytic functions (for details see [11]).

2 le 1 þ f c1 le 1 þ f c1  f c21 2   þ Oðf 3 Þ: ¼ þ Oðf Þ; ¼ l 1  jf c1 l 1  j f c1  f 2 c2 1

ð23Þ

Remarkably the first order approximation (23) coincides with the Hashin-Shtrikman upper bound lþ for l1 < l (take into account (17))

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

398

lþ

  1 þ j ll1  f 1  ll1  ; ¼l 1 þ j ll1 þ f j 1  ll1

l ¼ l1

  1 þ j1  f 1  ll1  ; ð1 þ j1 Þ ll1 þ f j1 1  ll1 ð24Þ

where the Hashin-Shtrikman lower bound completeness [18].

l is written for

Hereafter, the shear modulus of matrix is normalized to unity, i.e.,

l ¼ 1. The coefficients ak ðmÞ are written in terms of the Poisson ratio a1 ðmÞ ¼  a3 ðmÞ ¼ 

4. High order formulae The general iterative scheme to compute the effective shear modulus is outlined in the previous section. The results were writ4

ten in [11] up to Oðf Þ for the general 2D random locations and up

N1 X uk ðzÞ ¼ r2s ukðsÞ ðzÞ þ Oðr2N Þ;

ð25Þ

s¼0

The functions uk and wk ðs ¼ 0; 1; . . .Þ in the kth inclusion are presented by their Taylor polynomials ðsÞ

u

ðsÞ k ðzÞ

¼

M1 X

a

ðsÞ k;j ðz

ðsÞ

j

M

 ak Þ þ Oððz  ak Þ Þ;

ð26Þ

j¼0

ðsÞ

wk ðzÞ ¼

M1 X

ðsÞ

as z ! ak :

ðs1Þ

ðs1Þ

and bk;j

ðm þ 1Þ

7

a9 ðmÞ ¼

1 ðm þ 1Þ9

14 X k 15 ak ðmÞf þ Oðf Þ: k¼1

;

193:889m5 þ 594:385m4  803:573m3

418:671m6 þ 1654:74m5  2509:19m4

 2236:51m8 þ 8570:75m7  32965m6 þ 27094:3m5

1 ðm þ 1Þ

10

 4630:25m9 þ 26308:2m8  60019:5m7

þ 191379m6 111054m5 þ 523896m4  13463:8m3  þ 108465m2 þ 10131:4m þ 19491:6 a11 ðmÞ ¼

1 ðm þ 1Þ

11

 10363:1m10 þ 63337:9m9  179675m8

þ 451414m7 774552:m6 þ 833885m5  1:7462  106 m4 :  325402m3  470665m2 66131:5m  55740:6Þ a12 ðmÞ ¼

ðsÞ

ð28Þ

;

  105814m4 þ56288:1m3  23233:8m2  1160:52m  7102:14

1 ðm þ 1Þ12

 22724:3m11 þ 157024m10  478055m9

þ 1:24987  106 m8 2:5037  106 m7 þ 3:08447  106 m6

14

le ¼ 1 þ

ðm þ 1Þ4

ðm þ 1Þ  þ 15943:8m4  25066:9m3 þ9735:43m2 þ 1108:83m þ 2723:74 8

.

complex potentials uk and wk up to Oðf Þ. The integral (16) contains the function u in the form of the series in the concentration f and is calculated by the resides also within the package Mathematica in symbolic form. The limit n ! 1 is calculated by using of the formalism based on the Eisentein summation (for details see [11,26]) implemented in the symbolic algorithm. Analytical formulae include symbolically the elastic constants j; l j1 ; l1 and with geometrical parameters r; n; ak but in the form of huge expressions. These long formulae can be analyzed by substitution of numerical values. Below, the number of symbolic parameters is reduced including the special hexagonal geometry in order to present and check the obtained formulae for the effective constants. In further computations the physical properties of the inclusions are taken to be j1 ¼ 1 and l1 ¼ 0 that corresponds to empty fibers. Moreover, the constant j of matrix is represented m. by the Poisson ratio m as j ¼ 3 1þm First, we write the effective shear modulus

32ðm  1Þ3

 899:127m7 þ 4365:01m6  7910:82m5

1

of the coefficients ak;j and bk;j for s ¼ 1; 2; . . . ; 13 which yield the ðsÞ

a4 ðmÞ ¼ 

;

þ 6002:38m3  5550:21m2 778:073m  1026:87Þ

ðsÞ

We implemented the above scheme in the Mathematica computer algebra system, in order to obtain symbolic-numerical form



1

a7 ðmÞ ¼

result we get an iterative scheme on s when ak;j and bk;j are computed explicitly by ak;j

;

ðm þ 1Þ2

þ 2324:08m2 þ165:871m þ 345:94Þ;

j¼0

ðsÞ



ðm þ 1Þ6

ð27Þ

Substitute of (25)–(27) in (8) and (9). First, selecting the terms with the same powers of ðz  ak Þ we arrive to a system of linear algebraic equations on the coefficients of the functions uk and wk . Next, we select the terms of the same powers r 2s in these coefficients. As a

ðm þ 1Þ3

1

a6 ðmÞ ¼

a10 ðmÞ ¼ bk;j ðz  ak Þ j þ Oððz  ak ÞM Þ;

16ðm  1Þ2

8ðm  1Þ

ðm þ 1Þ5

a8 ðmÞ ¼

s¼0 N1 X ðsÞ wk ðzÞ ¼ r2s wk ðzÞ þ Oðr 2N Þ

a2 ðmÞ ¼ 

;

99:4786m4 þ 114:086m3  596:872m2 þ 114:086m  99:4786

a5 ðmÞ ¼

8

to Oðf Þ for the regular hexagonal array (see Supplementary to [11]). This general scheme can be improved by modification of the theoretical algorithm from [11] as follow. We are looking for the complex potentials uk and wk in the form

4

mþ1

 4:93632  106 m5 þ 3:90233  106 m4 þ1:00308  106 m3  þ 1:55173  106 m2 þ 264862m þ 158929 

1

a13 ðmÞ ¼

ðm þ 1Þ

13

52537:1m12 þ 349702:m11  1:4589  106 m10

þ 2:77236  106 m9  9:02833  106 m8 þ 8:8219  106 m7  1:77656  107 m6 þ1:71112  107 m5  7:22178  106 m4  1:37692  106 m3  4:28978  106 m2 830021:m  443663:Þ

a14 ðmÞ ¼



1 ðm þ 1:Þ

14

112950:m13 þ 912406:m12  3:17506  106 m11

þ 1:08494  107 m10  1:5547  107 m9 þ 5:4004  107 m8  2:23313  107 m7 þ1:03973  108 m6  2:49654  107 m5 þ 2:49512  107 m4 þ 2:96962  106 m3 þ1:17975  107 m2  þ 2:43867  106 m þ 1:22415  106 ð29Þ The final formulae (29) are written in terms of tigate the corresponding dependencies.

m in order to inves-

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

4.1. Negative

m of the matrix

Consider first another material case of circular holes embedded into the matrix with small negative Poisson coefficient. Foams can indeed possess the negative Poisson ratio [27,28], but other materials can also exhibit such property [29]. Such materials, known as auxetics, are of great practical interest because of their numerous potential applications, such as the design of fasteners, prostheses, piezo-composites, filters, earphones, seat cushions and superior dampers. Let us set l1 ¼ 0; m1 ¼ 0; m ¼ 0:05, leading to the following 16

series for the effective shear modulus computed up to Oðf Þ

le ðf Þ ¼ 1  4:21053f þ 9:30748f 2  20:5744f 3 þ 45:4803f 4 5

6

9

10

7

 137:881f þ 467:371f  1435:72f þ 4063:6f  11281:1f þ 32169:7f þ 276604f

12

 804073f

 6:69609  106 f

15

13

 94200:3f

and to the new series L1 again apply iterated roots Rð1Þ n , so that a new sequence of approximations for the critical index naturally arises, 2 ð1Þ T ð1Þ n ¼ T 1 lim ð1 þ z Rn ðzÞÞ;

11

where roots behave as Rð1Þ  z2 , as z ! 1, in order to extract a n constant/correction to T 1 . It turns that we can construct only two approximants

2:44524 ; ð1:09546z þ 1Þ2 0:166188 ð1Þ : R2 ðzÞ ¼  2 z  0:148903z  0:0679637 ð1Þ

R1 ðzÞ ¼

ð1Þ

þ 2:32021  106 f

14

ð2Þ

þ :

ð30Þ

Suppose that after the D  Log transformation applied to original series (30) and presented in terms of variable z, we are confronted with the problem of extraction the critical index from transformed series LðzÞ. Suppose also we would like to find the approximations to index through the sequence of iterated roots Rn ,

T n ¼ lim ðz Rn ðzÞÞ;

ð31Þ

z!1

which behaves as Rn  z1 , as z ! 1. Naturally one would expect to be able to construct any desired number of root approximants and obtain a well behaving sequence of approximations. But what to do if the sequence ends abruptly already at the first term R1 ðzÞ? This is precisely the case of series (30). We can only find

3:81852 ; 1 þ 2:19093z

no further improvement to the value of index, with T 2 ¼ 1:33475. Let us now look for the solution in the form of a factor approximant [33], also satisfying the two starting non-trivial terms from the series (30),

l0 ðf Þ ¼ ð1  1:10266f Þ1:33475 ð1 þ f Þ2:73876 ;

and corresponding index T 1 ¼ 1:74288. This estimate may be good, most often it is not. Nevertheless, one should attempt to perform more approximation steps.

ð2Þ

ð2Þ

T n ¼ T 2  lim ðz Pn;nþ1 ðzÞÞ:

1

ð37Þ

z!1

The corrected values for the critical index can be calculated readily and we have now reasonable estimate, T 3 ¼ 1:55176. Effective shear modulus can be reconstructed using the complete expression for the effective critical index, employing the approximant P 3;4 ðzÞ, but with additional requirement that the critical index equals precisely 32,

P3;4 ðzÞ ¼

22:6942z3  11:3617z2  0:18873z : 137:332z4  29:5068z3 þ 2:72942z2 þ 6:78983z þ 1 ð38Þ

and

l

 e ðf Þ

¼l

 0 ðf Þ exp

Let us construct the new series L1 ðzÞ ¼ RLðzÞ ,  ðzÞ

Z

f f c f

! P 3;4 ðzÞdz :

ð39Þ

0

The integral can be found analytically,

L1 ðzÞ ’ 1 þ b2 z þ b3 z þ b4 z þ b5 z þ    4

ð36Þ

so that our zero approximation for the critical index is equal to T 2 . Let us divide then the original series (30) by l0 ðf Þ, express the newly found series in term of variable z, then apply D  Log transformation and call the transformed series LðzÞ. Applying now the Padé approximants P n;nþ1 ðzÞ, one can obtain the following sequence of corrected approximations for the critical index [30–32],

ð32Þ

3

ð35Þ

T 2 ¼ 1:45324. This estimate for the index is in line with other estimates presented above. The procedure can be repeated, but there is

zfc ; to the original series. z ¼ f cff () f ¼ zþ1

2

ð34Þ

z!1

The first of them gives meaningless result, but the second approximant does better, leading to the new estimate for the index

8

The D-Log-Padé method basically fails in application to this series. Only through the method of corrected approximants applied to calculation of the critical index, we are able to achieve sensible results. The strategy pursued below can be tried for very difficult resummation problems. Let us introduce also the following transformation,

R1 ðzÞ ¼ 

399

5

le ðf Þ ¼ l0 ðf ÞFðf Þ;

¼ 1 þ 2:44524z2  5:35735z3 þ 47:7162z4  243:019z5 þ . . . :

ð33Þ

ð40Þ 0:16525

with positive power ð0:9069  f Þ ing in the following formulae,

coming from Fðf Þ, and result-

2 0:01256

0:16525

Fðf Þ ¼ ð0:9069  f Þ

ð0:04580 þ 0:3924f þ f Þ

2 0:09519

ð0:09019  0:4954f þ f Þ  0:0125649      ff þ0:392461Þþ0:0457968 2:7292 14:3374 0:824496 exp 0:1609 tan1 3:8827  0:9069f þ 0:07807 tan1 12:9196  0:9069f ð0:9069f Þ2  ;  0:0951904 f ðf 0:495435Þþ0:0901874 ð0:9069f Þ2

ð41Þ

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

400

approximant P 3;4 ðzÞ, but re-calculated with additional requirement that the critical index equals precisely 32,

or with explicitly extracted power-law

ðx2 þ 0:3925x þ 0:04580Þ

le ðf Þ ¼ 0:9394ð0:9069  xÞ2 3

0:01256 0:09519

ðx þ 1Þ2:7388 ðx2  0:4954x þ 0:09019Þ   2:7292  exp 0:1609tan1 3:8827  0:9069  x   14:3374 : ð42Þ þ 0:07807tan1 12:9196  0:9069  x

One also finds that the critical amplitude is equal to A ¼ 0:118849. Importantly, formula (40) satisfy the HashinShtrickman bounds. 4.2. Hexagonal array of circular holes. Positive

m of the matrix

Fðf Þ ¼ 1:02626

0:0511627  0:421789f þ f

2

0:02092 

19:7921z3  5:257z2  1:55789z :  59:4359z3  11:8029z2 þ 3:18517z þ 1

249:905z4

ð46Þ and

l

 e ðf Þ

¼l

Z

 0 ðf Þ exp

f f c f

! P 3;4 ðzÞ dz :

ð47Þ

0

The integral can be found analytically, so that

le ðf Þ ¼

Consider yet different material case. Set l1 ¼ 0; m1 ¼ 0; m ¼ 0:2. From the general expression we obtain the following series



P3;4 ðzÞ ¼

ð1  1:10266f Þ

1:5792

1:59202

ðf þ 1Þ

Fðf Þ;

ð48Þ

with

0:0447654 þ 0:328829f þ f

2

0:01868

0:0792



  exp 0:064 tan1 8:5119 þ

ð0:9069  f Þ    6:0061 8:7518  0:0169 tan1 8:0449 þ : f  0:9069 f  0:9069

ð49Þ

le ðf Þ ¼ 1  3:33333f þ 4:44444f 2  5:92593f 3 þ 7:90123f 4 5

6

7

 1545:79f þ 4290:08f

9

10

 12512:4f

11

þ 32169:6f

13

14

 443302f

15

þ 

 40:1005f þ 156:243f  379:549f þ 734:365f  74714:9f

þ 174931f

8 12

ð43Þ

Using the D  Log-Padé method one can obtain a two estimates for the critical index T 2 ¼ 1:5792; T 6 ¼ 1:65827. But one can not reconstruct the approximant corresponding to the first of estimates, because it is not a holomorphic function. One can use the result of D  Log-Padé method as an initial guess for some other method, i.e. corrected approximants. Let us look for the solution in the form of a factor approximant [33], also satisfying the two starting non-trivial terms from the series (43),



l0 ðf Þ ¼ ðf þ 1Þs1 1 

f fc

1:5792 ;

ð44Þ

so that our zero approximation for the critical index is equal to T ð0Þ ¼ 1:5792. Parameter S 1 ¼ 1:59202, as follows from the asymptotic equivalence with (43). Let us divide then the original series (43) by l0 , express the newly found series in term of variable z, then apply D  Log transformation and call the transformed series LðzÞ. Applying now the Padé approximants Pn;nþ1 ðzÞ, one can obtain the following sequence of corrected approximations for the critical index,

Tn¼T

ð0Þ

 lim ðzP n;nþ1 ðzÞÞ: z!1

ð45Þ

The ‘‘corrected” values for the critical index can be calculated readily and we have now even better estimates, T 2 ¼ 1:3722; T 3 ¼ 1:50746. The last estimate is generated by the unique among others, holomorphic approximant. Effective shear modulus can be reconstructed using the complete expression for the effective critical index, employing the

The critical amplitude A ¼ 0:393198 is calculated from (48). Also formula (48) fits perfectly within the Hashin-Shtrickman bounds. The Standard Padé technique does not work failing to produce a holomorphic approximant except in very high orders, signaling perhaps that much longer series are needed to synchronize performance of different techniques. 4.3. Holes in the incompressible matrix Let

us

set

yet

different

material

parameters

l1 ! 0; m1 ¼ 0; m ! 1. For this ‘‘material” case with holes within the incompressible matrix, the effective properties do not depend on the elastic properties of the components. For corresponding series

le ðf Þ ¼ 1  2f  17:7393f 5 þ 38:0128f 6  20:5148f 7 9

 157:341f þ 674:32f

10

þ 837:767f

12

 1637:25f

 25038:2f

15

þ 

 1113:33f 13

11

þ 8971:47f

14

ð50Þ

the coefficients increase pretty fast. Let us start with formulating the initial approximation, and recast it as



l0 ðf Þ ¼ 1 

f fc

T 0 Rðf Þ;

ð51Þ

where T 0 ¼ 2 is the initial guess for the critical index, and c Rðf Þ ¼ ð1 þ b1 f Þ 1 stands for the regular part. We set b1 ¼ 2 as the absolute value of a linear term in series, and discover from asymptotic equivalence that c1 ¼ 0:102658 [33]. In what follows we attempt to correct l0 ðf Þ, assuming instead of T 0 some functional dependence T ðf Þ, not unlike [32]. As f ! f c ; T ðf Þ ! T c , the corrected value for the critical index. The function T ðf Þ will be

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

designed in such a way, that it smoothly interpolates between the initial value T 0 and the sought value T c . The corrected functional form for the shear modulus is now



le ðf Þ ¼ 1 

f fc

T ðf Þ Rðf Þ:

ð52Þ

From (52) one can express T ðf Þ, but only formally since the exact form of le ðf Þ is not known. But we can use its asymptotic form (50), express T ðf Þ as a series and apply some resummation procedure (e.g. Padé technique). Finally calculate the limit of the approximants as f ! f c to obtain the critical index.

0:0682627zðz  0:638698Þðz þ 0:573096Þ ðz2  0:721815z þ 0:166311Þðz2 þ 0:406436z þ 0:0693321Þ ð57Þ  f  R and le ðf Þ ¼ l0 ðf Þ exp 0f c f P3;4 ðzÞdz . The integral can be found

P 3;4 ðzÞ ¼

analytically,

le ðf Þ ¼ 0:884629ð2f þ 1Þ0:1354 ð0:9069  f Þ0:06826 ð1  1:1027f Þ1:56826  0:03094 2  f  0:506464f þ 0:0724451  0:003196 2  f þ 0:366335f þ 0:0860216

e ðf Þ In what follows the ratio Cðf Þ ¼ lRðf , stands for an asymptotic Þ

form of the singular part of the solution, and as f ! 0

log ðCðf ÞÞ  ; T ðf Þ ’ log 1  ff

   exp 0:108285tan1 7:16694 þ

 4:776 f  0:9069   5:41642 :  0:0990339tan1 4:75874 þ f  0:9069

ð53Þ

c

which can be easily expanded in powers f, around the value of T 0 . 2

T ðf Þ ’ 2 þ 0:524941f þ 1:56668f þ 1:50835f 4

5

8

9

6

þ 78:3672f þ 267:403f þ 166:459f þ 729:203f

11

þ 1385:3f

12

þ :

7

4.4. General expression for the effective shear modulus for holes

10

ð54Þ

It appears that one can construct three meaningful Padé approximants, and find corresponding limits T 1 ðf c Þ ¼ 1:72105; T 2 ðf c Þ ¼ 1:56826; T 3 ðf c Þ ¼ 1:67526. Note that if a truncated polynomial of 8th order is used instead of the 15th order, there appear two more good estimates for the index, T 5 ðf c Þ ¼ 1:52692; T 6 ðf c Þ ¼ 1:47219. The 8th order polynomial is qualitatively similar to the 15th order but gives better value for the threshold, hence its better performance with regard to the critical index. The corrected index gets close to the exact 32 [34,15]. But corresponding Padé approximants are not holomorphic and do not represent the effective modulus correctly in the whole region of concentrations. Now we attempt to reconstruct a complete expression for the effective shear modulus. Employ the same idea as above suggested for inclusions. Let us start with different initial approximation,



l0 ðf Þ ¼ 1 

f fc

T 0 Rðf Þ;

ð55Þ

where T 0 ¼ 1:56826 is the new initial guess for the critical index borrowed from the previous stage of consideration. Here c Rðf Þ ¼ ð1 þ b1 f Þ 1 stands for the regular part. We again set b1 ¼ 2, and discover from asymptotic equivalence that c1 ¼ 0:135371. Let us divide then the original series (50) by l0 , express the newly found series in term of variable z, then apply D  Log transformation and call the transformed series LðzÞ. Applying now the Padé approximants Pn;nþ1 ðzÞ, one can obtain the following sequence of corrected approximations for the critical index,

T n ¼ T 0  limðz P n;nþ1 ðzÞÞ: z!1

ð58Þ

Formula (58) respects the Hashin-Shtrickman bounds (24) everywhere.

3

þ 20:2539f  6:61928f þ 26:11f þ 34:0555f

401

Consider an interpolation problem, when both lowconcentration series and high concentration formula are known. Amplitude A in the high-concentration limit is known [15] and the following analytical expression is valid,

AðmÞ ¼

4ðm þ 1Þ 32 pffiffiffi f c : 3p 3

The simplest root approximant satisfying the first non-trivial term from the series as well as the high-concentration limit is readily written in terms of variable zðf Þ ¼ f cff ,



2

2 ð2Þ l2 ðf ; mÞ ¼ ðpð2Þ 1 ðmÞzðf Þ þ 1Þ þ p2 ðmÞðzðf ÞÞ

where

2:4184 ð2Þ p1 ðmÞ ¼ ; mþ1

p2 ðmÞ ¼ ð2Þ

3=4

;

ð60Þ

  6:52172 ðm þ 1Þ2=3  0:896796 ðm þ 1Þ2

:

The next order iterated root approximant satisfying two non-trivial terms from the series as well as the high-concentration limit can be written following standard procedure,

1

ffi; l3 ðf ; mÞ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 2 2 2 3 ð3Þ ð3Þ ð3Þ þ p3 ðmÞðzðf ÞÞ ðp1 ðmÞzðf Þ þ 1Þ þ p2 ðmÞðzðf ÞÞ ð61Þ where p1 ðmÞ ¼ ð3Þ

ð56Þ

The corrected values for the critical index can be calculated readily and we have now two reasonably close estimates, T 2 ¼ 1:32585; T 3 ¼ 1:4837, appearing already in low orders. Effective shear modulus corresponding to T 3 can be reconstructed using the complete expression for the effective critical index, employing holomorphic approximant P3;4 ðzÞ, but corrected with additional requirement that the critical index equals precisely 32,

ð59Þ

¼

2:4184 ; mþ1 3:93618

p2 ðmÞ ¼ ð3Þ

3:93618m þ 1:01186

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3:93618mþ6:86051 ðmþ1Þ2

ðm þ 1Þ2

m  6:86051 ðm þ 1Þ

p3 ðmÞ ð3Þ

;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3:93618mþ6:86051 ðmþ1Þ2

þ 16:655

2

;

The higher order iterated root approximant satisfying three non-trivial terms from the series as well as the highconcentration limit can be written following standard procedure, l4 ðf ; mÞ ¼



2

ðp1 ðmÞzðf Þ þ 1Þ þ p2 ðmÞðzðf ÞÞ ð4Þ

ð4Þ

2

3=2

þ p3 ðmÞðzðf ÞÞ ð4Þ

4=3 3

þp4 ðmÞðzðf ÞÞ ð4Þ

4

3=8

;

ð62Þ

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

402

where

p1 ðmÞ ¼  ð4Þ

pa1 ðmÞ

pffiffiffi ; 3 3  pffiffiffi 1  2 ð4Þ p2 ðmÞ ¼ 5p ða1 ðmÞÞ2 þ 12 3pa1 ðmÞ  6p2 a2 ðmÞ ; 54 ð4Þ p3 ð

pffiffiffi 1  pffiffiffi 3 63 3p a1 ðmÞa2 ðmÞ  35 3p3 ða1 ðmÞÞ3 mÞ ¼ 972 pffiffiffi  378p2 ða1 ðmÞÞ2 324 3pa1 ðmÞ  pffiffiffi þ 324p2 a2 ðmÞ  27 3p3 a3 ðmÞ ;

ð4Þ



pffiffiffi 324a1 ðmÞ  126 3pða1 ðmÞÞ2  35p2 ða1 ðmÞÞ3 2916 2 pffiffiffi þ 108 3pa2 ðmÞ þ 63p2 a1 ðmÞa2 ðmÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffi þ 12 6pa1 ðmÞ 7pða1 ðmÞÞ2 þ 12 3a1 ðmÞ  6pa2 ðmÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q pffiffiffi pffiffiffi þ 7 2p3=2 ða1 ðmÞÞ2 7pa1 ðmÞ2 þ 12 3a1 ðmÞ  6pa2 ðmÞ 4=3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi pffiffiffi 6 2p3=2 a2 ðmÞ 7pða1 ðmÞÞ2 þ 12 3a1 ðmÞ  6pa2 ðmÞ  27p2 a3 ðmÞ

p4 ðmÞ ¼ 

þ

p4=3

2=3

144

p4 AðmÞ8=3

:

Good quality of interpolation, comparable to that of iterated roots, is achieved also by the following Padé approximant qffiffiffiffiffiffiffiffiffiffiffiffiffiffi f 0:723377 0:9069f þ1  qffiffiffiffiffiffiffiffiffiffiffi lp ðf ; mÞ ¼  : f qffiffiffiffiffiffiffiffiffiffiffiffiffiffi f 2:62412 3:6276 0:9069f f 2:95214f 2 þ 0:723377 0:9069f þ ðf 0:9069Þ2 ðmþ1Þ þ 1 ðf 0:9069Þðmþ1Þ ð63Þ

The technique of iterated roots could be extend to the case when the amplitude in the high-concentration regime is not known in advance, while the Padé approximants can ‘‘only” interpolate. In Fig. 4, iterated root and Padé approximants are compared with asymptotic formulae. Note that formula (63) anticipates an expansion in powers of ðf c  f Þ

1=2

in the vicinity of f c ,

lðf ; mÞ ’ AðmÞðf c  f Þ3=2 ð1 þ

A1 ðmÞ A2 ðmÞ 1=2 ðf  f Þ þ ðf  f Þ þ   Þ; AðmÞ c AðmÞ c

while iterated roots give the expansion in powers of ðf c  f Þ. All our formulae better reproduce the high-concentration limit, and thus expected to work well for large and moderately large concentrations, f > 0:6. In order to improve accuracy for smaller concentrations we develop below an approximants for lower and upper bounds and take their weighted average.

Mind that z ¼ zðf Þ ¼ f cff . Formula (62) does respect the lowconcentration and high-concentration limit-cases as shown in Fig. 4. Results for shear modulus for all concentrations and various Poisson’s coefficients, obtained with formulae (62) are shown in Fig. 2. When m ! 1, the material is referred as incompressible. When m ! 1=2, the material is highly compressible. The Poisson ratio is located in the interval 1=2 6 m 6 1 for plane strain. The most highly compressible materials such as critical fluids are characterized by m ! 1=2 [28]. Thus the dependencies demonstrating a decreasing effective modulus with decreasing Poisson’s ratio as shown in Fig. 2. are in agreement with intuition. In Fig. 3 the effective shear moduli of holes embedded into the matrix with m ¼ 1, are compared for different formulae. It is important to note that (62) adheres to both Hashin-Shtrikman bounds, while asymptotic expression with amplitude (59) only adheres to one of them.

4.5. Formula for the typical material cases For ceramics, metals and polymer matrices, the 2D Poisson ratio is typically in the range 0:33 < m < 0:54. The two relevant material cases studied in [34] belong to this interval and can be conveniently used for comparison. Upper bound is given by the simplest non-trivial additive selfsimilar approximant which respects the high-concentration limit and is correct as f ¼ 0,

lu ðf ; mÞ ¼ AðmÞðf c  f Þ3=2 þ ðf c  f Þ2 y2 ðmÞ; 3=2 1  f c AðmÞ y2 ðmÞ ¼ ; 2 fc

ð64Þ

while the lower bound is given by the next-order additive approximant, which also respects the first non-trivial term in the lowconcentration limit,

Fig. 2. Holes. Results for shear modulus for all concentrations and various Poisson’s coefficients, obtained with formula (62). The case of m ¼ 1 is shown with dotted line, and compared with m ¼ 0:5 (dashed), m ¼ 0 (dot-dashed), and m ¼  12 (solid).

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

403

Fig. 3. Case of m ¼  12. log ðl Þ of the effective shear modulus l ðf Þ for holes embedded into the matrix. The Hashin-Shtrikman bounds are shown with solid lines. The results from Eq. (62) are shown with dashed line. Eq. (59) is shown with dot-dashed line.

Fig. 4. Holes. Results for shear modulus for all concentrations and m ¼ 1, obtained with formula (62) are shown with dotted line; and corresponding results obtained with formula (63) are shown with dot-dashed line. They are compared with the original series (dashed) and with the high-concentration formula (59), shown with solid line.

ll ðf ; mÞ ¼ AðmÞðf c  f Þ3=2 þ u1 ðmÞðf c  f Þ2 þ u2 ðmÞðf c  f Þ5=2 ; 3=2 2f c a1 ðmÞ þ 2f c AðmÞ  5 ; u1 ðmÞ ¼  2 fc

u2 ðmÞ ¼ 

2f c a1 ðmÞ þ

3=2 f c ðAð 5=2 fc

mÞÞ þ 4

ð65Þ

:

When the subsequent self-similar approximants lk ðxÞ display substantial scatter it is proved effective to introduce the average form k X lk ðf Þ ¼ pi ðf Þli ðf Þ;

ð66Þ

i¼1

where the approximants

li ðf Þ are weighted with probabilities

jMi ðf Þj1 pi ðf Þ ¼ PN ; 1 j¼1 jMj ðf Þj

ð67Þ

in which N is the number of approximants involved and

Mi ðf Þ ¼

@ @f @ @f

li ðf Þ ; l0 ðf Þ

ð68Þ

are the so-called mapping multipliers [35,36]. The average is composed in such a way that it does not depend on the ‘‘initial condition” l0 ðf Þ. In our case the weighted average (66) of the lower and upper bounds, is given by 3=2 2 5=2 2ðAð mÞðf c f Þ þðf f c Þ u1 ðmÞþðf c f Þ u2 ðm ÞÞ

mÞðf c f Þ þðf f c Þ y2 ðmÞ p ffiffiffiffiffiffiffi þ Að 4ðf f c Þu1 ðmÞþpffiffiffiffiffiffiffi 3 f c f ð5ðf f c Þu2 ðmÞ3AðmÞÞ f c f AðmÞþ2ðf c f Þy2 ðmÞ

l ðf ; mÞ ¼

4ðf f c Þu1 ðmÞþ

3=2

2

2

þ pffiffiffiffiffiffiffi 1 3 f c f AðmÞþ2ðf c f Þy2 ðmÞ f c f ð5ðf f c Þu2 ðmÞ3AðmÞÞ

pffiffiffiffiffiffiffi2

:

2

ð69Þ Thus formula (69) provides the two-parametric formula for the effective shear modulus. The bounds and corresponding averages are shown in Fig. 5 for two typical values of the Poisson ratio. Generally speaking, the average formula is applicable for all m > 0:25. Neither of the bounds is good for all concentrations but their average is in much better overall agreement with the numerical results. For the lower bound we can also use the simplest root approximant which respects the first non-trivial term in the lowconcentration limit and general form of the critical behavior,

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

404

Fig. 5. Top: Upper bound (64) (dashed) compared to lower bound (65) (dotted) and weighted average (69) (solid) for m ¼ 0:33. Numerical results from [34] are shown as well. Bottom: Upper bound (64) (dashed) compared to lower bound (65) (dotted) and weighted average (69) (solid) for m ¼ 0:5. Numerical results from [34] are shown as well.

lr ðf ; mÞ ¼ 

1

ð70Þ

approximation requires 130 s. Every next approximation goes about 6 times longer than the previous one. Next, numerical computations are performed in order to determine geometrical param-

and it is straightforward to calculate the weighted average of

eters for the hexagonal array. The value of e3 is given explicitly by (22). Other similar sums are computed instantly for the regular hexagonal array. Computation of the sums for general locations of fibers are expensive and will be presented in a separate paper. The final formulae with the elastic constants and geometrical parameters in symbolic form are accessible as huge expressions. Their visualization in computer requires too much time. The retaining information of the local fields is too expensive. The package Mathematica contains analytical expression for local elastic fields and the effective constants in machine codes. We analyze these long formulae by substitution of numerical values, preparing of graphs and using other standard manipulations like integration. So, we can work with the obtained formulae knowing that they exist but we cannot look at them except their short fragments. In the present paper, we narrow the number of symbolic parameters and consider the special hexagonal geometry to present observable formulae for the effective constants. The mechanical properties of the inclusions are taken to be j1 ¼ 1 and l1 ¼ 0 in computations. Though the series in concentration is obtained for general elastic constants of fibers l1 and m1 the numerical computations for l1 > l show that the effective shear modulus in the form of series for the hexagonal array does not satisfy the Hashin-Shtrikman bounds for some f. Preliminary we can say that this case requires higher terms in f what will be discuss in a separate paper. It is

1

2:4184f ðf 0:9069Þðmþ1Þ

3=2 ;

Fig. 6. Top: Upper bound (64) (solid) compared to lower bound (70) (solid) and geometric average (71) (dotted) for m ¼ 0:33. Numerical results from [34] are shown as well. Bottom: Upper bound (64) (solid) compared to lower bound (70) (dotted) and geometric average (71) (dotted) for m ¼ 0:5. Numerical results from [34] are shown as well.

lu ðf ; mÞ and lr ðf ; mÞ. Such average is slightly accurate than (69) for

the two material cases from [34] as shown in Fig. 6, and is applicable for all m, but is cumbersome. Their geometric average,

la ðf ; mÞ ¼ ðlu ðf ; mÞlr ðf ; mÞÞ1=2 ; gives close values and is much simpler,

l

 a ðf ;

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u0:28372ðm þ 1Þð0:9069  f Þ3=2  0:297927ðf  0:9069Þ2 ðm  3:08105Þ mÞ ¼ u :  3=2 t 2:4184f 1  ðf 0:9069Þð mþ1Þ ð71Þ

Formula (71) is applicable for all m. For larger positive m it becomes close to (69). In the limiting case of rubber with m ¼ 1, (71), it is practically indistinguishable from (69). Same is true, e.g., for thallium with m ¼ 0:8182. For resin, with m ¼ 0:6129, the geometric average gives slightly higher results than weighted average. 5. Conclusion Symbolic computations of the effective constants are performed with the elastic constants j ; l j1 ; l1 and with geometrical parameters r; n; ak in the symbolic form. The obtained formulae for the local fields and the effective constants contains the geometrical parameters in the form of sums like (21). Computations up to 10

Oðf Þ requires about 3 h with a standard laptop. The eighth

ð1Þ

8

worth noting that expansions for random media up to Oðf Þ satisfy the Hashin-Shtrikman bounds [11]. The structure of series and

P. Drygas´ et al. / Computational Materials Science 139 (2017) 395–405

numerical examples for l1 < l show that high order approximations lie within the Hashin-Shtrikman bounds. This is the main reason why we restrict ourselves now by empty fibers when convergence to reasonable values is observed in resulting series. Formula (69) applies for the most typical cases of the matrices, while formula (71) can be applied for arbitrary elastic matrix. The final formula can be applied to the typical matrices made from ceramics, metals and polymers. It can be applied to other materials as well, e.g., for resin and thallium. Overall reasonable quality of the formulae can be improved at the expense of introducing terms of higher order obtained from the functional Eqs. (8) and (9). Acknowledgments PD, VM and WN were partially supported by Grant of the National Centre for Research and Development 2016/21/B/ ST8/01181. References [1] V.V. Vasiliev, E.V. Morozov, Advanced Mechanics of Composite Materials, Elsevier, Amsterdam etc., 2007. [2] Y. Fung, Biomechanics: Mechanical Properties of Living Tissues, Springer Science & Business Media, New York, 1993. [3] P.M. Adler, Porous Media: Geometry and Transports, Butterworth/Heinemann, Stoneham, 1992. [4] P.M. Adler, J.-F. Thovert, V.V. Mouvzenko, Fractured Porous Media, Oxford Univesity Press, Oxford, 2012. [5] B.D. Agarwal, L.J. Broutman, K. Chandrashekhara, Analysis and Performance of Fiber Composites, Wiley, NY, 2006. [6] N.I. Muskhelishvili, Some Basic Problems in the Mathematical Theory of Elasticity, Noordhoff, Groningen, 1954. [7] E.I. Grigolyuk, L.A. Filshtinsky, Perforated Plates and Shells, Nauka, Moscow, 1970 (in Russian). [8] E.I. Grigolyuk, L.A. Filshtinsky, Periodical Piece–Homogeneous Elastic Structures, Nauka, Moscow, 1991 (in Russian). [9] E.I. Grigolyuk, L.A. Filshtinsky, Regular Piece-Homogeneous Structures with Defects, Fiziko-Matematicheskaja Literatura, Moscow, 1994 (in Russian). [10] J. Helsing, An integral equation method for elastostatics of periodic composites, J. Mech. Phys. Solids 6 (1995) 815–828. [11] P. Drygas´, V. Mityushev, Effective elastic properties of random twodimensional composites, Int. J. Solids Struct. 97–98 (2016) 543–553. [12] I.V. Andrianov, L.I. Manevitch, Asymptotology: Ideas, Methods and Applications, Kluwer Academic Publishers, Dordrecht, Boston, London, 2002. [13] B. Budiansky, On the elastic moduli of some heterogeneous material, J. Mech. Phys. Solids 13 (1965) 223–227.

405

[14] T.S. Chow, J.J. Hermans, The elastic constants of fiber reinforced materials, J. Compos. Mater. 3 (1969) 382–396. [15] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer Verlag, New York, 2002. [16] V. Mityushev, S.V. Rogosin, Constructive Methods for Linear and Nonlinear Boundary Value Problems for Analytic Functions. Theory and Applications, Chapman & Hall CRC, London, etc., 1999. [17] V. Mityushev, Random 2D Composites and the Generalized Method of Schwarz, Adv. Math. Phys. 2015 (2015) 535128. [18] G.W. Milton, The Theory of Composites, Cambridge University Press, 2002. [19] P.G. De Gennes, On a relation between percolation theory and the elasticity of gels, Le J. Phys. Lett. 37 (1976) L1–L2. [20] A.R. Day, A. Snyder, E.J. Garboczi, M.F. Thorpe, The elastic moduli of a sheet containing circular holes, J. Mech. Phys. Solids 40 (5) (1992) 1031–1051. [21] A.V. Cherkaev, K.A. Lurie, G.W. Milton, Invariant properties of the stress in plane elasticity and equivalence classes of composites, Proc. Roy. Soc. A Lond. 438 (1992) 519–529. [22] M.F. Thorpe, I. Jasiuk, New results in the theory of elasticity for twodimensional composites, Proc. Roy. Soc. A Lond. 438 (1992) 531–544. [23] V. Mityushev, Thermoelastic plane problem for material with circular inclusions, Arch. Mech. 52 (6) (2010) 915–932. [24] P. Drygas´, Functional–differential equations in a class of analytic functions and its application to elastic composites, Complex Variables Elliptic Equat. 61 (8) (2016) 1145–1156. [25] S. Yakubovich, P. Drygas´, V. Mityushev, Closed-form evaluation of 2D static lattice sums, Proc. Roy. Soc. A 472 (2016) 20160510. [26] P. Drygas´, Generalized Eisenstein functions, J. Math. Anal. Appl. 444 (2) (2016) 1321–1331. [27] R. Lakes, Foam structures with a negative Poisson’s ratio, Science 235 (1987) 1038–1040. [28] G.N. Greaves, A.L. Greer, R.S. Lakes, T. Rouxe, Poisson’s ratio and modern materials, Nat. Mater. 10 (2011) 823–837. [29] J.W. Jiang, H.S. Park, Negative Poisson’s ratio in single-layer black phosphorus, Nat. Commun. 5 (2014) 47271. [30] S. Gluzman, V. Mityushev, W. Nawalaniec, Cross-properties of the effective conductivity of the regular array of ideal conductors, Arch. Mech. 66 (2014) 287–301. [31] S. Gluzman, V. Mityushev, Series, index and threshold for random 2D composite, Arch. Mech. 67 (2015) 75–93. [32] S. Gluzman, V. Mityushev, W. Nawalaniec, G.A. Starushenko, Effective conductivity and critical properties of a hexagonal array of superconducting cylinders, in: P.M. Pardalos, T.M. Rassias (Eds.), Contributions in Mathematics and Engineering: In Honor of Constantin Caratheodory, Springer, 2016, pp. 255–297. [33] S. Gluzman, V.I. Yukalov, D. Sornette, Self-similar factor approximants, Phys. Rev. E 67 (2003) 026109. [34] J.W. Eichen, S. Torquato, Determining elastic behavior of composites by the boundary element method, J. Appl. Phys. 74 (1993) 159–170. [35] V.I. Yukalov, S. Gluzman, Weighted fixed points in self-similar analysis of time series, Int. J. Mod. Phys. B 13 (1999) 1463–1476. [36] V.I. Yukalov, S. Gluzman, Extrapolation of power series by self-similar factor and root approximants, Int. J. Mod. Phys. B 18 (2004) 3027–3046.