The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite

The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite

Journal Pre-proof

The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite Fubin Tu, Yuyong Jiao, Xiaoyong Zhou, Yi Cheng, Fei Tan PII: DOI: Reference:

S0020-7683(19)30490-1 https://doi.org/10.1016/j.ijsolstr.2019.12.006 SAS 10557

To appear in:

International Journal of Solids and Structures

Received date: Revised date: Accepted date:

6 November 2019 3 December 2019 6 December 2019

Please cite this article as: Fubin Tu, Yuyong Jiao, Xiaoyong Zhou, Yi Cheng, Fei Tan, The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite, International Journal of Solids and Structures (2019), doi: https://doi.org/10.1016/j.ijsolstr.2019.12.006

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier Ltd.

Highlights • B-splines were implemented to discretize the polarized field in the Hashin and Shtrikman variational principle. • Only when the polarized field is sufficiently smooth, high order B-spline improves the numerical accuracy and efficiency. • Computational accuracy is improved if the constitutive matrix of a composite voxel is calculated by the laminate mixing rule. • The summation to calculate the discrete Green operator for periodic boundary conditions converges faster when higher order B-spline is used.

1

The implementation of B-splines to Hashin and Shtrikman variational principle based FFT method for the homogenization of composite Fubin Tua , Yuyong Jiaoa,∗, Xiaoyong Zhoua , Yi Chenga , Fei Tana a Faculty

of Engineering, China University of Geosciences (Wuhan), No.388 Lumo Road, Wuhan, Hubei, 430074, China

Abstract Fast Fourier transform (FFT) has been successfully used to estimate the effective properties of composites with meso periodic structure for more than two decades. Numerous improvements have been done to make it adequate for infinite contrast, bound elastic energy and composite voxels. The Hashin and Shtrikman variational principle based FFT method handles these three problems simultaneously and is very attractive. In this paper, B-splines were adopted to discretize the polarized field. The approximate polarized field was substituted into the Hashin and Shtrikman variational principle to obtain the discretized system. The constitutive matrix of a composite voxel was calculated by both the energetically consistent way and the laminate mixing rule. Numerical example of a square matrix with a circular inclusion shows that the results are more accurate than the remaining cases when the first and second order B-splines, as well as the laminate mixing constitutive matrix, are used. For higher order B-spline, the summation to compute the discrete Green operator for periodic boundary conditions converges faster while the discretized polarized field needs more iterations to converge. Only when the polarized field is sufficiently smooth, high order B-spline improves the numerical accuracy and efficiency. Keywords: B-splines, Discrete Green operator for periodic boundary ∗ Corresponding

author Email address: [email protected] (Yuyong Jiao)

Preprint submitted to International Journal of Solids and Structures

December 6, 2019

conditions, Composite voxel, Fast Fourier transform, Composite homogenization

1. Introduction Fast Fourier transform has been widely used to solve the mechanical, thermal, electromagnetic and hydraulic responses for heterogeneous medium with meso periodic structure and further estimate the effective properties of compos5

ites thanks to its merits of high computational efficiency, natural applicability to incompressible medium and regularly arranged voxels [1][2][3][4][5][6][7][8]. Brisard et al. [9][10][11] introduced a Hashin and Shtrikman variational principle and Galerkin approximation based FFT method which has the advantages of a rigorous elastic energy bound, discrete Green operator for periodic bound-

10

ary conditions with infinite frequencies, convergence for any phase contrast and an energetically consistent way for treating composite voxels. Nevertheless, the authors only used the piece-wise constant function for the discretization of polarized field. Other kinds of interpolation functions such as B-splines for discretizing the polarized field in Hashin and Shtrikman variational principle are

15

of great interest and remain to be implemented. It is well known that the minimum potential energy principle provides an elastic energy bound for finite element method (FEM) [12]. In an analogous way, the Hashin and Shtrikman variational principle [13] gives a rigorous bound on the elastic energy for heterogeneous problems [9][10][5]. Therefore, the FFT

20

method based on the Hashin and Shtrikman variational principle has a positive definite stiffness matrix for the discretized polarized field if the reference medium is stiffer than the stiffest phase. Since the appearance of FFT method, most of the researchers accounted for finite frequencies to compute the Green operator [1]. To alleviate the dis-

25

cretizing errors, Fourier expansion of the interpolation function with infinite frequencies and discrete Fourier transform of the polarized field at the center of voxels were combined to obtain a more accurate discrete Green operator 3

[9][10][11][5][14]. Brisard and Dormieux [10] named this kind of discrete Green operator as consistent discrete Green operator, which is effective to overcome 30

the Gibbs phenomenon[15]. Based on the consistent discrete Green operator, the polarized field is always convergent if the stiffest phase is chosen as the reference medium, even for infinite contrast, no extra treatment is required. In FEM, the FE mesh should be accommodative to the interface boundary to increase the estimated accuracy of the effective properties[16]. In FFT, due

35

to the regular arrangement of voxels, it is inevitable to meet a voxel with more than two material, which is called composite voxel. In the basic scheme, the constitutive matrix of the material at the center of a composite voxel is adopted to represent that of the composite voxel. However, this treatment may lead to unacceptable errors [9][17]. The Hashin and Shtrikman variational principle

40

and Galerkin approximation based FFT method naturally derives an equivalent and energetically consistent constitutive matrix for a composite voxel. Some other authors used a mixing rule named laminate mixing rule from the theory of composites [18] to construct an effective constitutive matrix for the composite voxel[17][19][20]. They took the normal of the material interface into account

45

and greatly reduced the error of stress or strain field near the material interface. B-splines and their derivants, such as non-uniformed rational B-splines (NURBS) and T-splines, are more attractive than Lagrange polynomials for Galerkin approximation where geometric accuracy and more than C 0 continuity are required [21][22]. Up to now, no attempt has been made to employ them to discretize the

50

polarized field in FFT method based on the Hashin and Shtrikman variational principle. Furthermore, for a heterogeneous medium with smoothly varying Young’s modulus, the polarized field is smooth everywhere. Thus, it makes sense to try to use high order B-splines to discretize the polarized field. Brisard et al. [9][10][11][5] only used the zero order B-spline (piece-wise constant

55

function) to discretize the polarized field. The absolute value of the Fourier expansion coefficient of the piece-wise constant function decreases slowly as the frequency grows, which requires more terms summation when calculating the converged discrete Green operator for periodic boundary conditions. For higher 4

order B-splines, the decreasing speed is faster since they are the convolutions of 60

piece-wise constant functions. In this paper, we devoted to implement B-splines to discretize the polarized field in the Hashin and Shtrikman variational principle and compare the results for different orders B-splines. We do not expect discretizations based on high order B-splines to improve on the accuracy and efficiency, unless the solution is expected to be sufficiently smooth. The importance

