Skip to main content

ORIGINAL RESEARCH article

Front. Mater., 21 April 2023
Sec. Computational Materials Science

A Green’s function-based approach to the concentration tensor fields in arbitrary elastic microstructures

  • 1Institute for Mechanics of Materials and Structures (IMWS), TU Wien (Vienna University of Technology), Vienna, Austria
  • 2Department of Materials Science, Polytechnic University of Madrid (UPM), Madrid, Spain

Computational homogenization based on FEM models is the gold standard when it comes to homogenization over a representative volume element (RVE), of so-called complex material microstructures, i.e., such which cannot be satisfactorily represented by an assemblage of homogeneous subdomains called phases. As a complement to the aforementioned models, which depend on the boundary conditions applied to the representative volume element and which, as a rule, do not give direct access to the macro-micro-relations in terms of concentration tensors, we here introduce a Green’s function-based homogenization method for arbitrary inhomogeneous microstructures: Inspired by the ideas underlying traditional phase-based homogenization schemes, such as the Mori-Tanaka or the self-consistent model, the new method rests on mapping, through the strain average rule, the microscopic strain fields associated with an auxiliary problem to the macroscopic strains subjected to the RVE. Thereby, the auxiliary problem is defined on a homogeneous infinite matrix subjected to homogeneous auxiliary strains and to inhomogeneous (fluctuating) polarization stresses representing the fluctuations of the microstiffness field, i.e., the complex microstructure within the RVE. The corresponding microscopic strains appear as the solution of a Fredholm integral equation, delivering a multilinear operator linking the homogeneous auxiliary strains to the microscopic strains. This operator, together with the aforementioned mapping, eventually allows for completing the model in terms of concentration tensor and homogenized stiffness quantification. This is illustrated by example of a sinusoidally fluctuating microstructure, whereby the corresponding singular convolution integrals are analytically evaluated from the solution of the Poisson’s equation, and this evaluation strategy is then analytically verified through a Cauchy principal value analysis, and numerically validated by a state-of-the-art FFT homogenization procedure. For the given example, the novel analytical method is several thousand times faster than an FTT-based computational homogenization procedure.

Highlights

• Macrostrains are imposed on elastic material volumes with arbitrary inhomogeneous microstructures.

• Corresponding microstrain distributions are determined in a semi-analytical fashion.

• It involves a homogeneous infinite domain with polarization stress distributions.

• This leads to a Fredholm integral equation involving the elastic Green’s function.

• The solution gives access to homogenized stiffness tensors of arbitrary inhomogeneous microstructures.

1 Introduction–Motivation and scope

The main problem in the wide field of micromechanics of materials (Hill, 1963; Zaoui, 2002) is to quantify the effect of mechanical property distribution throughout the microstructures filling a so-called representative volume element (RVE), on the overall mechanical properties of this RVE, i.e., the properties linking the macroscopic strains (being the average over the microscopic strains inside the RVE) to the macroscopic stresses (being the average over the microscopic stresses inside the RVE). Restricting the present contribution to the case of linear elasticity, the problem comprises the following mathematical relations (Zaoui, 2002):

• geometrical boundary conditions prescribed at the boundary of the RVE, SRVE, in the form proposed by Hashin (1983)

u̲x̲=Ex̲,x̲SRVE,(1)

with x̲ denoting the microscopic location vector, with u̲ denoting the microscopic displacement vector, and with E denoting the macroscopic strain tensor, which is independent of the location x̲, see also Table 1. The boundary conditions according to Eq. 1 imply the validity of the strain average rule (Hashin, 1963; Hashin, 1965; Hashin, 1983)

E=1VRVEVRVEεx̲dVx̲=ε,(2)

where ɛ denotes the microscopic linearized strain tensors, defined as the symmetric part of the microscopic gradient of the displacement field u(x̲), i.e.,

εx̲=12gradxu̲x̲+gradxu̲x̲T=gradxSu̲x̲;x̲VRVE,(3)

• the microscopic elastic law being a function of the microstructural position vector x̲

σx̲=cx̲:εx̲,x̲VRVE,(4)

with the microscopic stress tensor σ and the microscopic stiffness tensor c.

• equilibrium conditions

̲xσx̲+f̲x̲=0,x̲VRVE,(5)

where ̲x stands for the nabla operator and f̲ denotes the volume forces;

• the equivalence of macroscopic and microscopic expressions for virtual power densities of internal forces (Jiménez Segura et al., 2022), which, together with Eq. 5, yields the well-known stress average rule as (Hill, 1963; Zaoui, 2002)

Σ=1VRVEVRVEσx̲dVx̲=σ;(6)

• the strain concentration (or downscaling) relation linking, in a multilinear way, the macroscopic to the microscopic strain field (Hill, 1963; Zaoui, 2002)

εx̲=Ax̲:E,x̲VRVE,(7)

where A is the concentration (or downscaling) tensor;

• the macroscopic elastic law, which follows from Eq. 4 and Eqs. 6, 7 as

Σ=Chom:E,(8)

with the homogenized stiffness tensor reading as (Hill, 1963; Zaoui, 2002)

Chom=1VRVEVRVEcx̲:Ax̲dVx̲.(9)

The classical approach for making the problem of Eqs. 19 tractable is to restrict the discussion to Nr homogeneous subdomains or phases within the RVE. Accordingly, the general microstiffness distribution c(x̲) is replaced by a finite number of microstiffness tensors cr,r=1,,Nr, which characterize phases of different shapes, typically represented by means of ellipsoids. The strains in the latter are approximated from the solutions of Eshelby’s matrix-inhomogeneity problem (Eshelby, 1957), and combination of these solutions with the strain average rule specified for a finite number of phases leads to the well-known Mori-Tanaka or self-consistent models (Kröner, 1958; Mori and Tanaka, 1973; Benveniste, 1987; Benveniste et al., 1991), with many applications in a variety of disciplines, including construction and biomedical engineering (Bernard et al., 2003; Hellmich et al., 2004; Hellmich and Mang, 2005; Hofstetter et al., 2005; Fritsch and Hellmich, 2007). In this context, we note that composites with inclusions of different shapes and/or orientations may require additional symmetrization steps guaranteeing the existence of an elastic potential (Sevostianov and Kachanov, 2014; Jiménez Segura et al., 2023).

TABLE 1
www.frontiersin.org

TABLE 1. Mathematical symbols and abbreviations.

Besides this classical approach, it proved useful to introduce, within an RVE, infinitely many (non-spherical) phases, being associated with infinitely many space directions quantified through longitudinal and latitudinal Euler angles, and to associate infinitely many Eshelby problems to each of these directions (Fritsch et al., 2006). After an appropriate discretization of the involved integral expression, the microstrain state within the RVE can be represented sufficiently precisely, so as to allow for upscaling of brittle failure states from the phase-scale, up to the RVE-scale; and this has been shown again for construction and biomedical materials (Fritsch et al., 2009a; Fritsch et al., 2009b; Pichler and Hellmich, 2011; Pichler et al., 2013; Fritsch et al., 2013; Buchner et al., 2022).

Another extension of the classical composite mechanics estimates refers to coated inclusions being embedded in material matrices. A well-known way to approach this problem concerns the extension of Eshelby’s matrix-inclusion problem to a coated inclusion (where the coating may consist of numerous layers, forming an n-layered assemblage (Hervé and Zaoui, 1993; Lipinski et al., 2006)), and then resort to classical combination with the strain average rule, the latter being fed with the inclusion (or core) strains, the layer strains, and the matrix strains. Such models have been very helpful to decipher the mechanical behavior of bone-scaffold composites in tissue engineering (Bertrand and Hellmich, 2009) and of the interfacial transition zone in concrete (Königsberger et al., 2018). Recently, Xu and co-workers proposed a surprisingly simple mathematical alternative to the use of the multiply coated inclusion problem, namely, the repeated use of the Mori-Tanaka estimate, with sequential homogenization of, at a time, one inclusion and one layer playing the role of a matrix; and they successfully applied this strategy to polymer nanocomposites (Xu et al., 2017; Xu et al., 2018), mortar (Xu et al., 2019), and concrete (Xu et al., 2018; Guo et al., 2022).

Still, there is interest in homogenizing over non-uniform stiffness distributions c(x̲); or in other words, over micro-heterogeneous materials with complex microstructures, i.e., such microstructures which cannot be satisfactorily represented by an assemblage of phases as mentioned before. In this context, the most popular approach is based on the Finite Element Method - FEM (Zienkiewicz et al., 2005). It involves discretizing the RVE into very many finite elements, and subjecting it to suitable boundary conditions (Moës et al., 2003; Pahr and Zysset, 2008; Scheiner et al., 2009; Grimal et al., 2011). The latter may be homogeneous, as in Eq. 1, or periodic. This type of analysis, often referred to as computational homogenzation, has been applied to a variety of problems, including woven textures (Moës et al., 2003), bone microstructure (Pahr and Zysset, 2008; Grimal et al., 2011), tissue engineering scaffolds (Scheiner et al., 2009), and fiber-reinforced ultra-high performance concrete (Feng et al., 2022). Corresponding results depend on both the discretization level and the chosen boundary conditions, which requires careful sensitivity analyses to be carried out when aiming at quantitatively reliable results. Also, the computational effort increases with the square of the degrees of freedom, rendering a detailed representation of the microstructure as computationally very expensive. As a remedy to both the discretization and the CPU challenges, FFT-based homogenization schemes based on the Lippmann-Schwinger equation (Lippmann and Schwinger, 1950) have emerged as an interesting alternative to the FEM, in particular so when it comes to image-based computational homogenization (Moulinec and Suquet, 1998; Brisard and Dormieux, 2010; Cai et al., 2019). Such FFT methods are based on a voxel representation of the microstructure, with the elastic properties being constant over one voxel.

Yet, directing our attention back to Eqs. 19, we observe that both FEM and FFT-based homogenization techniques primarily focus on the homogenized stiffness tensor Chom, somewhat neglecting the concentration tensor field A. However, the latter quantity, giving access to microsopic stress and strain fields, is of great interest as well, in particular so as concerns upscaling of elasto-brittle material behavior (Fritsch et al., 2009a; Sanahuja et al., 2010; Fritsch et al., 2013; Königsberger et al., 2018; Wolfram et al., 2022), or of eigenstrains and eigenstresses (Levin, 1967; Rosen and Hashin, 1970; Wang et al., 2018).

This motivates the present paper, presenting a novel way to derive strain concentration tensor fields A(x̲), from a given microstiffness distribution c(x̲) characterizing an arbitrary inhomogeneous microstructure. For this purpose, we adopt a key idea underlying the classical phase-based homogenization approaches such as the Mori-Tanaka or the self-consistent scheme, namely, the introduction of an auxiliary problem defined on an infinite elastic domain, and the suitable combination of such an auxiliary problem with the strain average rule of Eq. 2. In this way, we can resort to fundamental elastic solutions in the form of Green’s functions, while also circumventing the rather awkward dependence of homogenization results on the chosen boundary conditions, as encountered with FEM-based computational homogenization approaches. Accordingly, the paper is organized as follows: Section 2 introduces an auxiliary problem on an infinite elastic domain, and its relation to the strain average rule reflecting the geometrical compatibility throughout the microscopically finite RVE. Section 3 covers a Green’s function-based solution to the auxiliary problem of Section 2. Section 4 presents the first Green’s function-based expression of the strain concentration tensor field in an arbitrary inhomogeneous microstructure. After an illustrative example for a microstructure with sinusoidally fluctuating bulk moduli, given in Section 5, and a Discussion in Section 6, the paper is concluded in Section 7.

2 An auxiliary problem on an infinite domain, and its relation to the RVE

Traditional phase-based micromechanical approaches, such as the Mori-Tanaka or the self-consistent estimates, are built on the solution of Eshelby’s matrix-inhomogeneity problem, where an ellipsoidal inhomogeneity of a specific stiffness is embedded into an infinite matrix of yet another stiffness. The latter matrix is remotely subjected to some auxiliary strains E0. The solution of Eshelby (1957) then relates the auxiliary strains to the homogeneous strains ɛI in the inhomogeneity, according to

εI=A0,I:E0,(10)

with A0,I as the concentration tensor associated with inhomogeneity I embedded in an infinite matrix subjected to E0. A0,I depends on the stiffness contrast between inhomogeneity and matrix, as well as on the shape of the inhomogeneity and its orientation with respect to the material directions of the matrix. In the traditional approach, ɛI is then associated to strains in an ellipsoidal phase inside the RVE, and different inhomogeneities which are all embedded in the same type of matrix are introduced so as to consider different phases within the RVE.

However, as we presently wish to go beyond phase assemblages, we need to extend the auxiliary problem of Eq. 10 beyond homogeneous ellipsoidal domain-related strain ɛI and, instead, introduce a general strain field of the form

εx̲=A0x̲:E0,(11)

with A0 as a concentration tensor field associated with an arbitrary inhomogeneous microstructure represented by a stiffness distribution c(x̲), spreading throughout an infinite domain, see Figure 1C. The stiffness distribution is considered as fluctuation around a homogeneous stiffness C0. The strains ε(x̲) then consist of two portions: (i) auxiliary strains E0 prevailing in the homogeneous infinite domain of stiffness C0, see Figure 1A, and (ii) fluctuations around E0 which arise from the fluctuations in the stiffness field, [c(x̲)C0], see Figure 1B. These strains can be derived from the Green’s functions known for elastic matrices, together with the concept of polarization stresses introduced by Eshelby, as will be detailed in Section 3. We are left with relating the new auxiliary problem of Eq. 11 to the RVE. Therefore, we consider a finite domain of c(x̲) which is statistically representative of the microstructure within the RVE, we identify the volume of this domain with the volume of the RVE, see Figure 1D, and we then apply the strain average rule of Eq. 2, which yields in combination with Eq. 11 that

E=1VRVEVRVEA0x̲dVx̲:E0=M1:E0M=1VRVEVRVEA0x̲dVx̲1,(12)

with M as the RVE-to-auxiliary strain conversion tensor. Multiplication from the left, of Eq. 12 with M, and insertion of the corresponding result into Eq. 11, yields

εx̲=A0x̲:1VRVEVRVEA0x̲dVx̲1:E,(13)

and comparison of Eq. 13 with Eq. 7 yields the strain concentration tensor as

Ax̲=A0x̲:1VRVEVRVEA0x̲dVx̲1=A0x̲:M,(14)

and when considering, in addition, Eq. 9, the homogenized stiffness tensor is eventually retrieved as

Chom=1VRVEVRVEcx̲:A0x̲dVx̲:1VRVEVRVEA0x̲dVx̲1,(15)

see Figures 1E, F. Eq. 15 can be reformulated by setting c(x̲)=C0+c(x̲)C0, which yields

Chom=C0+1VRVEVRVEcx̲C0:A0x̲dVx̲:1VRVEVRVEA0x̲dVx̲1.(16)

We are left with the determination of A0(x̲). Therefore, we will first link this property to the elastic Green’s function (Section 3), and provide a Green’s function-based expression of the strain concentration tensor field in an arbitrary inhomogeneous microstructure (Section 4), before giving an illustrative example (Section 5).

FIGURE 1
www.frontiersin.org

FIGURE 1. Illustration of elements making up the new homogenization scheme: (A) infinite homogeneous elastic matrix of stiffness C0 subjected to background (auxiliary) strains E0; (B) same matrix undergoing equilibrated (polarization) stresses which are equivalent to the effect of fluctuations (c(x̲)C0) in microstiffness, c(x̲), around the homogeneous stiffness C0, (C) sum of load cases (A) and (B); (D) selection of corresponding microstrains and microstiffnesses within a finite domain (RVE) characterizing the microheterogeneous material; (E) strain averaging over RVE, with =1VRVEVRVE(x̲)dV(x̲), and corresponding concentration and homogenization formulae; (F) macroscopic elastic representation of RVE.

3 Green’s function-based solution to the auxiliary problem of the infinite matrix with arbitrary inhomogeneous microelasticity distributions

In order to determine the concentration tensor field A0(x̲) of our new auxiliary problem, we use the method of Green’s functions (Fredholm, 1900; Ting and Lee, 1997). This method is only applicable for homogeneous elastic spaces, which motivates us to adapt a famous idea of Eshelby (1957), which concerns the equivalence of an inhomogeneous elasticity distribution to a homogeneous elastic space subjected to inhomogeneous polarization stresses τ. Mathematically speaking, the constitutive law of Eq. 4 is re-cast into the format (Willis, 1977)

σx̲=C0:εx̲+τx̲.(17)

Equating Eq. 17 with Eq. 4 yields

τx̲=cx̲C0:εx̲,(18)

see Figure 1B. Inserting the displacement-to-strain conversion relation of Eq. 3 into the equivalent constitutive law of Eq. 17, followed by inserting the corresponding result into the equilibrium condition of Eq. 5 yields

̲xC0:gradxSu̲x̲=f̲x̲+̲xτx̲.(19)

The solution of the linear partial differential Eq. 19, u̲(x̲), is the sum of the homogeneous solution u̲h and the particular solution u̲p,

u̲x̲=u̲hx̲+u̲px̲,(20)

whereby

• the homogeneous solution u̲h(x̲) satisfies the homogeneous linear partial differential equation

̲xC0:gradSu̲hx̲=0,(21)

with the inhomogeneous boundary conditions

u̲hx̲=E0X̲x̲,|x̲|;(22)

• and the particular solution u̲p(x̲) satisfies the inhomogeneous linear partial differential equation

̲xC0:gradSu̲px̲=f̲x̲+̲xτx̲,(23)

with homogeneous boundary conditions

u̲px̲=0,|x̲|.(24)

The homogeneous solution reads as

u̲hx̲=E0X̲x̲,x̲R3.(25)

The particular solution can be given in the form

u̲px̲=R3Gx̲y̲f̲y̲+̲yτy̲dVy̲,(26)

where G(x̲y̲) denotes the displacement-related Green’s function tensor, satisfying the differential equation (Kneer, 1965; Ting and Lee, 1997; Tonon et al., 2001; Xie et al., 2016)

̲xC0:gradxSGx̲y̲=1δx̲y̲,(27)

with boundary conditions

Gx̲y̲=0,|x̲|.(28)

In Eq. 27, 1 denotes the second-order identity tensor and δ denotes the Dirac function, with the following properties.

δx̲=0forx̲0,(29)
δx̲=forx̲=0,(30)
R3δx̲dVx̲=1.(31)

The last term in Eq. 26 can be transformed by means of the chain rule, the divergence theorem, and the boundary condition given through Eq. 28, yielding

R3Gx̲y̲̲yτy̲dVy̲=R3gradyGx̲y̲:τy̲dVy̲.(32)

Insertion of Eq. 32 into Eq. 26, adding the respective result to Eq. 25, and using the obtained expression u̲(x̲) in Eq. 3 yield the microscopic strain field throughout the RVE as

εx̲=E0+R3Gx̲y̲f̲y̲dVy̲R3Gx̲y̲:τy̲dVy̲.(33)

In Eq. 33, the following gradients of the Green’s functions were explicitly introduced

Gx̲y̲=gradxSGx̲y̲,(34)
Gx̲y̲=gradxSgradyGx̲y̲.(35)

Considering cases where, in Eq. 33, the effect of polarization stresses clearly outweighs that of volume forces, and inserting the polarization stress expression of Eq. 18 into Eq. 33 yields

εx̲=E0R3Gx̲y̲:cy̲C0:εy̲dVy̲.(36)

Eq. 36 is an implicit integral equation for the microscopic strain field, as the latter appears both as a separate term and part of an integrand in a volume integral over the infinite elastic domain. More precisely, the microstrains are the solution of a Fredholm equation of the second kind (Fredholm, 1900). The solution of Eq. 36 is found by means of an infinitely often repeated substitution process concerning ε(y̲) in the integral of Eq. 36. Accordingly, the term ε(y̲) in the integral of Eq. 36 is numbered in a way reflecting these insertion processes, namely, by ε(y̲(i)), i = 1, 2, , .

As a starting point, Eq. 36 is specified for y̲=y̲(1) in the integral of Eq. 36, yielding

εx̲=E0R3Gx̲y̲1:cy̲1C0:εy̲1dVy̲1.(37)

In order to come up with an expression for ε(y̲(1)) to be inserted into the integral in Eq. 37, we specify Eq. 37 for y̲(1)=y̲(2) and for x̲=y̲(1), yielding

εy̲1=E0R3Gy̲1y̲2:cy̲2C0:εy̲2dVy̲2.(38)

Insertion of Eq. 38 into Eq. 37 yields

εx̲=IR3Gx̲y̲1:cy̲1C0dVy̲1:E0+R3Gx̲y̲1:cy̲1C0:R3Gy̲1y̲2:cy̲2C0:εy̲2dVy̲2dVy̲1,(39)

where I is the symmetric fourth-order identity tensor. We now generalize this idea, in order to come up with the strain field for the (α)-th step of the substitution process, ε(y̲(α)). For this purpose, we specify Eq. 37 for x̲=y̲(α) and y̲(1)=y̲(α+1). Repeating this process over and over again gives access to a relation involving an auxiliary macrostrain concentration tensor field A0(x̲) and a residual term comprising the implicit strain field R[x̲,ε(y̲(N))], reading as