65

of considering the interface normal in a composite voxel through calculating its constitutive matrix by the laminate mixing rule was investigated. This paper is organized as follows. In the following section, backgrounds including the B-splines and the Lippmann-Schwinger equation are recalled. Then, the Hashin and Shtrikman variational principle and its B-spline discretization,

70

the algorithm for solving the discretized system and convergence detection are presented. In the fourth section, the discrete Green operators for periodic boundary conditions for different orders B-splines are discussed in detail. The laminate mixing rule for calculating the constitutive matrix of a composite voxel is given in the fifth section. Some numerical examples are performed to investi-

75

gate the effectiveness of the presented B-spline discretized Hashin and Shtrikman variational principle based FFT method in the sixth section. A short conclusion is given in the final section.

2. Backgrounds 2.1. B-splines B-splines can be defined by their convolutional or recursive property [23]. Using the latter property, we start with the piece-wise constant function as per

Nb,0 (ξ) =

  1  0

if ξb ≤ ξ ≤ ξb+1

(1)

otherwise.

where ξb denotes the coordinate of the bth knot. Obviously, the piece-wise constant function is the 0 order B-spline. For orders of p = 1, 2, 3, · · · , B-splines

5

are defined by Nb,p (ξ) =

80

ξ − ξb ξb+p+1 − ξ Nb,p−1 (ξ) + Nb+1,p−1 (ξ) ξb+p − ξb ξb+p+1 − ξb+1

(2)

Figure 1 shows the p = 0, 1, 2, 3 orders B-splines. They are positive everywhere and satisfy partition of unity [24][25]. The B-spline functions are of C p−1 continuity. They are compactly supported. The width of the B-spline’s support increases as the order grows. The 2D and 3D B-spline functions can be easily constructed by 2 and 3 times products of 1D B-splines. N1,0 N1,1 N1,2 N1,3

1

N1,p

0.8

0.6

0.4

0.2

0 0

1

2

3

4

ξ

Fig. 1. The p = 0, 1, 2, 3 orders B-splines.

The Fourier transform of B-splines is ˆp (ω) = sincp+1 ( ω ) N 2 85

where ω is the angular frequency, sinc stands for the function sincx =

(3) sinx x .

2.2. The Lippmann-Schwinger equation Let x be the coordinates of a point in a meso heterogeneous material with periodic structure, the constitutive matrix is a function of the coordinates and

6

denoted as C(x). To predict the effective properties of the heterogeneous material, the medium is submitted to a macroscopic strain E, which leads to displacements fluctuation u periodically fluctuates around the macroscopic counterpart E · x. The governing equations are   ∇ · σ = 0    σ = C(x) :       = E + 1 [∇ ⊗ u + (∇ ⊗ u)T ] 2

(4)

where ∇ is the Hamilton operator, ⊗ is the dyad operator, σ and  stand for the stress and strain. If the domain of a unit cell is Ω = {x ∈ Rd |0 ≤ xj ≤ Lj , j = 1, · · · d}, d = 2, 3

(5)

in which d stands for the dimension of the problem, then the fluctuation of the displacements satisfies u(x1 , · · · , xj + Lj , · · · , xd ) = u(x)

(6)

For the convenience of finding the solution, the heterogeneous problem governed by Equation 4 is usually transformed to a homogeneous problem by introducing a so called reference medium with constant constitutive matrix C0 in the whole cell. Then the mesoscopic constitutive equation becomes σ(x) = C0 : (x) + [C(x) − C0 ] : (x) = C0 : (x) + φ(x)

(7)

where φ(x) = [C(x) − C0 ] : (x) is defined as the polarized field. The transformed homogeneous problem can be solved by Green function method as (x) = E − (Γ0 ∗ φ)(x)

(8)

where Γ0 is the fourth rank continuous Green operator for periodic boundary

7

conditions for strains relies on C0 [26], ∗ denotes the convolution product. This equation can be rewritten by replacing the φ with  as (x) + {Γ0 ∗ [(C − C0 ) : ]}(x) = E

(9)

which is the well-known Lippmann-Schwinger equation. When the polarized field is taken as the unknown field, the Lippmann-Schwinger equation becomes [(C − C0 )−1 : φ](x) + (Γ0 ∗ φ)(x) = E

(10)

The convolution product in Ω space becomes arithmetic product in the Fourier space. For the sake of finding the solution, it’s convenient to write Equation 8 in the Fourier space as ˆ ˆ(ω) = −Γˆ0 (ω) : φ(ω)

(ω 6= 0),

ˆ(0) = E

(11)

where ω is the angular frequency in d dimensions, which is expressed as ω=(

2πα1 2παj 2παd T ,··· , ,··· , ) L1 Lj Ld

(12)

in which αj is the wave number in j direction. Detailed expressions of Γˆ0 can be referred to [27] for isotropic and anisotropic reference material. The continuous Green operator for isotropic case is as follows Γˆ0,jklm =

1 (δjl ωk ωm + δjm ωk ωl + δkl ωj ωm + δkm ωj ωl ) 4µ0 |ω|2 1 − ωj ωk ωl ωm 2µ0 (1 − ν0 )|ω|4

(13)

where µ0 and ν0 are the shear modulus and Poisson’s ratio of the reference material, respectively. µ0 can be expressed by the Young’s modulus E0 and Poisson’s ratio as µ0 =

E0 2(1 + ν0 )

8

(14)

3. The B-splines discretized Hashin and Shtrikman variational principle and algorithm for solving the discretized system 3.1. The Hashin and Shtrikman variational principle and its discretization by 90

B-splines Hashin and Shtrikman variational principle [13] was derived based on the principle of minimum potential energy, Legendre transformation and balance equation. Considering the fact that a variational principle is equivalent to a virtual displacement principle, we use the virtual polarized field δφ to multiply Equation 10 and then integrate the formula over the domain of a unit cell. Through this operation, we have Z

Ω

−1

δφ : (C − C0 )

: φdV +

Z

Ω

δφ : Γ0 ∗ φdV =

Z

δφ : EdV

(15)

Ω

Willis [28] has proved that the quadratic form of polarized field on the left side of this formula is positive definite when C0 ≥ C(x). This property leads to a discretized algebra system with a positive definite stiffness matrix. B-splines with orders of {0,1,2,3} are used to discretize the polarized field φ(x). The polarized field in a voxel va is discretized as

φva =

nva X

Nb (ξ)φb

(16)

b=1

where Nb (ξ) is the 2D or 3D B-spline interpolation function for control point b, φb is the polarized field on control point b, nva is the total number of control points whose interpolation functions’ supports have interaction with the va voxel. The virtual polarized field is discretized in the same way. Writing the discretized polarized field in a matrix form, we have   φ = N Φ

 δφ = N δΦ 9

(17)

where N denotes the shape function matrix, and Φ stands for the polarized field vector on all of the control points. Substituting Equation 17 into the first term of Equation 15, we obtain Z

Z δφ : (C − C0 )−1 : φdV = δΦT N T (C − C0 )−1 N ΦdV Ω Ω Z T = δΦ N T (C − C0 )−1 N dV Φ = δΦT K1 Φ

(18)

Ω

where K1 denotes the first part of the stiffness matrix. When integrating 95

the stiffness matrix of a composite voxel, a subdomain-wise Gauss integration method [29] was used. The composite voxel was first divided into several subvoxels. Then Gauss integration method was adopted for each sub-voxel and C(x) at every Gauss points could be taken into account. More details about the subdomain Gauss integration method can be referred to [30][31][32]. Regarding the term on the right side of Equation 15, we have Z

Z δφ : EdV = δΦT N T EdV Ω Ω Z T = δΦ N T EdV = δΦT F

(19)

Ω

100

where F denotes the macro strain vector that is exerted on the control points. The second term in Equation 15 is discretized as Z

Ω

{

∞ X

m1 =−∞

···

δφ : Γ0 ∗ φdV =

∞ X

sinc2(p+1)

md =−∞

M M 1 −1 d −1 X X L1 · · · Ld ˆ∗ · · · δφ α1 ,··· ,αd : (M1 · · · Md )2 α =0 α =0 1

d

π(α1 + m1 M1 ) π(αd + md Md ) (20) · · · sinc2(p+1) M1 Md

2π(α1 + m1 M1 ) 2π(αd + md Md ) ˆα ,··· ,α Γˆ0 [ ,··· , ]} : φ 1 d L1 Ld

The detailed derivation process can be referred to Appendix A. Then, we have

10

the discrete Green operator for periodic boundary conditions as Γˆ0;α1 ,··· ,αd =

∞ X

m1 =−∞

···

∞ X

sinc2(p+1)

md =−∞

π(α1 + m1 M1 ) ··· M1

(21)

π(αd + md Md ) ˆ 2π(α1 + m1 M1 ) 2π(αd + md Md ) sinc2(p+1) Γ0 [ ,··· , ] Md L1 Ld If we denote θˆα1 ,··· ,αd as verse discrete Fourier

L1 ···Ld ˆ ˆ (M1 ···Md )2 Γ0;α1 ,··· ,αd : φα1 ,··· ,αd transform θa1 ,··· ,ad = DFT−1 θˆα1 ,··· ,αd ,

and use the inthe second term

in Equation 15 becomes Z

Ω

δφ : Γ0 ∗ φdV = =

M 1 −1 X a1 =0

···

M 1 −1 X α1 =0

M d −1 X

···

M d −1 X

ˆ∗ ˆ δφ α1 ,··· ,αd : θα1 ,··· ,αd

αd =0

(22)

δφa1 ,··· ,ad : θa1 ,··· ,ad = δΦT Θ

ad =0

Finally, considering the arbitrariness of δΦ, Equation 15 can be rewritten as K1 Φ + Θ = F

(23)

The basic scheme was introduced by Moulinec and Suquet [33]. It was shown by Brisard and Dormieux [9] that the basic scheme could be seen as an asymptotically consistent approximation of the variational problem 18 (for p = 0) where Equation 21 is replaced with the following Equation 2πα1 2παd Γˆ0;α1 ,··· ,αd = Γˆ0 ( ,··· , ) L1 Ld

(24)

3.2. Algorithm for solving the discretized system and convergence detection Equation 23 is a positive definite linear algebra system with the unknowns Φ for linear elastic medium. It’s very difficult to have the convolution term Θ explicitly. As a consequence, a matrix-free and iterative solver was used. In this 105

paper, the conjugate gradient (CG) method was adopted to solve Equation 23. K1 Φ + Θ was carried out in the following four steps

11

ˆα ,··· ,α using the fast Fourier transform of φa ,··· ,a ; (1) compute the φ 1 1 d d ˆα ,··· ,α in the Fourier space; (2) compute θˆα1 ,··· ,αd = Γˆ0;α1 ,··· ,αd : φ 1 d (3) compute the θa1 ,··· ,ad using the inverse fast Fourier transform; 110

(4) compute K1 Φ + Θ in real space. Detailed algorithm associated with the conjugate gradient method can be referred to [34]. The convergence of Φ is detected by the following formula (Φn+1 − Φn ) : (Φn+1 − Φn ) ≤ crit Φn+1 : Φn+1

(25)

where Φn and Φn+1 are the polarized field vector at n and n + 1 time steps, crit is the critical error when Φ converges.

115

4. On the discrete Green operator for periodic boundary conditions Equation 21 accounts for infinite frequencies and Fourier expansion of the B-splines is employed herein. To rewrite the discrete Green operator expressed in Equation 21 in a more convenient form, we first simplify the following formula sinc2(p+1)

π(αj + mj Mj ) αj παj =( )2(p+1) sinc2(p+1) Mj αj + mj Mj Mj

(26)

Substituting Equation 26 into Equation 21, we obtain πα1 παd 2(p+1) 2(p+1) Γˆ0;α1 ,··· ,αd = α1 sinc2(p+1) · · · αd sinc2(p+1) M1 Md ∞ ∞ X X −2(p+1) −2(p+1) ··· (α1 + m1 M1 ) · · · (αd + md Md )

m1 =−∞

md =−∞

2π(α1 + m1 M1 ) 2π(αd + md Md ) Γˆ0 [ ,··· , ] L1 Ld

12

(27)

If all of αj = 0 except that αk 6= 0, then παk 2(p+1) Γˆ0;0,··· ,αk ,··· ,0 = αk sinc2(p+1) Mk

∞ X

mk =−∞

(28)

2π(αk + mk Mk ) −2(p+1) ˆ (αk + mk Mk ) Γ0 [0, · · · , , · · · , 0] Lk When all of the αj = 0, we have Γˆ0;0,··· ,0 = 0

(29)

Equations 27 and 28 are summations from −∞ to ∞, it has been proved that the convergence of the discrete Green operator for periodic boundary conditions 120

is very slow if we truncate the summation at a finite integer [9][14]. Brisard and Dormieux [9] used Riemann sums [35] to estimate the remainder terms when truncating the summation at a finite integer. However, in the literature [14], we find that only the first several terms of the consistent discrete Green operator have influence on the results. Therefore, in our paper, we truncate the

125