εx̲=A0x̲:E0+Rx̲,εy̲N.(40)

In more detail, the concentration tensor field reads as

A0x̲=I+n=1N1nAn0x̲,(41)

where

A10x̲=R3Gx̲y̲:cy̲C0dVy̲,(42)

and, for n > 1,

An0x̲=R3R3Gx̲y̲1:i=2ncy̲i1C0:Gy̲i1y̲i:cy̲nC0dVy̲ndVy̲1,n2,.(43)

The residual term after N iterations reads as

Rx̲,εy̲N=R3R3Gx̲y̲1:i=2Ncy̲i1C0:Gy̲i1y̲i:cy̲NC0:εy̲N×dVy̲NdVy̲1.(44)

Note that for N, R[x̲,ε(y̲(N))]0, provided the following requirement is met:

R3Gijkℓy̲α1y̲αckmny̲αCkmn0×Gmnpqy̲αy̲α+1dVy̲α<1,y̲α1,y̲α+1VRVE,(45)

for an arbitrary α.

4 Green’s function-based expression of the RVE-related strain concentration tensor field in an arbitrarily inhomogeneous microstructure

After obtaining an analytic expression for the auxiliary strain concentration tensor field A0, see Eq. 41, analytic derivations of the RVE-related quantities are formulated according to Section 2. First, the auxiliary-to-RVE strain concentration follows from Eq. 12. Substitution of x̲ by y̲(0), which is integrated over the volume of the RVE, VRVE, (while every other y̲(i) is integrated over the entire space, R3), allows for the following expression

M=I+n=1N1nVRVEVRVER3R3i=1nGy̲i1y̲i:cy̲iC0×dVy̲ndVy̲01.(46)

Thus, a comprehenive integral format for the strain concentration tensor field is obtained from inserting Eqs. 4143 and Eq. 46 into Eq. 14, yielding

Ax̲=IR3Gx̲y̲:cy̲C0dVy̲+n=2N1nR3R3Gx̲y̲1:i=2ncy̲i1C0:Gy̲i1y̲i:cy̲nC0dVy̲ndVy̲1:I+n=1N1nVRVEVRVER3R3i=1nGy̲i1y̲i:cy̲iC0×dVy̲ndVy̲01.(47)

Lastly, the homogenized stiffness tensor, in terms of Green’s functions, is obtained inserting Eq. 47 into Eq. 15. The resulting expression can be simplified by the substitution of x̲ by y̲(0), yielding

Chom=C0+1VRVEVRVEcy̲0C0dVy̲0+n=1N1nVRVEVRVER3R3i=1ncy̲i1C0:Gy̲i1y̲i:cy̲nC0dVy̲ndVy̲0:I+n=1N1nVRVEVRVER3R3i=1nGy̲i1y̲i:cy̲iC0×dVy̲ndVy̲01.(48)

5 Illustrative example: Microstructure with sinusoidally fluctuating bulk moduli

5.1 Complex microstructure with sinusoidally fluctuating microscopic bulk modulus

In order to illustrate the applicability of the novel integral expressions of Eqs. 4143, we resort to the Green’s function for an infinitely extended isotropic elastic body with bulk modulus k0 and Poisson’s ratio ν0; it reads as (Dvorak, 2012)

Gx̲y̲=1+ν06πk012ν011|x̲y̲|1+ν024πk01ν012ν0gradxgradx|x̲y̲|,(49)

where, with respect to the actual formula (4.5.12) given on page 107 of (Dvorak, 2012), we consider k0=2μ0(1+ν0)3(12ν0), with μ0 being the shear modulus. Moreover, we consider a sinusoidal stiffness distribution across this infinitely extended body, according to

cx̲=C0+3Δksin2πλ1x1sin2πλ2x2sin2πλ3x3Ivol,(50)

where x1, x2, and x3 are the components of location vector x̲ with respect to an orthonormal base frame e̲1, e̲2, e̲3, such that x̲=x1e̲1+x2e̲2+x3e̲3, and Ivol stands for the volumetric part of the symmetric fourth-order unity tensor I. The latter has the components Iijrs = 1/2 (δirδjs + δisδjr), while the former reads as Ivol=1/3(11), with components Iijkℓvol=1/3(δijδk), where 1 is the second-order unity tensor with the Kronecker delta δij as its components. Δk is a parameter which scales the stiffness variance and λi sets the size of one fluctuation in direction i.

5.2 Normal strain-related components of A0(x̲)

The auxiliary strain downscaling tensor A0(x̲) is the result of the sum of an infinite series, see Eq. 41. The first term of this series, A10(x̲), is defined through Eq. 42, so that consideration of Eq. 50 yields the component 1111 of A10 as

A1,11110x̲=h=13=13R3G11hℓhℓx̲y̲3Δksin2πλ1y1×sin2πλ2y2sin2πλ3y3Ih11voldVy̲=Δk=13R3G11x̲y̲sin2πλ1y1×sin2πλ2y2sin2πλ3y3dVy̲.(51)

In order to retrieve the components G1111, G1122, and G1133, we start with the general expression for the components of the Green’s function of Eq. 49, reading as

Gij=1+ν06πk012ν0δij|x̲y̲|1+ν024πk01ν012ν02xixj|x̲y̲|.(52)

The fourth-order Green’s function gradient according to Eq. 35 exhibits the following components

Gijkℓ=1+ν012πk012ν02xjyδik|x̲y̲|+2xiyδjk|x̲y̲|1+ν024πk01ν012ν04xixjxky|x̲y̲|.(53)

In order to evaluate Eq. 53, it is useful to recall the following properties of the spatial derivatives of the norm |x̲y̲|. The first-order derivative of the aforementioned norm reads as

xi|x̲y̲|=xiyi|x̲y̲|,(54)

revealing the interesting property

xi|x̲y̲|=yi|x̲y̲|.(55)

The first derivative of the inverse of the norm |x̲y̲| reads as

xi1|x̲y̲|=xiyi|x̲y̲|3,(56)

revealing the interesting property

xi1|x̲y̲|=yi1|x̲y̲|.(57)

Eq. 55 and Eq. 57 allow for re-writing Eq. 53 as

Gijkℓ=1+ν012πk012ν02xjxδik|x̲y̲|+2xixδjk|x̲y̲|+1+ν024πk01ν012ν04xixjxkx|x̲y̲|,(58)

so that the components occurring in Eq. 51 read as

G1111=1+ν06πk012ν02x121|x̲y̲|+1+ν024πk01ν012ν04x14|x̲y̲|,(59)
G1122=1+ν024πk01ν012ν04x12x22|x̲y̲|,(60)
G1133=1+ν024πk01ν012ν04x12x32|x̲y̲|.(61)

The sum of Eqs. 5961, to be calculated in Eq. 51, can be further simplified through an identity which follows from deriving Eq. 54 with respect to the components of the location vector,

2xi2|x̲y̲|=kixkyk2|x̲y̲|3.(62)

Namely, Eq. 62 entails the following identity

2x12|x̲y̲|+2x22|x̲y̲|+2x32|x̲y̲|=21|x̲y̲|,(63)

twofold derivation of which yields an identity comprising the derivatives of the norm |x̲y̲| occurring in Eqs. 5961, reading as

4x14|x̲y̲|+4x12x22|x̲y̲|+4x12x32|x̲y̲|=22x121|x̲y̲|.(64)

Insertion of Eqs. 5961 into Eq. 51, while considering Eq. 64, yields

A1,11110x̲=Δk1+ν012πk01ν0R32x121|x̲y̲|sin2πλ1y1×sin2πλ2y2sin2πλ3y3dVy̲.(65)

A change of variable to z̲=y̲x̲ results in

A1,11110x̲=Δk1+ν012πk01ν0R32z121|z̲|×sin2πλ1x1+z1sin2πλ2x2+z2×sin2πλ3x3+z3dVz̲.(66)

Transforming Eq. 66 by means of

sina+b=sinacosb+cosasinb,(67)

and considering the second derivative of the inverse of the norm of z̲ as

2z121|z̲|=2z12z22z32|z̲|5,(68)

as well as that the even and odd functions appearing as factors in the integral expression of Eq. 66 imply vanishing integrals

aa2z121|z̲|sin2πλizidzi=0,aR,i1,2,3,(69)

we arrive at

A1,11110x̲=Δk1+ν012πk01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3×R32z121|z̲|cos2πλ1z1×cos2πλ2z2cos2πλ3z3dVz̲.(70)

Noting that z11|z̲z̲|=z11|z̲z̲| and 2z121|z̲z̲|=2z121|z̲z̲|, the integral in Eq. 70 can be expressed by means of an auxiliary function in z̲, according to

R32z121|z̲|cos2πλ1z1cos2πλ2z2cos2πλ3z3dVz̲=2ϕz12z̲=0,(71)

where the auxiliary function ϕ(z̲) stands for

ϕz̲=R31|z̲z̲|cos2πλ1z1cos2πλ2z2cos2πλ3z3dVz̲.(72)

Eq. 72 exhibits two remarkable properties. Firstly, its format

ϕz̲=R31|z̲z̲|fz̲dVz̲(73)

is the solution of the Poisson’s equation

2ϕz̲=4πfz̲.(74)

Secondly, Eq. 72 is symmetric in the sense of

ϕz1,z2,z3=ϕz1,z3,z2=ϕz2,z3,z1=ϕz2,z1,z3=ϕz3,z1,z2=ϕz3,z2,z1ϕz1=ϕz2=ϕz32ϕz12=2ϕz22=2ϕz32.(75)

Use of the latter relations in Eq. 74, while accounting for the structure of f according to Eq. 72 and Eq. 73, yields

2ϕz̲=32ϕz12=4πcos2πλ1z1cos2πλ2z2cos2πλ3z3,(76)

where we made use of Eq. 75. Solving the equation for 2ϕz12(z̲=0) and inserting the corresponding result into Eq. 70 yield

A1,11110x̲=13Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3.(77)

The second term of the series for the auxiliary concentration tensor for the considered sinusoidal isotropic micro-stiffness distribution follows from insertion of Eq. 50 into Eq. 43, so that its component 1111 is obtained as

A2,11110x̲=h=13=13i=13j=13R3R3G11hℓx̲y̲×3Δksin2πλ1y1sin2πλ2y2sin2πλ3y3Ihijvol×m=13n=13Gijmny̲y̲×3Δksin2πλ1y1sin2πλ2y2×sin2πλ3y3Imn11voldVy̲dVy̲.(78)

Specification of Eq. 78 for Eq. 58, while considering the identity of Eq. 64, yields

A2,11110x̲=Δk1+ν012πk01ν02R32x121|x̲y̲|×sin2πλ1y1sin2πλ2y2sin2πλ3y3×2R31|y̲y̲|sin2πλ1y1×sin2πλ2y2sin2πλ3y3dVy̲dVy̲.(79)

Considering the last integral in Eq. 79 as the Green’s function solving the Poisson’s equation

2ϕy̲=4πsin2πλ1y1sin2πλ2y2sin2πλ3y3(80)

yields

A2,11110x̲=4πΔk1+ν012πk01ν02R32x121|x̲y̲|×sin22πλ1y1sin22πλ2y2sin22πλ3y3dVy̲.(81)

Recalling from Eq. 65 and Eq. 77 that

R32x121|x̲y̲|fy̲dVy̲=4π3fx̲,(82)

for any f(y̲) with the symmetry properties of Eq. 75, the integral in Eq. 81 can be solved, yielding

A2,11110x̲=13Δk1+ν03k01ν02sin22πλ1x1×sin22πλ2x2sin22πλ3x3.(83)

Repeating this derivation for the 1111-component of any other member of the series, i.e., for An,11110 with n > 2, one notices that

An,11110x̲=1n13Δk1+ν03k01ν0sin2πλ1x1×sin2πλ2x2sin2πλ3x3n.(84)

Moreover, because of the symmetry of the considered micro-stiffness distribution, An,11110(x̲)=An,22220(x̲)=An,33330(x̲) are equal and given by Eq. 84. Furthermore, it can be straightforwardly proved that Eq. 84 is the result of any component of the type An,ii0.

The components of the auxiliary strain downscaling tensor A0(x̲) are calculated as an infinite sum, see Eq. 41. Therefore, the explicit expression for the sum of an infinite geometric series, reading as

i=0γn=11γ,for|γ|<1,(85)

is applied to α=Δk(1+ν0)3k0(1ν0)sin2πλ1x1sin2πλ2x2sin2πλ3x3. This yields the normal strain concentration tensor components as

Aiiii0x̲=23+11+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3,(86)
Aiijj0x̲=13+11+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3.(87)

5.3 Shear strain-related components of A0(x̲)

Next, we turn to the shear-related components of the concentration tensor, i.e., to Aijkℓ with ij or k. The concentration tensor is driven by the micro-stiffness fluctuation around the base stiffness C0, see the expressions of Eqs. 4143, so that the chosen sinusoidal micro-stiffness distribution of Eq. 50, where the shear stiffness does not fluctuate around the basic contribution, implies that

An,ijk0x̲=0,nN,x̲VRVE,k.(88)

Combining Eq. 88 with Eq. 41 implies that

Aijij0x̲=1,x̲VRVE,ij,(89)

and

Aijkℓ0x̲=0,x̲VRVE,ij,k,ik,j.(90)

Let us now turn to the remaining non-vanishing shear-related components of the concentration tensor, i.e., to A1,ij0(x̲), with ij. For the microstiffness distribution of Eq. 50, the respective first term in the series of Eq. 41, defined by Eq. 42, reads as

A1,ij0x̲=h=13k=13R3G11hkx̲y̲3Δksin2πλ1y1×sin2πλ2y2sin2πλ3y3IhkvoldVy̲=Δkh=13R3G11hhx̲y̲sin2πλ1y1×sin2πλ2y2sin2πλ3y3dVy̲.(91)

Insertion of Eq. 58 into Eq. 91, while considering the identity resulting from derivation of Eq. 63 with respect to xi and xj, reading as

4xixjx12|x̲y̲|+4xixjx22|x̲y̲|+4xixjx32|x̲y̲|=22xixj1|x̲y̲|,(92)

yields

A1,ij0x̲=Δk1+ν012πk01ν0R32xixj1|x̲y̲|×sin2πλ1y1sin2πλ2y2sin2πλ3y3dVy̲.(93)

Introducing the variable z̲=y̲x̲, using the relation of Eq. 67, and considering

zjzij1|z̲|=3zizj|z̲|5,(94)

and the corresponding consequence of even and odd functions

aazi1|z̲|cos2πλizidzi=0,aR,(95)

Eq. 93 can be transformed to

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixicos2πλjxj×sin2πλkxkR33zizj|z̲|5sin2πλizi×sin2πλjzjcos2πλkzkdVz̲,(96)

with i,j,l = 1,2,3, whereby ij and ki,j. Introducing the microstructure-related variable change z̲=2πλ1z1,2πλ2z2,2πλ3z3 into Eq. 96 yields

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixi×cos2πλjxjsin2πλkxk×R33zizj|z̲|5sinzisinzjcoszkdVz̲,(97)

with i,j,l = 1,2,3, whereby ij and ki,j. Considering

zizi2+zj2+zk25/2sinzidzi=2K1zj2+zk23zj2+zk2,(98)

with K1 being the modified Bessel function of the second kind, Eq. 97 can be re-written as

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixicos2πλjxjsin2πλkxk×2zjzj2+zk2×K1zj2+zk2sinzjcoszkdzjdzk,(99)

with i,j,l = 1,2,3, whereby ij and ki,j. Applying a transformation towards polar coordinates, zj=ρcosθ and zk=ρsinθ, yields

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixicos2πλjxjsin2πλkxk×02π02ρcosθK1ρsinρcosθcosρsinθdρdθ,(100)

with i,j,l = 1,2,3, whereby ij and ki,j. Thanks to the trigonometric identity sin(a)cos(b)=12[sin(a+b)+sin(ab)], Eq. 100 can be reformulated as

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixicos2πλjxjsin2πλkxk×02πcosθ0ρK1ρsincosθ+sinθρ+sincosθsinθρdρdθ,(101)

with i,j,l = 1,2,3, whereby ij and ki,j. Thus, integrating over ρ yields

A1,ij0x̲=Δk1+ν012πk01ν0cos2πλixicos2πλjxjsin2πλkxk×02πcosθcosθ+sinθ2+2cosθsinθ+cosθsinθ22cosθsinθ+sinh1cosθ+sinθ2+2cosθsinθ3/2+sinh1cosθsinθ22cosθsinθ3/2dθ,(102)

with i,j,l = 1,2,3, whereby ij and ki,j. The integral remaining in Eq. 102 amounts to 4π3. Thus, the first element of the series is

A1,ij0x̲=4π3Δk1+ν012πk01ν0cos2πλixi×cos2πλjxjsin2πλkxk,(103)

with i,j,l = 1,2,3, whereby ij and ki,j.

The following element of the series A2,ijll0(x̲) reads, after the introduction of the corresponding Green’s functions from Eq. 58, as

A2,ij0x̲=Δk1+ν012πk01ν02R32xixj1|x̲y̲|×sin2πλ1y1sin2πλ2y2sin2πλ3y3×2R31|y̲y̲|sin2πλ1y1×sin2πλ2y2sin2πλ3y3dVy̲dVy̲.(104)

After application of Poisson’s Eq. 74 and Eq. 104 can be expressed as follows, when using sin2(a)=12[1cos(2a)]

A2,ij0x̲=Δk1+ν012πk01ν0218R32xixj1|x̲y̲|×1cos22πλ1y11cos22πλ2y2×1cos22πλ3y3dVy̲.(105)

Proceeding with an analogous process to the one carried out to obtain A1,ijll0(x̲), i.e., variable change to z̲=y̲x̲, partial derivation with respect to zi and zj, change of the variable z̲=2πλ1z1,2πλ2z2,2πλ3z3, integration over zi, change towards polar coordinate system (zj=ρcosθ and zk=ρsinθ), integration over r and, lastly, integration over θ yield

A2,ij0x̲=Δk1+ν012πk01ν02sin22πλixisin22πλjxj×π4π6cos22πλkxk,(106)

with i,j,l = 1,2,3, whereby ij and ki,j. Similarly, the third element was computed as

A3,ij0x̲=Δk1+ν012πk01ν039π16cos2πλixicos2πλjxjsin2πλkxk27π176cos32πλixicos2πλjxjsin2πλkxk27π176cos2πλixicos32πλjxjsin2πλkxk9π176cos2πλixicos2πλjxjsin32πλkxk+27π304cos32πλixicos32πλjxjsin2πλkxk+9π304cos32πλixicos2πλjxjsin32πλkxk+9π304cos2πλixicos32πλjxjsin32πλkxkπ48cos32πλixicos32πλjxjsin32πλkxk,(107)

with i,j,l = 1,2,3, whereby ij and ki,j.

For the following terms An,ijll0(x̲), corresponding substitution of Green’s function tensor components of Eq. 58 and reiterated application of Poisson’s Eq. 74 yield

An,ij0x̲=Δk1+ν012πk01ν0nR32xixj1|x̲y̲|×sinn2πλ1y1sinn2πλ2y2sinn2πλ3y3dVy̲.(108)

These terms can be computed in the same manner as the previous ones, converging rapidly due to the factor Δk(1+ν0)12πk0(1ν0)n.

5.4 Tensorial link between auxiliary and real macrostrains: Access to strain concentration tensor field

For the present case, the RVE is regarded as any assembly of a finite number of fluctuations which are periodically repeated, i.e.,

x̲VRVE=nλ3λ1λ2λ3nλλ1/2x1nλλ1/2,nλλ2/2x2nλλ2/2,nλλ3/2x3nλλ3/2,(109)

with nλ as the number of stiffness waves along edge directions of a box-shaped RVE with a sinusoidal microstructure. The concentration tensor associated with this microstructure, A(x̲), is related to the auxiliary concentration tensor calculated in the previous sections, A0(x̲), according to Eq. 14. Thus, this section is devoted to the derivation of the tensorial link M, see Eq. 12. For the sake of simplicity, the inverse of M, M1, will be calculated first, in order to obtain the RVE-to-auxiliary strain conversion tensor as

M=M11.(110)

Like in the previous sections, the components of tensor M1 will be obtained individually. The first components studied are M1iiii. Insertion of Eq. 86 into the inverse of Eq. 12 yields

(M1)iiii=231VRVEVRVEdVx̲+1VRVEVRVE11+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3dVx̲.(111)

Clearly, the term in square brackets is equal to 1, while we must focus on the other integral expression

I=1VRVEnλλ1/2nλλ1/2nλλ2/2nλλ2/2nλλ3/2nλλ3/211+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3×dx1dx2dx3=1λ1λ2λ3λ1/2λ1/2λ2/2λ2/2λ3/2λ3/211+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3×dx1dx2dx3.(112)

Proceeding with a change of variable x̲=2πλ1x1,2πλ2x2,2πλ3x3 yields