summation at a finite integer for calculating the discrete Green operator. It is worth mentioning that the consistent discrete Green operator in [14] is only conditionally convergent and further theoretical justification should be offered when using it. Figure 2 draws the sinc2(p+1) ( ω2 ) with p = {0, 1, 2, 3}. Figure 2 shows that

130

the frequency band shrinks narrower and the amplitude decays faster as p grows. In other words, the convergence of the summation to compute the discrete Green operator may be faster and it is possible to truncate the summation at a smaller integer for larger p.

5. On the composite voxel 135

Equation 18 is derived from Hashin and Shtrikman variational principle so that it’s named as energetically consistent stiffness. In this formula, the polarized field is assumed to be a constant in a composite voxel, which conflicts 13

p=0 p=1 p=2 p=3

1

sinc2(p+1) ( ω2 )

0.8

0.6

0.4

0.2

0 −16

−12

−8

−4

0

4

8

12

16

ω

Fig. 2. The sinc2(p+1) ( ω ) with p = {0, 1, 2, 3}. 2

with the real case. To overcome this shortcoming, Kabel et al. [17] introduced a superior laminate mixing rule that accounts for the normal of the material 140

interface. If we denote the elastic matrix of a composite voxel as Clam , it is obtained by solving the following system [17][18] −1 −1

(P + λ(Clam − λI)

)

=

R

Ωva

(P + λ(C − λI)−1 )−1 dV Vva

(30)

where I is the identity matrix, Ωva is the domain occupied by ath voxel, Vva is the volume of ath voxel, λ is a parameter which should be larger than the largest eigenvalue of the elastic matrix C for all x. P is a matrix relies on the material interface normal. The tensor form of P is expressed as Pjklm =

1 (nj δkl nm + nj δkm nl + nk δjl nm + nk δjm nl ) − nj nk nl nm 2

(31)

where nj is the normal of the material interface in j direction. The approximate method for computing nj can be referred to [17].

14

1m

r Ei υ i Em υm 1m

Fig. 3. Schematic diagram of a unit cell.

Replacing C in Equation 18 by Clam , we have K1 =

Z

Ω

N T (Clam − C0 )−1 N dV

(32)

6. Numerical example 6.1. Square matrix with a circular inclusion 145

Regarding the numerical example, a 2D (plane stress elasticity) periodic structure with a circular inclusion was investigated. A unit cell of the periodic structure is shown in Figure 3. The size of the cell was set to be 1m × 1m. The radius of the circular inclusion is r. r = 0.05m in the following 4 subsubsections and r is a variable relies on the fiber volume fraction in the final subsubsection.

150

Ei and Em are the Young’s moduli of the inclusion and matrix. Em was fixed to be 20GPa. νi and νm are the Poisson’s ratios of the inclusion and matrix, respectively. νi = νm = 0.3 was considered. Ex = Ey = 0.005 was imposed on the periodic structure. For Hashin and Shtrikman variational principle based scheme, the C0 of the

155

reference medium was set to 1.01 times of that of the stiffest phase to avoid sin15

gularity when calculating (C − C0 )−1 . While for the basic scheme, the Young’s modulus of the reference medium was E0 =

Ei +Em 2

according to Michel et al.

[1]. Square voxels were adopted, i.e., M1 = M2 for 2D problem. 8 × 8, 16 × 16, 32 × 32, · · · , 256 × 256 grids were used for FFT analyses. The critical error crit

160

to detect the convergence of the polarized field was set to be 1.0 × 10−8 . λ was chosen as the Young’s modulus of the stiffest phase.

6.1.1. Summation truncation for the discrete Green operator for periodic boundary conditions The Young’s modulus of the circular inclusion was set as 200GPa and 256 × 256 165

grids were used. Figure 2 shows that the Fourier transform of B-splines has the largest remainder term at ωj =

Mj 2

− 1. Therefore, Γˆ0;xxxx;127,127 was chosen

as a representative to investigate the reasonableness for truncating the summation at a finite integer to calculate the discrete Green operator. To analyze the errors, we assumed that the summation that was truncated at 10,000 could be 170

used to represent the infinite summation or as a reference value. The bounds of mj were -10,000 and 10,000 when truncating the summation at 10,000. To reduce the truncation errors from computer, we first ordered the terms by their magnitudes. Then, we performed the summations by order of increasing magnitude.

175

Figure 4 shows the relative error of Γˆ0;xxxx;127,127 versus the upper bound of mj for B-splines with p = {0, 1, 2, 3}. For all orders of B-splines, Γˆ0;xxxx;127,127 gradually converges to the summation truncated at 10,000 as the bound of mj grows. For p = 0, the convergence rate is very slow. The higher the B-spline order, the faster the convergence rate of Γˆ0;xxxx;127,127 . In the following sub-

180

subsection, we truncate the summation to calculate the discrete Green operator when the relative error of Γˆ0;xxxx;127,127 < 1.0 × 10−4 , the corresponding upper bounds of mj are {2165, 6, 2, 2} for p = {0, 1, 2, 3}.

16

Fig. 4. The relative error of Γˆ0;xxxx;127,127 versus the upper bound of mj for B-splines with p = {0, 1, 2, 3}.

6.1.2. Strain energy When investigating the strain energy, the Young’s modulus of the circular 185

inclusion was set as 200GPa. Finite element results were also provided for comparison. Figure 5 shows the strain energies versus the voxel numbers for FEM, basic scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme. For the basic scheme, the strain energy is not monotonic as the voxel number increases. The basic scheme does not provide a strain

190

energy bound. The strain energies monotonically approach to the FE results for all orders of B-splines as the voxel number increases. For voxel number less than 64, the results for Hashin and Shtrikman variational principle based schemes are deviatoric from the FE results more than the basic scheme. The errors may result from the assumption that the polarized field is constant in a

195

composite voxel. The strain energy for the higher order B-spline is a bit closer to the FE result than that for the lower order B-spline. To decrease the numerical errors for small voxel number, the elastic matrix for a composite voxel was calculated by the laminate mixing rule. Figure 6 shows the strain energies versus the voxel numbers for FEM, basic scheme and

17

7.8

FEM Basic scheme B-spline p = 0 B-spline p = 1 B-spline p = 2 B-spline p = 3

Strain energy (105 J)

7.7 7.6 7.5 7.4 7.3 7.2 7.1 4

8

16

32

64

128

256

512

Mj

Fig. 5. Strain energies versus the voxel numbers for FEM, basic scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme.

200

B-splines discretized Hashin and Shtrikman variational principle based scheme combined with constitutive matrix computed by the laminate mixing rule. The order of the B-spline has little influence on the result. For voxel number less than 64, the strain energy for laminate mixing constitutive matrix is closer to the FE result than that for energetically consistent constitutive matrix.

205

6.1.3. Number of iterations to convergence To investigate the number of iterations to convergence, we used 128 × 128 voxels while varied the phase contrast. The Young’s modulus of the inclusion was varied from 0.00002GPa to 20,000,000GPa. The remaining parameters were the same as those for the numerical example in the last subsubsection. Here we

210

only investigated the basic scheme, B-splines of 0 and 2 orders. Figure 7 shows the iteration numbers to convergence versus

Ei Em

for the basic

scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme. The basic scheme diverges as converges as 215

Ei Em

Ei Em

approaches to ∞. However, it

→ 0, the same conclusion as that in [36]. For the Hashin and

Shtrikman variational principle based scheme discretized with B-splines, the

18

7.8 B-spline B-spline B-spline B-spline

Strain energy (105 J)

7.7

FEM Basic scheme p = 0, Laminate p = 1, Laminate p = 2, Laminate p = 3, Laminate

7.6 7.5 7.4 7.3 7.2 7.1 4

8

16

32

64

128

256

512

Mj

Fig. 6. Strain energies versus the voxel numbers for FEM, basic scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme combined with constitutive matrix computed by the laminate mixing rule.

iteration number to convergence increases as

Ei Em

grows. The iteration number

to convergence approaches a bound value for infinite contrast. The iterations to reach convergence are more for higher order B-spline. Figure 8 shows the iteration numbers to convergence versus 220

Ei Em

for B-splines

discretized Hashin and Shtrikman variational principle based scheme without and with constitutive matrix computed by the laminate mixing rule. When the constitutive matrix is computed by the laminate mixing rule, it needs more iterations to convergence. For p = 0 and

Ei Em

< 1, the iteration numbers to conver-

gence are similar for the cases without and with constitutive matrix computed 225

by the laminate mixing rule. However, the trend is contrary for p = 2, the iteration numbers to convergence are similar for the cases without and with constitutive matrix computed by the laminate mixing rule as

Ei Em

> 1.

6.1.4. Strain fields To compare the strain fields for different cases, 256 × 256 voxels were used to 230

discretize the unit cell. 512 × 512 elements were used to obtain the FE result for comparison. The Young’s modulus of the circular inclusion was set as 200GPa.

19

Fig. 7. Iteration numbers to convergence versus EEi for the basic scheme and B-splines m discretized Hashin and Shtrikman variational principle based scheme.

Fig. 8. Iteration numbers to convergence versus EEi for B-splines discretized Hashin and m Shtrikman variational principle based scheme without and with constitutive matrix computed by the laminate mixing rule.

20

Figure 9 shows the strain fields x for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme. Due to the symmetry, only a quarter of the contour for the unit cell is presented. 235

Compared with the FE result, the basic scheme overestimates the maximum strain and underestimates the minimum strain near the inclusion, i.e., the white zone near the fiber exceeds the upper and lower bounds of the color bar. The results calculated by the B-splines discretized Hashin and Shtrikman variational principle based scheme are closer to the FE result than the basic scheme. The

240

maximum x is underestimated a bit for p = 0. As expected, no significant improvement is observed when high order B-splines were used to discretize the polarized field in Hashin and Shtrikman variational principle. The strain fields x for the B-splines discretized Hashin and Shtrikman variational principle based scheme with constitutive matrix computed by the lam-

245

inate mixing rule are shown in Figure 10. When the constitutive matrix is computed by the laminate mixing rule, the strain fields are closer to the FE result. The maximum x is still underestimated a bit for p = 0. 6.1.5. Effective bulk modulus The macro homogeneous properties of the composite with meso periodic structure can be obtained by FFT analysis. In the example we investigated, we can have the effective bulk modulus κ by taking the volume average of Equation 7 as follows R σx +σy R φx +φy ( 2 )dV ( 2 )dV E0 Ω κ= R + Ω = 2(1 − ν0 ) Ex + Ey ( + y )dV Ω x

(33)

The volume fraction of the fiber was varied while the remaining parameters were 250

the same as those for the example in the last subsubsection. Figure 11 shows the variation of the effective bulk modulus as the volume fraction of fiber changes for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme. The effective bulk modulus varies nonlinearly as the fiber volume fraction increases. The ba-

21

FEM

Basic scheme 0.009

0.6

0.009

0.6

0.008

0.008 0.58

0.007 0.006

0.56

0.005 0.004

0.54

y (m)

y (m)

0.58

0.007 0.006

0.56

0.005 0.004

0.54

0.003 0.52

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

0.58

0.002 0.001

0.5

0.6

0.5

0.52

0.54

0.56

x (m)

x (m)

(a)

(b)

B-spline p = 0

0.58

0.6

B-spline p = 1 0.009

0.6

0.009

0.6

0.008

0.008 0.58

0.007 0.006

0.56

0.005 0.004

0.54

y (m)

y (m)

0.58

0.007 0.006

0.56

0.005 0.004

0.54

0.003 0.52

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

0.58

0.002 0.001

0.5

0.6

0.5

0.52

0.54

0.56

x (m)

x (m)

(c)

(d)

B-spline p = 2

0.58

0.6

B-spline p = 3 0.009

0.6

0.009

0.6

0.008

0.008 0.58

0.007 0.006

0.56

0.005 0.004

0.54

y (m)

y (m)

0.58

0.007 0.006

0.56

0.005 0.004

0.54

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

0.58

0.6

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

x (m)

x (m)

(e)

(f )

0.58

0.6

Fig. 9. The strain fields x for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme: (a) FEM; (b) Basic scheme; (c) B-spline, p = 0; (d) B-spline, p = 1; (e) B-spline, p = 2; (f) B-spline, p = 3.

22

B-spline p = 0, Laminate

B-spline p = 1, Laminate 0.009

0.6

0.009

0.6

0.008

0.008 0.58

0.007 0.006

0.56

0.005 0.004

0.54

y (m)

y (m)

0.58

0.007 0.006

0.56

0.005 0.004

0.54

0.003 0.52

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

0.58

0.002 0.001

0.5

0.6

0.5

0.52

0.54

0.56

x (m)

x (m)

(a)

(b)

B-spline p = 2, Laminate

0.58

0.6

B-spline p = 3, Laminate

0.6

0.6

0.009

0.009

0.008

0.008 0.58

0.007 0.006

0.56

0.005 0.004

0.54

y (m)

y (m)

0.58

0.007 0.006

0.56

0.005 0.004

0.54

0.003 0.52

0.002 0.001

0.5 0.5

Fig. 10. principle B-spline, B-spline,

0.52

0.54

0.56

0.58

0.6

0.003 0.52

0.002 0.001

0.5 0.5

0.52

0.54