I=12π3ππππππ11+Δk1+ν03k01ν0sinx1sinx2sinx3×dx1dx2dx3.(113)

Solving the integral for x1 yields

I=12π3ππππ2tan1ψx2,x3+tanπ/21ψ2x2,x31ψ2x2,x3dx2dx3ππππ2tan1ψx2,x3+tanπ/21ψ2x2,x31ψ2x2,x3dx2dx3=12π2ππππ11ψ2x2,x3dx2dx3,(114)

where ψ(x2,x3)=Δk(1+ν0)3k0(1ν0)sin(x2)sin(x3). Integrating I with respect to x2 yields

I=22π2ππKχ2sin2x3+11χ2sin2x3Kχ2sin2x31χ2sin2x3dx3,(115)

where K(a) is the complete elliptic integral of the first kind with parameter a, and

χ=Δk1+ν03k01ν0(116)

is characteristic for each microstructure. The value of I, see Eq. 115, has been obtained numerically by means of different integration methods, including the trapezoidal or Simpson’s rule (Whittaker and Robinson, 1967; Horwitz, 2001), for several values of χ. The value of I is equal to 1 for χ = 0, and increases non-linearly with increasing χ, see Figure 2. Thus, from Eq. 111,

FIGURE 2
www.frontiersin.org

FIGURE 2. Numerical value of I, see Eq. 115, obtained for several values of χ=Δk(1+ν0)3k0(1ν0).

M1iiii=23+Iχ.(117)

The next components to be considered are Miijj1, with ij. They read, from Eq. 87 and Eq. 115, as

M1iijj=13+Iχ,ij.(118)

The shear components read as

M1ijij=1,ij,(119)

and

M1ijkℓ=0,ij,ik.(120)

Thus, the auxiliary-to-RVE tensor reads as

M=13I1/3+2I1/3I1/3I0001/3I1/3+2I1/3I0001/3I1/3I1/3+2I0000003I0000003I0000003I.(121)

Lastly, the real strain concentration tensor field is computed according to Eq. 14 as

Ax̲=AiiiiAiijjAiijj000AiijjAiiiiAiijj000AiijjAiijjAiiii000000100000010000001,(122)

where

Aiiiix̲=23+13Iχ11+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3,(123)
Aiijjx̲=13+13Iχ11+Δk1+ν03k01ν0sin2πλ1x1sin2πλ2x2sin2πλ3x3,(124)

see Eq. 115 and Figure 2 for I(χ).

5.5 Homogenized stiffness of sinusoidally fluctuating microstructure

From Eq. 16, the difference between the homogenized stiffness and the background stiffness C0 reads as

ΔChom=ChomC0=1VRVEVRVEcx̲C0:Ax̲dVx̲.(125)

One more time, the components of this tensor will be obtained individually. The only non-vanishing components ΔChom,ijkℓ are those with i = j and k = , due to Eq. 122 and

cijkℓx̲Cijkℓ0=Δksin2πλ1x1sin2πλ2x2sin2πλ3x3δijδk.(126)

Therefore, the components ΔChom,iikk read as

ΔChom,iikk=1VRVEVRVEΔksin2πλ1x1sin2πλ2x2×sin2πλ3x3p=13Appkkx̲dVx̲.(127)

Thus, inserting Eq. 123 and Eq. 124 into Eq. 127 and applying x̲=2πλx̲ yields

ΔChom,iikk=ΔkI2π3VRVEsinx1sinx2sinx31+Δk1+ν03k01ν0sinx1sinx2sinx3dVx̲.(128)

Integrating Eq. 128 with respect to x1 from −π to π yields

ΔChom,iik=3k01ν0I1+ν02π2ππππdx2dx3ππππ11ψ2x2,x3dx2dx3=3k01ν01+ν0Iχ1Iχ,(129)

whereby we have made use of Eq. 114. Clearly, for the stiffness field of Eq. 50, the remaining components vanish, considering Eq. 122 and Eq. 126.

6 Discussion

It is worthwhile to discuss key characteristics of the convolution integral-type mathematical relations for the concentration tensor fields introduced in the present paper. We start by mentioning that Eq. 48 is, to our best knowledge, the first-ever explicit integral formulation for the concentration tensor fields arising from a general microstiffness field. Actually, in the pertinent literature, see, e.g., the review of Zaoui (2002), the mathematical existence of such a field is mentioned, while any corresponding explicit expression is missing.

Our new expression, Eq. 48, provides a common basis for the treatment for virtually any type of microstructure, be they general harmonic microstructures, with microelasticity distributions which can be represented by means of Fourier series, as described in Section 6.5, or even distributions with discontinuities, allowing for the consideration of classical composite material morphologies, arising from the assembly of a finite number of microstructural domains with uniform stiffnesses, normally referred to as “material phases”. Accordingly, the novel method is also apt for large volume fractions of one of the aforementioned domains, i.e., to the so-called “large concentration composition”.

As it is well known that the Green kernel occurring in the aforementioned integral expressions is singular at the point x̲=y̲, the corresponding convolutions need to be carried out with care, and it is therefore interesting to compare the solution strategy based on the Poisson’s equation, as applied throughout Section 5, with the more traditional way of evaluating such integrals, namely, by introducing an infinitesimally small sphere around the singularity, and by transforming the volume integral within that sphere to a surface integral across the sphere’s surface. This will be covered in the first subsection of the present Discussion section.

Alternatively, one may wish a numerical confirmation of our new analytical approach to strain concentration tensor fields, and we provide such a confirmation in terms of FFT-based computational homogenization, in the second subsection of the present Discussion section.

It is also instructive to compare our approach to earlier suggestions for the use of a Fredholm integral equation similar to Eq. 36, often referred to as the Lippmann-Schwinger equation; in particular so concerning the domain over which the convolution integral is evaluated, the type of polarization field considered, and the relation of the Fredholm integral equation to the macroscopic strain associated with the RVE. This is the topic of the third subsection of the present Discussion section.

Finally, we discuss the range of validity of the Fredholm integral Eq. 36, and the practical evaluation of the concentration tensor expressions in the case of microstructures which are more general than that with the sinusoidally fluctuating microscopic bulk moduli covered in Section 5. Corresponding deliberations conclude the present Discussion section.

6.1 Singular convolution integrals–Analytical validation by means of Cauchy principal value

All integral expressions defining concentration tensor fields, such as Eq. 14, Eqs. 4143, Eq. 51, Eq. 65, and Eq. 70, exhibit singularities at x̲=y̲, i.e., at x̲y̲=z̲=0. In Section 5, we circumvented a direct treatment of this singularity, when evaluating Eq. 70 from the solution of the Poisson’s equation, in combination with an auxiliary function in z̲. In order to check the relevance of this strategy, we here evaluate the integral in Eq. 70 by an alternative approach, sometimes referred to as the Cauchy principal value analysis. Therefore, the integral in Eq. 70 is split into two portions associated with two integration domains: The first one is a sphere around the singular point (with a variable radius ϵ, eventually tending towards zero), and the second one is the remaining (unbounded) three-dimensional space.

Denoting the small spherical domain as Vϵ, the integral in Eq. 70 can be recast as, see, e.g., Buryachenko (2007), p. 54,

Vϵ2z121|z̲|fz̲dVz̲=Vϵ2z121|z̲|fz̲f0dVz̲+Vϵ2z121|z̲|dVz̲f0,(130)

whereby

fz̲=cos2πλz1cos2πλz2cos2πλz3,(131)

is fully in line with the developments of Eq. 72 and Eq. 73. As stated before, we are interested in the limit case of ϵ → 0 where

limϵ0Vϵ2z121|z̲|fz̲f0dVz̲=0,(132)

so that

limϵ0Vϵ2z121|z̲|fz̲dVz̲=limϵ0Vϵ2z121|z̲|dVz̲f0=limϵ0Sϵz11|z̲|n1dSz̲f0,(133)

whereby we made use of the divergence theorem, with n̲ standing for the outward normal onto the spherical surface. Notably, the surface integral in Eq. 133 does not exhibit any singularity any more, since the radius ϵ, however small it may become, never actually reaches zero, so that the integrand in the last integral of Eq. 133 always stays finite. Let us evaluate the latter in more detail: Realizing that

ϵ=|x̲y̲|=|z̲|=|z̲|=z12+z22+z321/2,(134)

the integrand in the surface integral of Eq. 133 can be transformed to

z1z12+z22+z321/2=z1z12+z22+z323/2=z1ϵ3.(135)

Then, the surface integral in Eq. 133 is preferably evaluated in spherical coordinates, where.

z1=ϵcosϕsinθ,(136)
n1=cosϕsinθ,(137)
dS=ϵ2sinθdϕdθ,(138)

so that use of Eqs. 135138 in the surface integral of Eq. 133 yields an expression which becomes independent of ϵ, and hence, of the limiting process. In mathematical detail, we have

limϵ0Sϵz11|z̲|n1dSz̲f0=θ=0πsin3θdθϕ=02πcos2ϕdϕf0=4π3f0,(139)

and hence, since f (0) = 1,

Vϵ2z121|z̲|fz̲dVz̲=4π3(140)

The result of Eq. 140 is fully equivalent with evaluating Eq. 76 for z̲=0 and inserting the corresponding result into Eq. 71. This proves the solution strategy for singular integrals, as given through Eqs. 7177. Accordingly, the small spherical integration domain yields the solution of the entire volume integral (spanning also the entire three-dimensional space outside the small spherical domain); hence, the integral of Eq. 65, when evaluated over the three-dimensional space expect for the small sphere enclosing the singularity at the origin z̲=0, vanishes. This last statement can also be found in the book of Buryachenko (2007), namely, as the last equation of (3.29) in the aforementioned reference.

The situation is totally different when it comes to the shear-related concentration tensor components according to Eq. 93. There, the Cauchy principal value needs to be multiplied by sin (0) = 0 so that the integration over the small sphere delivers zero, and it is the domain outside the small sphere, which solely contributes to the integral in Eq. 93. A procedure for solving this regular integral was presented, see Eqs. 94103.

6.2 Strain concentration tensor fields and homogenized stiffness: Numerical confirmation by means of FFT homogenization

In order to gain further confidence into our novel method, we evaluate the analytically defined strain concentration tensor fields of Eq. 123 and Eq. 124 for a particular numerical choice of sinusoidal microelasticity field, see Table 2 for this choice, and then compare this evaluation to suitably chosen microstrain fields computed by means of the classical FFT homogenization methods proposed by Moulinec and Suquet (1994), Moulinec and Suquet (1998), and widely expanded and used thereafter (Lucarini et al., 2021).

TABLE 2
www.frontiersin.org