0.56

x (m)

x (m)

(c)

(d)

0.58

0.6

The strain fields x for the B-splines discretized Hashin and Shtrikman variational based scheme with constitutive matrix computed by the laminate mixing rule: (a) p = 0, Laminate; (b) B-spline, p = 1, Laminate; (c) B-spline, p = 2, Laminate; (d) p = 3, Laminate.

23

255

sic scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme predict the effective bulk modulus well compared with the FE result except the B-spline with order p = 3 for high fiber volume fraction.

Fig. 11. The variation of the effective bulk modulus as the volume fraction of fiber changes for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme.

Figure 12 shows the variation of the effective bulk modulus as the volume fraction of fiber changes for the FEM and the B-splines discretized Hashin and 260

Shtrikman variational principle based scheme with constitutive matrix computed by the laminate mixing rule. Even though the constitutive matrix is computed by the laminate mixing rule, B-spline with order p = 3 fails to estimate the effective bulk modulus when the fiber volume fraction is greater than 20%. The reason for this phenomenon is maybe because the support of the

265

third order B-spline is too wide. B-splines discretized Hashin and Shtrikman variational principle based scheme with constitutive matrix computed by the laminate mixing rule predicts the FE results well when the B-spline order is less than 3. 6.2. Square matrix with smoothly varying Young’s modulus The merit of the high order B-spline is that it is good at representing smoothly distributing field. To obtain smoothly distributing polarized field, the 24

Fig. 12. The variation of the effective bulk modulus as the volume fraction of fiber changes for the FEM and the B-splines discretized Hashin and Shtrikman variational principle based scheme with constitutive matrix computed by the laminate mixing rule.

material parameters of the meso structure should be smoothly varied. Therefore, a 2D (plane stress elasticity) periodic structure with smoothly varying Young’s modulus was investigated. The size of the unit cell was the same as that in the last subsection. The Poisson’s ratio of the whole matrix was set as ν = 0.3. Ex = Ey = 0.005 was imposed on the periodic structure. Emax and Emin are the maximum and minimum Young’s moduli of the matrix. Emin was fixed to be 20GPa. The Young’s modulus was varied as per E = (Emax − Emin ) 270

1 + sin(2πx) sin(2πy) + Emin 2

(34)

When Emax = 200GPa, Figure 13 shows the contour of Young’s modulus. In this example, 128 × 128 grid was used for FFT analysis. A FE simulation with 1024 × 1024 elements was carried out for comparison. The remaining parameters are the same as those in the last subsection. 6.2.1. Strain energies and strain fields

275

In this subsubsection, the maximum Young’s modulus was set to be 200GPa. Table 1 lists the strain energies for the FEM, the basic scheme and the B-splines

25

2 × 1011

1 0.9

1.8 × 1011

0.8

1.6 × 1011

y (m)

0.7

1.4 × 1011

0.6

1.2 × 1011

0.5

1 × 1011

0.4 0.3

8 × 1010

0.2

6 × 1010

0.1

4 × 1010

0

2 × 1010

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

x (m) Fig. 13. The contour of Young’s modulus when Emax = 200GPa.

26

Table 1 Strain energies for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme.

Item FEM Basic scheme B-spline, p = 0

Strain energy (106 J) 3.43284 3.4328 3.4462

Item B-spline, p = 1 B-spline, p = 2 B-spline, p = 3

Strain energy (106 J) 3.4336 3.4328 3.4335

discretized Hashin and Shtrikman variational principle based scheme. Except the Hashin and Shtrikman variational principle based scheme with 0 order Bspline, the strain energies for the basic scheme and B-splines with orders of 280

{1, 2, 3} are close to that for FEM. Figure 14 shows the strain fields x for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme. The Hashin and Shtrikman variational principle based scheme with 0 order B-spline overestimates the maximum strain and underestimates the minimum

285

strain (i.e., the white zone exceeds the upper and lower bounds of the color bar). The results obtained from the basic scheme and the Hashin and Shtrikman variational principle based scheme with B-splines of orders {1, 2, 3} are close to the FE result. To quantitatively analyze the errors of the strain fields, we define the relative error of strain x by the square of an approximate L2 norm as per

err =

127 P

127 P

2 FEM (FFT x;a1 ,a2 − x;a1 ,a2 ) a1 =0 a2 =0 127 P 127 P FEM 2 (x;a1 ,a2 ) a1 =0 a2 =0

(35)

FEM where FFT x;a1 ,a2 and x;a1 ,a2 are the strains obtained from FFT and FEM analy290

ses. Figure 15 shows the relative error of strain x for the B-splines discretized Hashin and Shtrikman variational principle based scheme versus the order of the B-spline. For heterogeneous medium with smoothly varying Young’s modulus, the relative error of strain x deceases as the order of B-spline increases. Considering the fact that only L2 regularity is required for the polarized field

27

FEM

Basic scheme

1

1

0.012

0.012

0.011

0.011 0.8

0.01 0.009

0.6

0.008 0.007

0.4

y (m)

y (m)

0.8

0.01 0.009

0.6

0.008 0.007

0.4

0.006

0.006

0.005

0.2

0.005

0.2

0.004

0.004

0.003

0 0

0.2

0.4

0.6

0.8

0.003

0

1

0

0.2

0.4

0.6

x (m)

x (m)

(a)

(b)

B-spline p = 0

0.8

1

B-spline p = 1

1

1

0.012

0.012

0.011

0.011 0.8

0.01 0.009

0.6

0.008 0.007

0.4

y (m)

y (m)

0.8

0.01 0.009

0.6

0.008 0.007

0.4

0.006

0.006

0.005

0.2

0.005

0.2

0.004

0.004

0.003

0 0

0.2

0.4

0.6

0.8

0.003

0

1

0

0.2

0.4

0.6

x (m)

x (m)

(c)

(d)

B-spline p = 2

0.8

1

B-spline p = 3

1

1

0.012

0.012

0.011

0.011 0.8

0.01 0.009

0.6

0.008 0.007

0.4

y (m)

y (m)

0.8

0.01 0.009

0.6

0.008 0.007

0.4

0.006 0.005

0.2

0.006 0.005

0.2

0.004 0.003

0 0

0.2

0.4

0.6

0.8

1

0.004 0.003

0 0

0.2

0.4

0.6

x (m)

x (m)

(e)

(f )

0.8

1

Fig. 14. The strain fields x for the FEM, the basic scheme and the B-splines discretized Hashin and Shtrikman variational principle based scheme: (a) FEM; (b) Basic scheme; (c) B-spline, p = 0; (d) B-spline, p = 1; (e) B-spline, p = 2; (f) B-spline, p = 3.

28

295