TABLE 2. Mechanical properties and quantities associated to the particular microelastic material, see Eq. 50.

The key idea of FFT homogenization is to start with an estimate for the microstrain fields, to compute a corresponding estimate for the polarization stress field, then test the latter estimate on the Fourier transform of Eq. 33 with f̲0, reading as

ε̂k̲=δk̲E0Ĝk̲:τ̂k̲,(141)

with the wave vector k̲, and the hat-symbol indicating the Fourier transform, according to

ε̂k̲=12π3R3εx̲expik̲x̲dVx̲.(142)

Namely, prescribing the polarization stress on the right-hand side of Eq. 141 yields a new estimate of microstrain field in the wave domain, which needs to be back-transformed into the traditional location space domain. This new estimate allows for computation of an improved estimate of the polarization stress fields, as described before. The sequential estimation of the polarization stress field (in the location space) and the microstrain field (in the wave space) is repeated until two successive strain estimates differ only very slightly from each other. This algorithm becomes especially appealing if both the Fourier transform and the inverse Fourier transform are realized discretely, on a finite number of locations given by the voxels making up an image. As such a discrete transform can be set up to deliver results in a particularly fast fashion, it is called Fast Fourier transformation, abbreviated as FFT (Cochran et al., 1967), and accordingly, the aforementioned algorithm is referred to as FFT homogenization.

In order to get access to the strain concentration tensor field (which is not the standard target of an FFT homogenization study), we identify the component Aiiii, with i = 1, 2, 3, as the microstrain field component ɛii, with i = 1, 2, 3, which arises from a macroscopic strain tensor of the format E=e̲ie̲i, with e̲i denoting the base vector pointing into the ith direction of an orthonormal base frame. In other words, the only non-vanishing component of the aforementioned macroscopic strain tensor is Eii = 1. Analogously, component Aiijj, with i, j = 1, 2, 3 and ij, is the microstrain field component ɛii, with i = 1, 2, 3, arising from a macroscopic strain tensor of the format E=e̲je̲j.

Different, increasingly fine, FFT discretizations of the investigated sinusoidal microelasticity distribution deliver strain tensor concentration fields which very satisfactorily converge to the analytical solution of Eq. 123 and Eq. 124, see Figure 3 and Figure 4. Accordingly, the FFT-determined homogenized stiffness properties agree virtually perfectly with the analytical results according to Eq. 129, see Table 3. At the same time, our new analytical solution strategy is remarkably efficient from the viewpoint of CPU time: on a Core i5-1035G1 processor in a Lenovo Ideapad S340-15IIL computer, the numerical evaluation of the integral I(χ) according to Eq. 115, the only operation needing non-negligible computer time, lasts for only 0.004 s, while FFT solution processes for 1003 voxels last by a factor of over 6,500 longer, namely, for 26.7 s.

FIGURE 3
www.frontiersin.org

FIGURE 3. Strain concentration field Aiiii(x̲), i = 1, 2, 3, arising from the sinusoidal microelasticity distribution of Eq. 50, at the plane xi = λi/4, i = 1, 2, 3, obtained (A) analytically from Eq. 123, and numerically by means of FFT (Lucarini et al., 2021) with (B) 323 voxels, (C) 523 voxels, and (D) 1003 voxels.

FIGURE 4
www.frontiersin.org

FIGURE 4. Strain concentration field Aiijj(x̲), i, j, = 1, 2, 3, at the plane xi = λi/4, i, j, = 1, 2, 3, obtained (A) analytically from Eq. 124, and numerically by means of FFT (Lucarini et al., 2021) with (B) 323 voxels, (C) 523 voxels, and (D) 1003 voxels.

TABLE 3
www.frontiersin.org

TABLE 3. Numerical value of stiffness tensor components of the benchmark example calculated by means of the proposed analytical model, see Eq. 129, and by FFT with different discretizations.

6.3 Use of the Lippmann-Schwinger equation: Auxiliary problems, integration domains, and macroscopic strains associated with the RVE

From a terminological viewpoint, we note that equations of the format of Eq. 36, irrespective of the chosen integration domains or the format of the polarization stresses, are often referred to as the Lippmann-Schwinger equation, as a similar equation has been proposed by Lippmann and Schwinger (1950) in the field of quantum mechanics. In this context, it is interesting to compare our present contribution to earlier micromechanical applications of the Lippmann-Schwinger equation. A form which is virtually identical to Eq. 36 appears as Eq. 9 in (Molinari and El Mouden, 1996); except for a sign change stemming from the definition of the fourth-order Green operator as the twofold gradient with respect to x̲ – this differs from our definition (35). Mathematically speaking, the aforementioned sign change is due to

gradxgradxGx̲y̲=gradxSgradyGx̲y̲.(143)

However, the actual use of Eq. 9 in (Molinari and El Mouden, 1996) is quite different from our present use of Eq. 36, as Molinari and El Mouden (1996) introduce an infinite number of uniform subfields of microscopic stiffnesses, representing strongly interacting “elastic inclusions” within an RVE of a composite material. At the same time, as a certain commonality of the approach of Molinari and El Mouden (1996) and our present contribution, we note that the latter authors’ strain ɛ0 plays exactly the role of our auxiliary strain E0: it is the strain applied to the auxiliary homogeneous, infinite matrix which undergoes polarization stresses. However, different from our approach to solve this auxiliary problem so as to provide the microscopic strains as a function of the auxiliary strains, Molinari and El Mouden (1996) apply the strain average rule directly to the Fredholm integral equation, i.e., to their Eq. 9, and this leads to their Eq. 16 linking auxiliary and RVE-related strains. On this basis, they then discuss explicit solutions for finite numbers of inclusions within a periodically repeating cubic cell. In this context, Molinari and El Mouden (1996) apply the strain average rule to an infinite domain, as can be seen from their Eq. 15, while our Eq. 12 is clearly related to the (finite) RVE, and hence to the average rule of Eq. 2, which arises from the Hashin displacement boundary conditions imposed onto the RVE according to Eq. 1. We also note that neither Eq. 9 nor Eq. 16 in (Molinari and El Mouden, 1996) give access to the concentration tensor fields - so that our expression according to Eq. 14, together with Eqs. 4143, turn out as an interesting original aspect of the present paper.

Eq. 36 of the present paper is also reminiscent of Eq. 2.28 in (Torquato, 1997). However, different from our approach, Torquato (1997) restricts a non-vanishing polarization field to a finite domain within his infinitely large auxiliary problem subjected to some auxiliary strain ɛ0, the role of which is comparable to our auxiliary strain E0. In this context, he notes that the result of the corresponding convolution integral depends on the shape of the aforementioned finite domain, a situation which does not occur in our anaylsis in which the convolution integrals are evaluated throughout the unbounded auxiliary matrix. Eventually, Torquato (1997) lets his finite polarization domain coincide with the RVE of an anisotropic two-phase composite, for which he identifies stiffness series expansions in powers of “elastic polarizabilities”.

A further difference appears between Eq. 36 of the present contribution and the formally similar Eq. 4 of the famous paper of Moulinec and Suquet (1994), see also (Moulinec and Suquet, 1998). The latter authors introduce the convolution integral directly on the (finite) RVE, noting that a corresponding explicit Green’s function can only be given in the case of periodic displacement boundary conditions imposed onto the RVE, and of corresponding microscopic strains which fluctuate periodically around their average (i.e., around the macroscopic strain). Namely, it is under this periodicity condition, that an explicit solution for the convolution problem exists in the Fourier space, which, in turn, allows for the development of a very efficient algorithm for the mechanical treatment of images made up of pixels or voxels, with the polarization stress being constant throughout one pixel or voxel, respectively.

Green’s operators in convolution integrals over a finite volume (i.e., differing from our present integration over an infinite auxiliary domain) have been already introduced in the 1970s: In this context, Zeller and Dederichs (1973) noted that the corresponding Green’s functions read as G(x̲,y̲)=G(y̲,x̲), rather than G(x̲y̲)=G(y̲x̲), an aspect which was somewhat overlooked by Korringa (1973). However, explicit expressions for the aforementioned Green’s functions are not available, so that Zeller and Dederichs (1973) restricted their analysis to series expansions for small stiffness fluctations, while Kröner (1977) uses convolution integrals over finite volumes for the derivation of bounds for the effective elastic moduli of disordered materials.

Our iterative scheme for solving the Fredholm integral Eq. 36 also bears some similarities with earlier contributions in the field: Kröner (1977) presents an iterative solution for the Lippmann-Schwinger equation formulated directly on the RVE, and Torquato (1997) proposes an iterative scheme which finally delivers the polarization stress as a function of the homogeneous auxiliary strain ɛ0.

6.4 Range of validity of Lippmann-Schwinger equation

The practical relevance of the case where the polarization stresses in Eq. 33 outweigh the effect of the microscopic volume forces, which is the prerequisite for the Lippmann-Schwinger Eq. 36 to hold, deserves further discussion: Within the RVE, the microscopic stresses σ fluctuate around their spatial average, which is the macroscopic stress Σ, and the characteristic length scale d of this fluctuation is scale-separated from the length of the RVE, RVE, which reads mathematically as

σΣx̲=σx̲,withσσx̲=d.(144)

Due to the mathematical structure of the microscopic equilibrium conditions, see Eq. 5, any microscopic volume forces leading to microscopic stress fluctuations are required to change their sign, i.e., their direction, over distances as small as d. Practically speaking, this is an exceptional case: Even in composites with high contrast in mass density, the corresponding gravitational forces of varying magnitude would always share the same direction; or in other words, practically relevant force fields are often parallel within the RVE. We note in passing, that such micro-parallel force field, directly implying the validity of Eq. 36, even fulfill a force field average rule (Jiménez Segura et al., 2022).

6.5 Practical note concerning generally harmonic microstiffness fluctuations

Finally, we discuss which aspects of Section 5 hold beyond the restriction to sinusoidally fluctuating microstiffnesses, and how the semi-analytical solutions presented in this section may be generalized to generally harmonic microstiffness fluctuations. In this context, the key generalization step would be the representation of any, arbitrarily general continuous microstiffness distribution across a finite RVE by a three-dimensional Fourier series. Generalizing, in this way, the example distribution of Eq. 50 to an arbitrarily inhomogenous bulk modulus distribution Δk (x1, x2, x3) yields

cx̲=C0+3klmcklmexp2πikx1λ1+lx2λ2+mx3λ3,(145)

with i now standing for the imaginary unit, and with the Fourier coefficients cklm being obtained from

cklm=aaaaaaexp2πikx1λ1+lx2λ2+mx3λ3×Δkx1,x2,x3dx1dx2dx3.(146)