in Hashin and Shtrikman variational principle, the numerical accuracy has very little improvement for B-splines with p > 1 compared to the result for B-spline with p = 1.

Fig. 15. The relative error of strain x for the B-splines discretized Hashin and Shtrikman variational principle based scheme versus the order of the B-spline.

6.2.2. Number of iterations to convergence The maximum Young’s modulus Emax was varied from 20GPa to 20,000,000GPa. 300

The remaining parameters were the same as those for the numerical example in the last subsubsection. Here we only investigated the basic scheme, B-splines of 0 and 2 orders. Figure 16 shows the iteration numbers to convergence versus

Emax Emin

for the

basic scheme and B-splines discretized Hashin and Shtrikman variational princi305

ple based scheme. The iteration number to convergence approaches to a bound value for infinite contrast for the basic scheme and B-splines discretized Hashin and Shtrikman variational principle based scheme. Iterations to convergence for B-splines discretized Hashin and Shtrikman variational principle based scheme with 0 order B-spline are more than that for the basic scheme. For smoothly

310

varying Young’s modulus, the iterations to reach convergence for B-splines discretized Hashin and Shtrikman variational principle based scheme with 2 order 29

B-spline are less than that with 0 order B-spline, as well as the basic scheme when

Emax Emin

≥ 1000. High order B-spline discretized Hashin and Shtrikman vari-

ational principle based scheme improves the numerical accuracy and efficiency 315

a bit for the problem whose Young’s modulus is smoothly varied.

max Fig. 16. Iteration numbers to convergence versus E for the basic scheme and B-splines Emin discretized Hashin and Shtrikman variational principle based scheme.

7. Conclusion In this paper, different orders of B-splines were used to discretize the polarized field in FFT method for the homogenization of composites with meso periodic structure based on the Hashin and Shtrikman variational principle. The 320

constitutive matrix of a composite voxel was treated by both the energetically consistent method and the laminate mixing rule. Through numerical examples, we draw the following conclusions: (1) The higher the B-spline order, the faster the convergence of the discrete Green operator for periodic boundary conditions.

325

(2) The iteration number to reach convergence approaches to a finite value for infinite contrast when the Hashin and Shtrikman variational principle based scheme is used. 30

(3) When the constitutive matrix is computed by the laminate mixing rule, it needs more iterations to reach convergence. 330

(4) Numerical accuracy is improved if the constitutive matrix of a composite voxel is computed by the laminate mixing rule. (5) Only when the polarized field is sufficiently smooth, high order B-spline improves the numerical accuracy and efficiency.

Acknowledgments 335

This work was supported by the Natural Science Foundation of China (Nos. 11902296, 41731284, 41920104007, 11672360, 41702314 & 51879245), the Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan), (Nos. CUG170682) and the Director Foundation for Key Laboratory of the Ministry of Education (No. CUG2019ZR12).

340

Appendix A. Discretization of the convolution term in Hashin and Shtrikman variational principle Referring to the derivation procedure in Brisard and Dormieux [9] and utilizing the convolution theorem, the Fourier expansion of Equation 16 is as follows ˆ 1 , · · · , ωd ) = sincp+1 πα1 · · · sincp+1 παd φ ˆα ,··· ,α φ(ω 1 d M1 Md

(A.1)

ˆα ,··· ,α is the discrete Fourier transform of φa ,··· ,a , which can be where φ 1 1 d d carried out by fast Fourier transform. Mj is the total voxel number in j direction. The j coordinate at the center of the aj th voxel is xj = angular frequency is ωj =

2π Lj αj .

Lj Mj aj .

The αj th

ix

Due to the periodic property of e , we have

ˆα ,··· ,α ,··· ,α = φ ˆα ,··· ,α +M ,··· ,α φ 1 j 1 j j d d

31

(A.2)

The convolution product Γ0 ∗ φ is expressed as a Fourier series Γ0 ∗ φ(x) =

∞ X

···

α1 =−∞

∞ X

αd =−∞

ˆ 1 , · · · , ωd )ei(ω1 x1 +···+ωd xd ) Γˆ0 (ω1 , · · · , ωd ) : φ(ω (A.3)

Substituting Equation A.3 into the second term of Equation 15, we obtain Z

Ω

δφ : Γ0 ∗ φdV =

∞ X

α1 =−∞

···

Z ∞ X

αd =−∞

δφ :

Ω

(A.4)

ˆ 1 , · · · , ωd )ei(ω1 x1 +···+ωd xd ) dV Γˆ0 (ω1 , · · · , ωd ) : φ(ω Using the Plancherel theorem, we have Z

Ω

δφ : Γ0 ∗ φdV =

L1 · · · Ld (M1 · · · Md )2

∞ X

α1 =−∞

···

∞ X

(A.5)

αd =−∞

ˆ∗ (ω1 , · · · , ωd ) : Γˆ0 (ω1 , · · · , ωd ) : φ(ω ˆ 1 , · · · , ωd ) δφ ˆ∗ is the conjugate of φ. ˆ Substituting Equation A.1 into Equation A.5, where φ we have Z

Ω

δφ : Γ0 ∗ φdV =

L1 · · · Ld (M1 · · · Md )2

∞ X

α1 =−∞

···

∞ X

αd =−∞

πα1 παd ˆ∗ ˆα ,··· ,α sinc2(p+1) · · · sinc2(p+1) δφ : Γˆ0 (ω1 , · · · , ωd ) : φ 1 d M1 Md α1 ,··· ,αd

(A.6)

ˆ 1 , · · · , ωd ), Equation A.6 can be rewritten Accounting for the periodicity of φ(ω as Z

Ω

{

∞ X

m1 =−∞

···

δφ : Γ0 ∗ φdV =

∞ X

md

M M 1 −1 d −1 X X L1 · · · Ld ˆ∗ · · · δφ α1 ,··· ,αd : (M1 · · · Md )2 α =0 α =0 1

d

π(α1 + m1 M1 ) π(αd + md Md ) sinc2(p+1) · · · sinc2(p+1) M1 Md =−∞

2π(α1 + m1 M1 ) 2π(αd + md Md ) ˆα ,··· ,α Γˆ0 [ ,··· , ]} : φ 1 d L1 Ld (A.7)

32

References [1] J.-C. Michel, H. Moulinec, P. Suquet, Effective properties of composite materials with periodic microstructure: a computational approach, Computer 345

methods in applied mechanics and engineering 172 (1-4) (1999) 109–143. [2] V. Vinogradov, G. Milton, An accelerated FFT algorithm for thermoelastic and non-linear composites, International Journal for Numerical Methods in Engineering 76 (11) (2008) 1678–1695. [3] F. Willot, L. Gillibert, D. Jeulin, Microstructure-induced hotspots in the

350

thermal and elastic responses of granular media, International Journal of Solids and Structures 50 (10) (2013) 1699–1709. [4] R. A. Lebensohn, A. K. Kanjarla, P. Eisenlohr, An elasto-viscoplastic formulation based on fast fourier transforms for the prediction of micromechanical fields in polycrystalline materials, International Journal of Plas-

355

ticity 32 (2012) 59–69. [5] F. Bignonnet, L. Dormieux, FFT-based bounds on the permeability of complex microstructures, International Journal for Numerical and Analytical Methods in Geomechanics 38 (16) (2014) 1707–1723. [6] F. Willot, B. Abdallah, Y.-P. Pellegrini, Fourier-based schemes with mod-

360

ified green operator for computing the electrical response of heterogeneous media with accurate local fields, International Journal for Numerical Methods in Engineering 98 (7) (2014) 518–533. [7] V. Monchiet, FFT based iterative schemes for composites conductors with non-overlapping fibers and kapitza interface resistance, International Jour-

365

nal of Solids and Structures 135 (2018) 14–25. [8] K. S. Djaka, S. Berbenni, V. Taupin, R. A. Lebensohn, A FFT-based numerical implementation of mesoscale field dislocation mechanics: Application to two-phase laminates, International Journal of Solids and Structures. 33

[9] S. Brisard, L. Dormieux, FFT-based methods for the mechanics of com370

posites: A general variational framework, Computational Materials Science 49 (3) (2010) 663–671. [10] S. Brisard, L. Dormieux, Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites, Computer Methods

375

in Applied Mechanics and Engineering 217 (2012) 197–212. [11] S. Brisard, Reconstructing displacements from the solution to the periodic Lippmann-Schwinger equation discretized on a uniform grid, International Journal for Numerical Methods in Engineering 109 (4) (2017) 459–486. [12] T. H. Pian, P. Tong, Basis of finite element methods for solid continua,

380

International Journal for Numerical Methods in Engineering 1 (1) (1969) 3–28. [13] Z. Hashin, S. Shtrikman, On some variational principles in anisotropic and nonhomogeneous elasticity, Journal of the Mechanics and Physics of Solids 10 (4) (1962) 335–342.

385

[14] K. S. Eloh, A. Jacques, S. Berbenni, Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains, International Journal of Plasticity 116 (2019) 1–23. [15] A. J. Jerri, The Gibbs phenomenon in Fourier analysis, splines and wavelet approximations, Springer Science & Business Media, 2013.

390

[16] J. Yvonnet, Q.-C. He, C. Toulemonde, Numerical modelling of the effective conductivities of composites with arbitrarily shaped inclusions and highly conducting interface, Composites Science and Technology 68 (13) (2008) 2818–2825. [17] M. Kabel, D. Merkert, M. Schneider, Use of composite voxels in FFT-based

395

homogenization, Computer Methods in Applied Mechanics and Engineering 294 (2015) 168–188. 34

[18] G. W. Milton, The theory of composites, Materials and Technology 117 (1995) 483–493. [19] M. Kabel, A. Fink, M. Schneider, The composite voxel technique for inelas400

tic problems, Computer Methods in Applied Mechanics and Engineering 322 (2017) 396–418. [20] R. Chari`ere, A. Marano, L. G´el´ebart, Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites, International Journal of Solids and Structures 182-183 (2020) 1–14.

405

[21] T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (39-41) (2005) 4135– 4195. [22] M. A. Scott, M. J. Borden, C. V. Verhoosel, T. W. Sederberg, T. J. Hughes,

410

Isogeometric finite element data structures based on B´ezier extraction of tsplines, International Journal for Numerical Methods in Engineering 88 (2) (2011) 126–156. [23] C. De Boor, C. De Boor, E.-U. Math´ematicien, C. De Boor, C. De Boor, A practical guide to splines, springer-verlag New York, 1978.

415

[24] J. M. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and applications, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 289–314. [25] I. Babuˇska, J. M. Melenk, The partition of unity method, International journal for numerical methods in engineering 40 (4) (1997) 727–758.

420

[26] E. Kr¨ oner, Bounds for effective elastic moduli of disordered materials, Journal of the Mechanics and Physics of Solids 25 (2) (1977) 137–155. [27] T. Mura, Micromechanics of defects in solids, Springer Science & Business Media, 2013. 35

[28] J. Willis, Bounds and self-consistent estimates for the overall properties 425

of anisotropic composites, Journal of the Mechanics and Physics of Solids 25 (3) (1977) 185–202. [29] J. H. Lim, D. Sohn, J. H. Lee, S. Im, Variable-node finite elements with smoothed integration techniques and their applications for multiscale mechanics problems, Computers & structures 88 (7-8) (2010) 413–425.

430

[30] D. Ling, L. Bu, F. Tu, Q. Yang, Y. Chen, A finite element method with mesh-separation-based approximation technique and its application in modeling crack propagation with adaptive mesh refinement, International Journal for Numerical Methods in Engineering 99 (7) (2014) 487–521. [31] J.-H. Lv, Y.-Y. Jiao, P. Wriggers, T. Rabczuk, X.-T. Feng, F. Tan, Ef-

435

ficient integration of crack singularities in the extended finite element method: Duffy-distance transformation and conformal preconditioning strategy, Computer Methods in Applied Mechanics and Engineering 340 (2018) 559–576. [32] J.-H. Lv, Y.-Y. Jiao, X.-T. Feng, P. Wriggers, X.-Y. Zhuang, T. Rabczuk,

440

A series of Duffy-distance transformation for integrating 2D and 3D vertex singularities, International Journal for Numerical Methods in Engineering 118 (1) (2019) 38–60. [33] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Computer

445

methods in applied mechanics and engineering 157 (1-2) (1998) 69–94. [34] R. Barrett, M. W. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst, Templates for the solution of linear systems: building blocks for iterative methods, SIAM, 1994.

450

[35] Q.-D. To, M.-T. Nguyen, G. Bonnet, V. Monchiet, V.-T. To, Overall elastic

36

properties of composites from optimal strong contrast expansion, International Journal of Solids and Structures 120 (2017) 245–256. [36] V. Monchiet, G. Bonnet, A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary con455

trast, International Journal for Numerical Methods in Engineering 89 (11) (2012) 1419–1436.

37

Declaration of Interest Statement We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the paper submitted.

Author Contributions The program and manuscript are mainly written by the first author Fubin Tu. The idea comes from the first author Fubin Tu and the second author Yuyong Jiao. Xiaoyong Zhou, Yi Cheng and Fei Tan help to derive the formula, implement the numerical example and draw the figures. All authors have given approval to the final version of the manuscript.