In other words, sums of products of any three trigonometric functions, be they sine or cosine,—rather than the product of three sine or three cosine terms—would occur throughout the convolution integrals. However, the effect on the corresponding modifications of Eq. 66 and Eq. 93 are only minor: Thanks to Eq. 67 and

cosa+b=cosacosbsinasinb,(147)

the structure of the integrals in Eq. 70 and Eq. 96 stays unaffected. This shows the considerable potential of our method for material investigation based on Fourier-representation of images (Chung et al., 2007)—as an interesting complement to the popular voxel-based FFT schemes.

7 Conclusion

In the present paper, it is for the first time that explicit integral expressions for the concentration tensor fields arising from generally non-uniform microelasticity distributions have been given. More precisely, the effects of elastic behavior at the microscale are represented by means of the Green’s function formalism, leading to Fredholm integral equations which provide novel, series-type integral expressions for the concentration tensor field. The latter may be analytically solved in cases where the involved integral expressions can be formulated as derivatives of the solutions of the Poisson equation, then providing an unprecedented direct access to macro-to-micro scale transition relations, as expressed by the concentration tensor expressions of Eq. 47, and Eqs. 122124. This opens new avenues for exploring the mechanical effect of eigenstrains in hierarchical material systems with complex morphologies, as an interesting alternative to classical computational homogenization. The new approach also provides semi-analytical access to the homogenized stiffness, such as that calculated for a microstructure with sinusoidally fluctuating bulk moduli, see Eq. 16 and Eq. 129. Since I(χ) ≥ 1, see Figure 2, the resulting homogenized stiffness, Chom, is smaller than or equal to the average stiffness, C0. This is fully consistent with the famous result of Voigt (1889) that the average over the microstiffness is larger than the homogenized stiffness, in the sense that (Zaoui, 2002)

E:AT:c:AChom:E0,E.(148)

Besides its obvious fundamental knowledge-related and conceptual merits, the new methods appears as extremely advantegeous from a computational viewpoint, in particular if the microstructure can be represented by just a few elements of the series given in Section 6.5, because relevant computational power is needed only for the one-time determination of the Fourier coefficients and the value of the integral occurring in the function graphed in Figure 2; whereby the latter my be even represented by some suitably chosen fitting function.

Nevertheless, as regards classical composite microstructures, well-known homogenization methods based on the Eshelby-Laws matrix-inhomogeneity problem, such as suitable generalizations of the Mori-Tanaka method (Benveniste, 1987), accounting for symmetrization strategies if required (Sevostianov and Kachanov, 2014; Jiménez Segura et al., 2023), may turn out as more efficient than an approach starting from the general expression given in Eq. 48 and Eq. 49. This underlines the sustained success of advanced composite mechanics in the classical sense.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

NJ: investigation, methodology, formal analysis, computation, visualization, writing—original draft, review and editing. BP: supervision, conceptualization, funding acquisition, methodology, formal analysis, writing—original draft, writing—review and editing. CH: supervision, conceptualization, funding acquisition, methodology, formal analysis, writing—review and editing.

Funding

The authors gratefully acknowledge financial support in the framework of the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764691. The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Benveniste, Y. (1987). A new approach to the application of Mori-Tanaka’s theory in composite materials. Mech. Mater. 6, 147–157. doi:10.1016/0167-6636(87)90005-6

CrossRef Full Text | Google Scholar

Benveniste, Y., Dvorak, G. J., and Chen, T. (1991). On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media. J. Mech. Phys. Solids 39, 927–946. doi:10.1016/0022-5096(91)90012-d

CrossRef Full Text | Google Scholar

Bernard, O., Ulm, F.-J., and Lemarchand, E. (2003). A multiscale micromechanics-hydration model for the early-age elastic properties of cement-based materials. Cem. Concr. Res. 33, 1293–1309. doi:10.1016/s0008-8846(03)00039-5

CrossRef Full Text | Google Scholar

Bertrand, E., and Hellmich, C. (2009). Multiscale elasticity of tissue engineering scaffolds with tissue-engineered bone: A continuum micromechanics approach. J. Eng. Mech. 135, 395–412. doi:10.1061/(asce)0733-9399(2009)135:5(395)

CrossRef Full Text | Google Scholar

Brisard, S., and Dormieux, L. (2010). FFT-based methods for the mechanics of composites: A general variational framework. Comput. Mater. Sci. 49, 663–671. doi:10.1016/j.commatsci.2010.06.009

CrossRef Full Text | Google Scholar

Buchner, T., Königsberger, M., Jäger, A., and Füssl, J. (2022). A validated multiscale model linking microstructural features of fired clay brick to its macroscopic multiaxial strength. Mech. Mater. 170, 104334. doi:10.1016/j.mechmat.2022.104334

CrossRef Full Text | Google Scholar

Buryachenko, V. (2007). Micromechanics of heterogeneous materials. Springer Science and Business Media.

Google Scholar

Cai, X., Brenner, R., Peralta, L., Olivier, C., Gouttenoire, P.-J., Chappard, C., et al. (2019). Homogenization of cortical bone reveals that the organization and shape of pores marginally affect elasticity. J. R. Soc. Interface 16, 20180911. doi:10.1098/rsif.2018.0911

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, M. K., Dalton, K. M., Shen, L., Evans, A. C., and Davidson, R. J. (2007). Weighted Fourier series representation and its application to quantifying the amount of gray matter. IEEE Trans. Med. Imaging 26, 566–581. doi:10.1109/tmi.2007.892519

PubMed Abstract | CrossRef Full Text | Google Scholar

Cochran, W., Cooley, J., Favin, D., Helms, H., Kaenel, R., Lang, W., et al. (1967). What is the fast Fourier transform? Proc. IEEE 55, 1664–1674. doi:10.1109/PROC.1967.5957

CrossRef Full Text | Google Scholar

Dvorak, G. J. (2012). Micromechanics of composite materials. Springer Science and Business Media.

Google Scholar

Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 241, 376–396.

Google Scholar

Feng, T., Jia, M., Xu, W., Wang, F., Li, P., Wang, X., et al. (2022). Three-dimensional mesoscopic investigation of the compression mechanical properties of ultra-high performance concrete containing coarse aggregates. Cem. Concr. Compos. 133, 104678. doi:10.1016/j.cemconcomp.2022.104678

CrossRef Full Text | Google Scholar

Fredholm, I. (1900). Sur les équations de l’équilibre d’un corps solide élastique. Acta Math. 23, 1–42. doi:10.1007/bf02418668

CrossRef Full Text | Google Scholar

Fritsch, A., Dormieux, L., and Hellmich, C. (2006). Porous polycrystals built up by uniformly and axisymmetrically oriented needles: Homogenization of elastic properties. Comptes Rendus Mécanique 334, 151–157. doi:10.1016/j.crme.2006.01.008

CrossRef Full Text | Google Scholar

Fritsch, A., Dormieux, L., Hellmich, C., and Sanahuja, J. (2009a). Mechanical behavior of hydroxyapatite biomaterials: An experimentally validated micromechanical model for elasticity and strength. J. Biomed. Mater. Res. Part A 88A, 149–161. doi:10.1002/jbm.a.31727

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritsch, A., Hellmich, C., and Dormieux, L. (2009b). Ductile sliding between mineral crystals followed by rupture of collagen crosslinks: Experimentally supported micromechanical explanation of bone strength. J. Theor. Biol. 260, 230–252. doi:10.1016/j.jtbi.2009.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritsch, A., Hellmich, C., and Young, P. (2013). Micromechanics-derived scaling relations for poroelasticity and strength of brittle porous polycrystals. J. Appl. Mech. 80, 020905. doi:10.1115/1.4007922

CrossRef Full Text | Google Scholar

Fritsch, A., and Hellmich, C. (2007). ‘Universal’ microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: Micromechanics-based prediction of anisotropic elasticity. J. Theor. Biol. 244, 597–620. doi:10.1016/j.jtbi.2006.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimal, Q., Raum, K., Gerisch, A., and Laugier, P. (2011). A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties. Biomechanics Model. Mechanobiol. 10, 925–937. doi:10.1007/s10237-010-0284-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, W., Han, F., Jiang, J., and Xu, W. (2022). A micromechanical framework for thermo-elastic properties of multiphase cementitious composites with different saturation. Int. J. Mech. Sci. 224, 107313. doi:10.1016/j.ijmecsci.2022.107313

CrossRef Full Text | Google Scholar

Hashin, Z. (1983). Analysis of composite materials—a survey. J. Appl. Mech. 50, 481–505. doi:10.1115/1.3167081

CrossRef Full Text | Google Scholar

Hashin, Z. (1963). Theory of mechanical behavior of heterogeneous media. Philidelphia, Pa: Towne School of Civil and Mechanical Engineering, University of Pennsylvania.

Google Scholar

Hashin, Z. (1965). Viscoelastic behavior of heterogeneous media. J. Appl. Mech. 32, 630–636. doi:10.1115/1.3627270

CrossRef Full Text | Google Scholar

Hellmich, C., Barthélémy, J.-F., and Dormieux, L. (2004). Mineral–collagen interactions in elasticity of bone ultrastructure – A continuum micromechanics approach. Eur. J. Mech. - A/Solids 23, 783–810. doi:10.1016/j.euromechsol.2004.05.004

CrossRef Full Text | Google Scholar

Hellmich, C., and Mang, H. (2005). Shotcrete elasticity revisited in the framework of continuum micromechanics: From submicron to meter level. J. Mater. Civ. Eng. 17, 246–256. doi:10.1061/(asce)0899-1561(2005)17:3(246)

CrossRef Full Text | Google Scholar

Hervé, E., and Zaoui, A. (1993). inclusion-based micromechanical modelling. Int. J. Eng. Sci. 31, 1–10. doi:10.1016/0020-7225(93)90059-4

CrossRef Full Text | Google Scholar

Hill, R. (1963). Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids 11, 357–372. doi:10.1016/0022-5096(63)90036-x

CrossRef Full Text | Google Scholar

Hofstetter, K., Hellmich, C., and Eberhardsteiner, J. (2005). Development and experimental validation of a continuum micromechanics model for the elasticity of wood. Eur. J. Mech. - A/Solids 24, 1030–1053. doi:10.1016/j.euromechsol.2005.05.006

CrossRef Full Text | Google Scholar

Horwitz, A. (2001). A version of Simpson’s rule for multiple integrals. J. Comput. Appl. Math. 134, 1–11. doi:10.1016/s0377-0427(00)00444-1

CrossRef Full Text | Google Scholar

Jiménez Segura, N., Pichler, B. L., and Hellmich, C. (2023). Concentration tensors preserving elastic symmetry of multiphase composites. Mech. Mater. 178, 104555. doi:10.1016/j.mechmat.2023.104555

CrossRef Full Text | Google Scholar

Jiménez Segura, N., Pichler, B. L., and Hellmich, C. (2022). Stress average rule derived through the principle of virtual power. ZAMM - J. Appl. Math. Mech./Zeitschrift für Angewandte Math. und Mech. 102, e202200091. doi:10.1002/zamm.202200091

CrossRef Full Text | Google Scholar

Kneer, G. (1965). Über die Berechnung der Elastizitätsmoduln vielkristalliner Aggregate mit Textur. Phys. Status Solidi (b) 9, 825–838. doi:10.1002/pssb.19650090319

CrossRef Full Text | Google Scholar

Königsberger, M., Hlobil, M., Delsaute, B., Staquet, S., Hellmich, C., and Pichler, B. (2018). Hydrate failure in itz governs concrete strength: A micro-to-macro validated engineering mechanics model. Cem. Concr. Res. 103, 77–94. doi:10.1016/j.cemconres.2017.10.002

CrossRef Full Text | Google Scholar

Korringa, J. (1973). Theory of elastic constants of heterogeneous media. J. Math. Phys. 14, 509–513. doi:10.1063/1.1666346

CrossRef Full Text | Google Scholar

Kröner, E. (1958). Berechnung der elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls [Calculation of the elastic constant of the multi-crystal from the constants of the single crystals]. Z. für Phys. 151, 504–518.

Google Scholar

Kröner, E. (1977). Bounds for effective elastic moduli of disordered materials. J. Mech. Phys. Solids 25, 137–155. doi:10.1016/0022-5096(77)90009-6

CrossRef Full Text | Google Scholar

Levin, V. M. (1967). Thermal expansion coefficient of heterogeneous materials. Mekhanika Tverd. Tela 2, 83–94.

Google Scholar

Lipinski, P., Barhdadi, E. H., and Cherkaoui, M. (2006). Micromechanical modelling of an arbitrary ellipsoidal multi-coated inclusion. Philos. Mag. 86, 1305–1326. doi:10.1080/14786430500343868

CrossRef Full Text | Google Scholar

Lippmann, B. A., and Schwinger, J. (1950). Variational principles for scattering processes. I. Phys. Rev. 79, 469–480. doi:10.1103/physrev.79.469

CrossRef Full Text | Google Scholar

Lucarini, S., Upadhyay, M., and Segurado, J. (2021). Fft based approaches in micromechanics: Fundamentals, methods and applications. Model. Simul. Mater. Sci. Eng. 30, 023002. doi:10.1088/1361-651x/ac34e1

CrossRef Full Text | Google Scholar

Moës, N., Cloirec, M., Cartraud, P., and Remacle, J. F. (2003). A computational approach to handle complex microstructure geometries. Comput. Methods Appl. Mech. Eng. 192, 3163–3177. doi:10.1016/s0045-7825(03)00346-3

CrossRef Full Text | Google Scholar

Molinari, A., and El Mouden, M. (1996). The problem of elastic inclusions at finite concentration. Int. J. Solids Struct. 33, 3131–3150. doi:10.1016/0020-7683(95)00275-8

CrossRef Full Text | Google Scholar

Mori, T., and Tanaka, K. (1973). Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 21, 571–574. doi:10.1016/0001-6160(73)90064-3

CrossRef Full Text | Google Scholar

Moulinec, H., and Suquet, P. (1994). A fast numerical method for computing the linear and nonlinear properties of composites. Comptes-Rendus l’Académie Sci. Série II 318, 1417–1423.

Google Scholar

Moulinec, H., and Suquet, P. (1998). A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Methods Appl. Mech. Eng. 157, 69–94. doi:10.1016/s0045-7825(97)00218-1

CrossRef Full Text | Google Scholar

Pahr, D. H., and Zysset, P. K. (2008). Influence of boundary conditions on computed apparent elastic properties of cancellous bone. Biomechanics Model. Mechanobiol. 7, 463–476. doi:10.1007/s10237-007-0109-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Pichler, B., Hellmich, C., Eberhardsteiner, J., Wasserbauer, J., Termkhajornkit, P., Barbarulo, R., et al. (2013). Effect of gel–space ratio and microstructure on strength of hydrating cementitious materials: An engineering micromechanics approach. Cem. Concr. Res. 45, 55–68. doi:10.1016/j.cemconres.2012.10.019

CrossRef Full Text | Google Scholar

Pichler, B., and Hellmich, C. (2011). Upscaling quasi-brittle strength of cement paste and mortar: A multi-scale engineering mechanics model. Cem. Concr. Res. 41, 467–476. doi:10.1016/j.cemconres.2011.01.010

CrossRef Full Text | Google Scholar

Rosen, B. W., and Hashin, Z. (1970). Effective thermal expansion coefficients and specific heats of composite materials. Int. J. Eng. Sci. 8, 157–173. doi:10.1016/0020-7225(70)90066-2

CrossRef Full Text | Google Scholar

Sanahuja, J., Dormieux, L., Meille, S., Hellmich, C., and Fritsch, A. (2010). Micromechanical explanation of elasticity and strength of gypsum: From elongated anisotropic crystals to isotropic porous polycrystals. J. Eng. Mech. 136, 239–253. doi:10.1061/(asce)em.1943-7889.0000072

CrossRef Full Text | Google Scholar

Scheiner, S., Sinibaldi, R., Pichler, B., Komlev, V., Renghini, C., Vitale-Brovarone, C., et al. (2009). Micromechanics of bone tissue-engineering scaffolds, based on resolution error-cleared computer tomography. Biomaterials 30, 2411–2419. doi:10.1016/j.biomaterials.2008.12.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Sevostianov, I., and Kachanov, M. (2014). On some controversial issues in effective field approaches to the problem of the overall elastic properties. Mech. Mater. 69, 93–105. doi:10.1016/j.mechmat.2013.09.010

CrossRef Full Text | Google Scholar

Ting, T. C. T., and Lee, V.-G. (1997). The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids. Q. J. Mech. Appl. Math. 50, 407–426. doi:10.1093/qjmam/50.3.407

CrossRef Full Text | Google Scholar

Tonon, F., Pan, E., and Amadei, B. (2001). Green’s functions and boundary element method formulation for 3D anisotropic media. Comput. Struct. 79, 469–482. doi:10.1016/s0045-7949(00)00163-2

CrossRef Full Text | Google Scholar

Torquato, S. (1997). Effective stiffness tensor of composite media—I. Exact series expansions. J. Mech. Phys. Solids 45, 1421–1448. doi:10.1016/s0022-5096(97)00019-7

CrossRef Full Text | Google Scholar

Voigt, W. (1889). Über die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper [On the relation between the elasticity constants of isotropic bodies]. Ann. Phys. 274, 573–587. doi:10.1002/andp.18892741206

CrossRef Full Text | Google Scholar

Wang, H., Hellmich, C., Yuan, Y., Mang, H., and Pichler, B. (2018). May reversible water uptake/release by hydrates explain the thermal expansion of cement paste? — Arguments from an inverse multiscale analysis. Cem. Concr. Res. 113, 13–26. doi:10.1016/j.cemconres.2018.05.008

CrossRef Full Text | Google Scholar

Whittaker, E., and Robinson, G. (1967). “The trapezoidal and parabolic rules,”in The calculus of observations: A treatise on numerical mathematics (Dover, New York: Cope Press), 156–158.

Google Scholar

Willis, J. R. (1977). Bounds and self-consistent estimates for the overall properties of anisotropic composites. J. Mech. Phys. Solids 25, 185–202. doi:10.1016/0022-5096(77)90022-9

CrossRef Full Text | Google Scholar

Wolfram, U., Peña Fernández, M., McPhee, S., Smith, E., Beck, R. J., Shephard, J. D., et al. (2022). Multiscale mechanical consequences of ocean acidification for cold-water corals. Sci. Rep. 12, 8052. doi:10.1038/s41598-022-11266-w

CrossRef Full Text | Google Scholar

Xie, L., Zhang, C., Sladek, J., and Sladek, V. (2016). Unified analytical expressions of the three-dimensional fundamental solutions and their derivatives for linear elastic anisotropic materials. Proc. R. Soc. A Math. Phys. Eng. Sci. 472, 20150272. doi:10.1098/rspa.2015.0272

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W., Wu, F., Jiao, Y., and Liu, M. (2017). A general micromechanical framework of effective moduli for the design of nonspherical nano- and micro-particle reinforced composites with interface properties. Mater. Des. 127, 162–172. doi:10.1016/j.matdes.2017.04.075

CrossRef Full Text | Google Scholar

Xu, W., Wu, Y., and Gou, X. (2019). Effective elastic moduli of nonspherical particle-reinforced composites with inhomogeneous interphase considering graded evolutions of elastic modulus and porosity. Comput. Methods Appl. Mech. Eng. 350, 535–553. doi:10.1016/j.cma.2019.03.021

CrossRef Full Text | Google Scholar

Xu, W., Wu, Y., and Jia, M. (2018). Elastic dependence of particle-reinforced composites on anisotropic particle geometries and reinforced/weak interphase microstructures at nano- and micro-scales. Compos. Struct. 203, 124–131. doi:10.1016/j.compstruct.2018.07.009

CrossRef Full Text | Google Scholar

Zaoui, A. (2002). Continuum micromechanics: Survey. J. Eng. Mech. 128, 808–816. doi:10.1061/(asce)0733-9399(2002)128:8(808)

CrossRef Full Text | Google Scholar

Zeller, R., and Dederichs, P. (1973). Elastic constants of polycrystals. Phys. Status Solidi (b) 55, 831–842. doi:10.1002/pssb.2220550241

CrossRef Full Text | Google Scholar

Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2005). The finite element method: Its basis and fundamentals. Elsevier.

Google Scholar

Keywords: complex microstructure, Green’s function, concentration tensor, homogenized stiffness, Fredholm integral

Citation: Jiménez Segura N, Pichler BLA and Hellmich C (2023) A Green’s function-based approach to the concentration tensor fields in arbitrary elastic microstructures. Front. Mater. 10:1137057. doi: 10.3389/fmats.2023.1137057

Received: 03 January 2023; Accepted: 13 March 2023;
Published: 21 April 2023.

Edited by:

Marek Pawlikowski, Warsaw University of Technology, Poland

Reviewed by:

Wenxiang Xu, Hohai University, China
Reinaldo Rodriguez-Ramos, University of Havana, Cuba

Copyright © 2023 Jiménez Segura, Pichler and Hellmich. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Christian Hellmich, christian.hellmich@tuwien.ac.at

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.