Front. Mater. Frontiers in Materials Front. Mater. 2296-8016 Frontiers Media S.A. 1137057 10.3389/fmats.2023.1137057 Materials Original Research A Green’s function-based approach to the concentration tensor fields in arbitrary elastic microstructures Jiménez Segura et al. 10.3389/fmats.2023.1137057 Jiménez Segura Nabor 1 2 Pichler Bernhard L. A. 1 Hellmich Christian 1 * 1 Institute for Mechanics of Materials and Structures (IMWS), TU Wien (Vienna University of Technology), Vienna, Austria 2 Department of Materials Science, Polytechnic University of Madrid (UPM), Madrid, Spain

Edited by: Marek Pawlikowski, Warsaw University of Technology, Poland

Reviewed by: Wenxiang Xu, Hohai University, China

Reinaldo Rodriguez-Ramos, University of Havana, Cuba

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

This article was submitted to Computational Materials Science, a section of the journal Frontiers in Materials

 21 04 2023 2023 10 1137057 03 01 2023 13 03 2023 Copyright © 2023 Jiménez Segura, Pichler and Hellmich. 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.

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.

 complex microstructure Green’s function concentration tensor homogenized stiffness Fredholm integral 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, S RVE, in the form proposed by Hashin (1983)

u ̲ x ̲ = E ⋅ x ̲ , ∀ x ̲ ∈ S RVE ,

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 = 1 V RVE ∫ V RVE ε x ̲ d V x ̲ = ⟨ ε ⟩ ,

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 ̲ = 1 2 grad x u ̲ x ̲ + grad x u ̲ x ̲ T = grad x S u ̲ x ̲ ; ∀ x ̲ ∈ V RVE ,

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

σ x ̲ = c x ̲ : ε x ̲ , ∀ x ̲ ∈ V RVE ,

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

• equilibrium conditions

 ∇ ̲ x ⋅ σ x ̲ + f ̲ x ̲ = 0 , ∀ x ̲ ∈ V RVE ,

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)

 Σ = 1 V RVE ∫ V RVE σ x ̲ d V x ̲ = ⟨ σ ⟩ ;

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

 ε x ̲ = A x ̲ : E , ∀ x ̲ ∈ V RVE ,

where A is the concentration (or downscaling) tensor;

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

 Σ = C hom : E ,

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

 C hom = 1 V RVE ∫ V RVE c x ̲ : A x ̲ d V x ̲ . The classical approach for making the problem of Eqs. 1–9 tractable is to restrict the discussion to N r homogeneous subdomains or phases within the RVE. Accordingly, the general microstiffness distribution c ( x ̲ ) is replaced by a finite number of microstiffness tensors c r , r = 1 , … , N r , 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).

Mathematical symbols and abbreviations.

Variables
1 = second-order identity tensor
A = strain concentration tensor
A 0 = strain concentration tensor field associated with an arbitrary inhomogeneous microstructure filling an infinite medium
A n 0 = nth term of the series defining the auxiliary strain concentration tensor field in the infinite elastic domain
A 0 , I = Eshelby-problem-related strain concentration tensor associated with inhomogeneity I embedded in an infinite elastic domain
c = microscopic stiffness tensor
C hom = macroscopic homogenized stiffness tensor
c ijkℓ = component ijkℓ of c
C 0 = stiffness tensor of homogeneous elastic space
CPU = Central Processing Unit
E = macroscopic strain subjected to RVE
e ̲ i = orthonormal base vector
E 0 = auxiliary strain subjected to infinite homogeneous elastic matrix
f ̲ = microscopic volume force
FEM = Finite Element Method
FFT = Fast Fourier Transform
G = Green’s function
G = third-order tensor denoting the symmetric gradient of Green’s function
G = fourth-order tensor denoting the symmetric gradient of the gradient of Green’s function
G ijkℓ = component ijkℓ of G
I = symmetric fourth-order identity tensor
I ijkℓ vol = component ijkℓ of I vol
k 0 = bulk modulus of homogeneous elastic space
M = fourth-order RVE-to-remote strain conversion tensor
R = residual term in solution for implicit integral
n = integer numbering member in series expression for A 0
n λ = number of stiffness waves along edge directions of a box-shaped RVE with sinusoidal microstructure
RVE = representative volume element
S RVE = surface of the RVE
u ̲ = microscopic displacement field
u ̲ h = homogeneous solution to differential equation for u ̲
u ̲ p = particular solution to differential equation for u ̲
V RVE = volume of the RVE
x ̲ = microscopic position vector
x i = ith component of x ̲
y ̲ = microscopic position vector in convolution integral formulation
z ̲ = microscopic position vector in convolution integral formulation
α = integer labeling iteration step in solution of Fredholm integral equation
γ = variable in geometric series
δ = Dirac delta function
δ ij = Kronecker delta
Δk = microscopic stiffness fluctuation
ɛ = microscopic strain tensor field
θ = polar coordinate
λ = wavelength of microscopic stiffness fluctuation
ν 0 = Poisson’s ratio of homogeneous elastic space
ρ = radial coordinate
σ = microscopic stress tensor
Σ = macroscopic stress associated with the RVE
τ = polarization stress
ϕ = solution of Poisson’s equation
Operators
grad x = microscopic gradient with respect to variable x ̲
grad x S = microscopic symmetric gradient with respect to variable x ̲
• T = transpose operator, acting on second-order tensor as • i j T = • j i and on fourth-order tensor as • ijkℓ T = • k ℓ i j
∏ r = 1 3 • r = double contraction-type product of a bunch of tensors: •1:•2:•3
∑ r = 1 N • r = summation over variable r from 1 to N
∇ ̲ x ⋅ • = microscopic divergence with respect to variable x ̲
⋅ = dot product or contraction product
: = double contraction product
⊗ = dyadic product
⟨•⟩ = volume average over the RVE, ⟨ • ⟩ = 1 V RVE ∫ V RVE • ( x ̲ ) d V ( x ̲ )

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. 1–9, we observe that both FEM and FFT-based homogenization techniques primarily focus on the homogenized stiffness tensor C hom , 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 E 0. The solution of Eshelby (1957) then relates the auxiliary strains to the homogeneous strains ɛ I in the inhomogeneity, according to ε I = A 0 , I : E 0 , with A 0 , I as the concentration tensor associated with inhomogeneity I embedded in an infinite matrix subjected to E 0. A 0 , 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 ̲ = A 0 x ̲ : E 0 , with A 0 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 C 0 . The strains ε ( x ̲ ) then consist of two portions: (i) auxiliary strains E 0 prevailing in the homogeneous infinite domain of stiffness C 0 , see Figure 1A, and (ii) fluctuations around E 0 which arise from the fluctuations in the stiffness field, [ c ( x ̲ ) − C 0 ] , 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 = 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ : E 0 = M − 1 : E 0 ⇒ M = 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ − 1 , 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 ̲ = A 0 x ̲ : 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ − 1 : E , and comparison of Eq. 13 with Eq. 7 yields the strain concentration tensor as A x ̲ = A 0 x ̲ : 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ − 1 = A 0 x ̲ : M , and when considering, in addition, Eq. 9, the homogenized stiffness tensor is eventually retrieved as C hom = 1 V RVE ∫ V RVE c x ̲ : A 0 x ̲ d V x ̲ : 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ − 1 , see Figures 1E, F. Eq. 15 can be reformulated by setting c ( x ̲ ) = C 0 + c ( x ̲ ) − C 0 , which yields C hom = C 0 + 1 V RVE ∫ V RVE c x ̲ − C 0 : A 0 x ̲ d V x ̲ : 1 V RVE ∫ V RVE A 0 x ̲ d V x ̲ − 1 . We are left with the determination of A 0 ( 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).

Illustration of elements making up the new homogenization scheme: (A) infinite homogeneous elastic matrix of stiffness C 0 subjected to background (auxiliary) strains E 0; (B) same matrix undergoing equilibrated (polarization) stresses which are equivalent to the effect of fluctuations ( c ( x ̲ ) − C 0 ) in microstiffness, c ( x ̲ ) , around the homogeneous stiffness C 0 , (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 ⟨ • ⟩ = 1 V RVE ∫ V RVE • ( x ̲ ) d V ( 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 A 0 ( 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 ̲ = C 0 : ε x ̲ + τ x ̲ . Equating Eq. 17 with Eq. 4 yields τ x ̲ = c x ̲ − C 0 : ε x ̲ , 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 ∇ ̲ x ⋅ C 0 : grad x S u ̲ x ̲ = − f ̲ x ̲ + ∇ ̲ x ⋅ τ x ̲ .

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 ̲ h x ̲ + u ̲ p x ̲ , whereby

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

∇ ̲ x ⋅ C 0 : grad S u ̲ h x ̲ = 0 ,

with the inhomogeneous boundary conditions u ̲ h x ̲ = E 0 X ̲ ⋅ x ̲ , ∀ | x ̲ | → ∞ ;

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

∇ ̲ x ⋅ C 0 : grad S u ̲ p x ̲ = − f ̲ x ̲ + ∇ ̲ x ⋅ τ x ̲ ,

with homogeneous boundary conditions u ̲ p x ̲ = 0 , ∀ | x ̲ | → ∞ . The homogeneous solution reads as u ̲ h x ̲ = E 0 X ̲ ⋅ x ̲ , ∀ x ̲ ∈ R 3 . The particular solution can be given in the form u ̲ p x ̲ = ∫ R 3 G x ̲ − y ̲ ⋅ f ̲ y ̲ + ∇ ̲ y ⋅ τ y ̲ d V y ̲ , 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) ∇ ̲ x ⋅ C 0 : grad x S G x ̲ − y ̲ = − 1 δ x ̲ − y ̲ , with boundary conditions G x ̲ − y ̲ = 0 , ∀ | x ̲ | → ∞ . In Eq. 27, 1 denotes the second-order identity tensor and δ denotes the Dirac function, with the following properties. δ x ̲ = 0 for x ̲ ≠ 0 , δ x ̲ = ∞ for x ̲ = 0 , ∫ R 3 δ x ̲ d V x ̲ = 1 . 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 ∫ R 3 G x ̲ − y ̲ ⋅ ∇ ̲ y ⋅ τ y ̲ d V y ̲ = − ∫ R 3 grad y G x ̲ − y ̲ : τ y ̲ d V y ̲ . 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 ̲ = E 0 + ∫ R 3 G x ̲ − y ̲ ⋅ f ̲ y ̲ d V y ̲ − ∫ R 3 G x ̲ − y ̲ : τ y ̲ d V y ̲ . In Eq. 33, the following gradients of the Green’s functions were explicitly introduced G x ̲ − y ̲ = grad x S G x ̲ − y ̲ , G x ̲ − y ̲ = grad x S grad y G x ̲ − y ̲ . 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 ̲ = E 0 − ∫ R 3 G x ̲ − y ̲ : c y ̲ − C 0 : ε y ̲ d V y ̲ . 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 ̲ = E 0 − ∫ R 3 G x ̲ − y ̲ 1 : c y ̲ 1 − C 0 : ε y ̲ 1 d V y ̲ 1 . 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 = E 0 − ∫ R 3 G y ̲ 1 − y ̲ 2 : c y ̲ 2 − C 0 : ε y ̲ 2 d V y ̲ 2 . Insertion of Eq. 38 into Eq. 37 yields ε x ̲ = I − ∫ R 3 G x ̲ − y ̲ 1 : c y ̲ 1 − C 0 d V y ̲ 1 : E 0 + ∫ R 3 G x ̲ − y ̲ 1 : c y ̲ 1 − C 0 : ∫ R 3 G y ̲ 1 − y ̲ 2 : c y ̲ 2 − C 0 : ε y ̲ 2 d V y ̲ 2 d V y ̲ 1 ,

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 A 0 ( x ̲ ) and a residual term comprising the implicit strain field R [ x ̲ , ε ( y ̲ ( N ) ) ] , reading as ε x ̲ = A 0 x ̲ : E 0 + R x ̲ , ε y ̲ N . In more detail, the concentration tensor field reads as A 0 x ̲ = I + ∑ n = 1 N − 1 n A n 0 x ̲ , where A 1 0 x ̲ = ∫ R 3 G x ̲ − y ̲ : c y ̲ − C 0 d V y ̲ , and, for n > 1, A n 0 x ̲ = ∫ R 3 … ∫ R 3 G x ̲ − y ̲ 1 : ∏ i = 2 n c y ̲ i − 1 − C 0 : G y ̲ i − 1 − y ̲ i : c y ̲ n − C 0 d V y ̲ n … d V y ̲ 1 , ∀ n ∈ 2 , ∞ .

The residual term after N iterations reads as R x ̲ , ε y ̲ N = ∫ R 3 … ∫ R 3 G x ̲ − y ̲ 1 : ∏ i = 2 N c y ̲ i − 1 − C 0 : G y ̲ i − 1 − y ̲ i : c y ̲ N − C 0 : ε y ̲ N × d V y ̲ N … d V y ̲ 1 .

Note that for N → ∞, R [ x ̲ , ε ( y ̲ ( N ) ) ] → 0 , provided the following requirement is met: ∫ R 3 G ijkℓ y ̲ α − 1 − y ̲ α c k ℓ m n y ̲ α − C k ℓ m n 0 × G mnpq y ̲ α − y ̲ α + 1 d V y ̲ α < 1 , ∀ y ̲ α − 1 , y ̲ α + 1 ∈ V RVE , 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 A 0 , 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, V RVE, (while every other y ̲ ( i ) is integrated over the entire space, R 3 ), allows for the following expression M = I + ∑ n = 1 N − 1 n V RVE ∫ V RVE ∫ R 3 … ∫ R 3 ∏ i = 1 n G y ̲ i − 1 − y ̲ i : c y ̲ i − C 0 × d V y ̲ n … d V y ̲ 0 − 1 . Thus, a comprehenive integral format for the strain concentration tensor field is obtained from inserting Eqs. 41–43 and Eq. 46 into Eq. 14, yielding A x ̲ = I − ∫ R 3 G x ̲ − y ̲ : c y ̲ − C 0 d V y ̲ + ∑ n = 2 N − 1 n ∫ R 3 … ∫ R 3 G x ̲ − y ̲ 1 : ∏ i = 2 n c y ̲ i − 1 − C 0 : G y ̲ i − 1 − y ̲ i : c y ̲ n − C 0 d V y ̲ n … d V y ̲ 1 : I + ∑ n = 1 N − 1 n V RVE ∫ V RVE ∫ R 3 … ∫ R 3 ∏ i = 1 n G y ̲ i − 1 − y ̲ i : c y ̲ i − C 0 × d V y ̲ n … d V y ̲ 0 − 1 .

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 C hom = C 0 + 1 V RVE ∫ V RVE c y ̲ 0 − C 0 d V y ̲ 0 + ∑ n = 1 N − 1 n V RVE ∫ V RVE ∫ R 3 … ∫ R 3 ∏ i = 1 n c y ̲ i − 1 − C 0 : G y ̲ i − 1 − y ̲ i : c y ̲ n − C 0 d V y ̲ n … d V y ̲ 0 : I + ∑ n = 1 N − 1 n V RVE ∫ V RVE ∫ R 3 … ∫ R 3 ∏ i = 1 n G y ̲ i − 1 − y ̲ i : c y ̲ i − C 0 × d V y ̲ n … d V y ̲ 0 − 1 .

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. 41–43, we resort to the Green’s function for an infinitely extended isotropic elastic body with bulk modulus k 0 and Poisson’s ratio ν 0; it reads as (Dvorak, 2012) G x ̲ − y ̲ = 1 + ν 0 6 π k 0 1 − 2 ν 0 1 1 | x ̲ − y ̲ | − 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 grad x grad x | x ̲ − y ̲ | , where, with respect to the actual formula (4.5.12) given on page 107 of (Dvorak, 2012), we consider k 0 = 2 μ 0 ( 1 + ν 0 ) 3 ( 1 − 2 ν 0 ) , with μ 0 being the shear modulus. Moreover, we consider a sinusoidal stiffness distribution across this infinitely extended body, according to c x ̲ = C 0 + 3 Δ k sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 I vol , where x 1, x 2, and x 3 are the components of location vector x ̲ with respect to an orthonormal base frame e ̲ 1 , e ̲ 2 , e ̲ 3 , such that x ̲ = x 1 e ̲ 1 + x 2 e ̲ 2 + x 3 e ̲ 3 , and I vol stands for the volumetric part of the symmetric fourth-order unity tensor I . The latter has the components I ijrs = 1/2 (δ ir δ js + δ is δ jr ), while the former reads as I vol = 1 / 3 ( 1 ⊗ 1 ) , with components I ijkℓ vol = 1 / 3 ( δ i j δ 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 <inline-formula id="inf74"> <mml:math id="m124"> <mml:msup> <mml:mrow> <mml:mi mathvariant="double-struck">A</mml:mi> </mml:mrow> <mml:mrow> <mml:mn>0</mml:mn> </mml:mrow> </mml:msup> <mml:mrow> <mml:mo stretchy="false">(</mml:mo> <mml:mrow> <mml:munder accentunder="false"> <mml:mrow> <mml:mi>x</mml:mi> </mml:mrow> <mml:mo accent="true">̲</mml:mo> </mml:munder> </mml:mrow> <mml:mo stretchy="false">)</mml:mo> </mml:mrow> </mml:math> </inline-formula>

The auxiliary strain downscaling tensor A 0 ( x ̲ ) is the result of the sum of an infinite series, see Eq. 41. The first term of this series, A 1 0 ( x ̲ ) , is defined through Eq. 42, so that consideration of Eq. 50 yields the component 1111 of A 1 0 as A 1,1111 0 x ̲ = ∑ h = 1 3 ∑ ℓ = 1 3 ∫ R 3 G 11 hℓ hℓ x ̲ − y ̲ 3 Δ k sin 2 π λ 1 y 1 × sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 I h ℓ 11 vol d V y ̲ = Δ k ∑ ℓ = 1 3 ∫ R 3 G 11 ℓ ℓ x ̲ − y ̲ sin 2 π λ 1 y 1 × sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 d V y ̲ . In order to retrieve the components G 1111, G 1122, and G 1133, we start with the general expression for the components of the Green’s function of Eq. 49, reading as G i j = 1 + ν 0 6 π k 0 1 − 2 ν 0 δ i j | x ̲ − y ̲ | − 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 2 ∂ x i ∂ x j | x ̲ − y ̲ | . The fourth-order Green’s function gradient according to Eq. 35 exhibits the following components G ijkℓ = 1 + ν 0 12 π k 0 1 − 2 ν 0 ∂ 2 ∂ x j ∂ y ℓ δ i k | x ̲ − y ̲ | + ∂ 2 ∂ x i ∂ y ℓ δ j k | x ̲ − y ̲ | − 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 4 ∂ x i ∂ x j ∂ x k ∂ y ℓ | x ̲ − y ̲ | .

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 ∂ ∂ x i | x ̲ − y ̲ | = x i − y i | x ̲ − y ̲ | , revealing the interesting property ∂ ∂ x i | x ̲ − y ̲ | = − ∂ ∂ y i | x ̲ − y ̲ | . The first derivative of the inverse of the norm  | x ̲ − y ̲ | reads as ∂ ∂ x i 1 | x ̲ − y ̲ | = − x i − y i | x ̲ − y ̲ | 3 , revealing the interesting property ∂ ∂ x i 1 | x ̲ − y ̲ | = − ∂ ∂ y i 1 | x ̲ − y ̲ | . Eq. 55 and Eq. 57 allow for re-writing Eq. 53 as G ijkℓ = − 1 + ν 0 12 π k 0 1 − 2 ν 0 ∂ 2 ∂ x j ∂ x ℓ δ i k | x ̲ − y ̲ | + ∂ 2 ∂ x i ∂ x ℓ δ j k | x ̲ − y ̲ | + 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 4 ∂ x i ∂ x j ∂ x k ∂ x ℓ | x ̲ − y ̲ | ,

so that the components occurring in Eq. 51 read as G 1111 = − 1 + ν 0 6 π k 0 1 − 2 ν 0 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | + 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 4 ∂ x 1 4 | x ̲ − y ̲ | , G 1122 = 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 4 ∂ x 1 2 ∂ x 2 2 | x ̲ − y ̲ | , G 1133 = 1 + ν 0 24 π k 0 1 − ν 0 1 − 2 ν 0 ∂ 4 ∂ x 1 2 ∂ x 3 2 | x ̲ − y ̲ | . The sum of Eqs. 59–61, 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, ∂ 2 ∂ x i 2 | x ̲ − y ̲ | = ∑ k ≠ i x k − y k 2 | x ̲ − y ̲ | 3 . Namely, Eq. 62 entails the following identity ∂ 2 ∂ x 1 2 | x ̲ − y ̲ | + ∂ 2 ∂ x 2 2 | x ̲ − y ̲ | + ∂ 2 ∂ x 3 2 | x ̲ − y ̲ | = 2 1 | x ̲ − y ̲ | , twofold derivation of which yields an identity comprising the derivatives of the norm  | x ̲ − y ̲ | occurring in Eqs. 59–61, reading as ∂ 4 ∂ x 1 4 | x ̲ − y ̲ | + ∂ 4 ∂ x 1 2 ∂ x 2 2 | x ̲ − y ̲ | + ∂ 4 ∂ x 1 2 ∂ x 3 2 | x ̲ − y ̲ | = 2 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | . Insertion of Eqs. 59–61 into Eq. 51, while considering Eq. 64, yields A 1,1111 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 ∫ R 3 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | sin 2 π λ 1 y 1 × sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 d V y ̲ . A change of variable to z ̲ = y ̲ − x ̲ results in A 1,1111 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 ∫ R 3 ∂ 2 ∂ z 1 2 1 | − z ̲ | × sin 2 π λ 1 x 1 + z 1 sin 2 π λ 2 x 2 + z 2 × sin 2 π λ 3 x 3 + z 3 d V z ̲ .

Transforming Eq. 66 by means of sin a + b = sin a cos b + cos a sin b , and considering the second derivative of the inverse of the norm of − z ̲ as ∂ 2 ∂ z 1 2 1 | − z ̲ | = 2 − z 1 2 − − z 2 2 − − z 3 2 | − z ̲ | 5 , as well as that the even and odd functions appearing as factors in the integral expression of Eq. 66 imply vanishing integrals ∫ − a a ∂ 2 ∂ z 1 2 1 | − z ̲ | sin 2 π λ i z i d z i = 0 , ∀ a ∈ R , i ∈ 1,2,3 , we arrive at A 1,1111 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 × ∫ R 3 ∂ 2 ∂ z 1 2 1 | − z ̲ | cos 2 π λ 1 z 1 × cos 2 π λ 2 z 2 cos 2 π λ 3 z 3 d V z ̲ .

Noting that ∂ ∂ z 1 1 | z ′ ̲ − z ̲ | = − ∂ ∂ z 1 ′ 1 | z ′ ̲ − z ̲ | and ∂ 2 ∂ z 1 2 1 | z ′ ̲ − z ̲ | = ∂ 2 ∂ z 1 ′ 2 1 | z ′ ̲ − z ̲ | , the integral in Eq. 70 can be expressed by means of an auxiliary function in z ̲ ′ , according to ∫ R 3 ∂ 2 ∂ z 1 2 1 | − z ̲ | cos 2 π λ 1 z 1 cos 2 π λ 2 z 2 cos 2 π λ 3 z 3 d V z ̲ = ∂ 2 ϕ ∂ z 1 ′ 2 z ̲ ′ = 0 , where the auxiliary function ϕ ( z ̲ ′ ) stands for ϕ z ̲ ′ = ∫ R 3 1 | z ′ ̲ − z ̲ | cos 2 π λ 1 z 1 cos 2 π λ 2 z 2 cos 2 π λ 3 z 3 d V z ̲ . Eq. 72 exhibits two remarkable properties. Firstly, its format ϕ z ̲ ′ = ∫ R 3 1 | z ′ ̲ − z ̲ | f z ̲ d V z ̲ is the solution of the Poisson’s equation ∇ 2 ϕ z ̲ ′ = − 4 π f z ̲ ′ . Secondly, Eq. 72 is symmetric in the sense of ϕ z 1 ′ , z 2 ′ , z 3 ′ = ϕ z 1 ′ , z 3 ′ , z 2 ′ = ϕ z 2 ′ , z 3 ′ , z 1 ′ = ϕ z 2 ′ , z 1 ′ , z 3 ′ = ϕ z 3 ′ , z 1 ′ , z 2 ′ = ϕ z 3 ′ , z 2 ′ , z 1 ′ → ∂ ϕ ∂ z 1 ′ = ∂ ϕ ∂ z 2 ′ = ∂ ϕ ∂ z 3 ′ → ∂ 2 ϕ ∂ z 1 ′ 2 = ∂ 2 ϕ ∂ z 2 ′ 2 = ∂ 2 ϕ ∂ z 3 ′ 2 .

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 ̲ ′ = 3 ∂ 2 ϕ ∂ z 1 ′ 2 = − 4 π cos 2 π λ 1 z 1 ′ cos 2 π λ 2 z 2 ′ cos 2 π λ 3 z 3 ′ , where we made use of Eq. 75. Solving the equation for ∂ 2 ϕ ∂ z 1 ′ 2 ( z ̲ ′ = 0 ) and inserting the corresponding result into Eq. 70 yield A 1,1111 0 x ̲ = 1 3 Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 . 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 A 2,1111 0 x ̲ = ∑ h = 1 3 ∑ ℓ = 1 3 ∑ i = 1 3 ∑ j = 1 3 ∫ R 3 ∫ R 3 G 11hℓ x ̲ − y ̲ × 3 Δ k sin 2 π λ 1 y 1 sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 I h ℓ i j vol × ∑ m = 1 3 ∑ n = 1 3 G ijmn y ̲ − y ′ ̲ × 3 Δ k sin 2 π λ 1 y 1 ′ sin 2 π λ 2 y 2 ′ × sin 2 π λ 3 y 3 ′ I m n 11 vol d V y ′ ̲ d V y ̲ .

Specification of Eq. 78 for Eq. 58, while considering the identity of Eq. 64, yields A 2,1111 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 2 ∫ R 3 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | × sin 2 π λ 1 y 1 sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 × ∇ 2 ∫ R 3 1 | y ̲ − y ̲ ′ | sin 2 π λ 1 y 1 ′ × sin 2 π λ 2 y 2 ′ sin 2 π λ 3 y 3 ′ d V y ̲ ′ d V y ̲ .

Considering the last integral in Eq. 79 as the Green’s function solving the Poisson’s equation ∇ 2 ϕ y ̲ = − 4 π sin 2 π λ 1 y 1 sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 yields A 2,1111 0 x ̲ = − 4 π − Δ k 1 + ν 0 12 π k 0 1 − ν 0 2 ∫ R 3 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | × sin 2 2 π λ 1 y 1 sin 2 2 π λ 2 y 2 sin 2 2 π λ 3 y 3 d V y ̲ . Recalling from Eq. 65 and Eq. 77 that ∫ R 3 ∂ 2 ∂ x 1 2 1 | x ̲ − y ̲ | f y ̲ d V y ̲ = − 4 π 3 f x ̲ , for any f ( y ̲ ) with the symmetry properties of Eq. 75, the integral in Eq. 81 can be solved, yielding A 2,1111 0 x ̲ = 1 3 Δ k 1 + ν 0 3 k 0 1 − ν 0 2 ⁡ sin 2 2 π λ 1 x 1 × sin 2 2 π λ 2 x 2 sin 2 2 π λ 3 x 3 . Repeating this derivation for the 1111-component of any other member of the series, i.e., for A n , 1111 0 with n > 2, one notices that A n , 1111 0 x ̲ = − 1 n 1 3 − Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 × sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 n .

Moreover, because of the symmetry of the considered micro-stiffness distribution, A n , 1111 0 ( x ̲ ) = A n , 2222 0 ( x ̲ ) = A n , 3333 0 ( 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 A n , i i ℓ ℓ 0 .

The components of the auxiliary strain downscaling tensor A 0 ( 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 = 1 1 − γ , f o r | γ | < 1 , is applied to α = − Δ k ( 1 + ν 0 ) 3 k 0 ( 1 − ν 0 ) sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 . This yields the normal strain concentration tensor components as A iiii 0 x ̲ = 2 3 + 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 , A iijj 0 x ̲ = − 1 3 + 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 .

5.3 Shear strain-related components of <inline-formula id="inf94"> <mml:math id="m181"> <mml:msup> <mml:mrow> <mml:mi mathvariant="double-struck">A</mml:mi> </mml:mrow> <mml:mrow> <mml:mn>0</mml:mn> </mml:mrow> </mml:msup> <mml:mrow> <mml:mo stretchy="false">(</mml:mo> <mml:mrow> <mml:munder accentunder="false"> <mml:mrow> <mml:mi>x</mml:mi> </mml:mrow> <mml:mo accent="true">̲</mml:mo> </mml:munder> </mml:mrow> <mml:mo stretchy="false">)</mml:mo> </mml:mrow> </mml:math> </inline-formula>

Next, we turn to the shear-related components of the concentration tensor, i.e., to A ijkℓ with i ≠ j or k ≠ ℓ. The concentration tensor is driven by the micro-stiffness fluctuation around the base stiffness C 0 , see the expressions of Eqs. 41–43, so that the chosen sinusoidal micro-stiffness distribution of Eq. 50, where the shear stiffness does not fluctuate around the basic contribution, implies that A n , i j k ℓ 0 x ̲ = 0 , ∀ n ∈ N , x ̲ ∈ V RVE , k ≠ ℓ . Combining Eq. 88 with Eq. 41 implies that A ijij 0 x ̲ = 1 , ∀ x ̲ ∈ V RVE , i ≠ j , and A ijkℓ 0 x ̲ = 0 , ∀ x ̲ ∈ V RVE , i ≠ j , k ≠ ℓ , i ≠ k , j ≠ ℓ .

Let us now turn to the remaining non-vanishing shear-related components of the concentration tensor, i.e., to A 1 , i j ℓ ℓ 0 ( x ̲ ) , with i ≠ j. For the microstiffness distribution of Eq. 50, the respective first term in the series of Eq. 41, defined by Eq. 42, reads as A 1 , i j ℓ ℓ 0 x ̲ = ∑ h = 1 3 ∑ k = 1 3 ∫ R 3 G 11hk x ̲ − y ̲ 3 Δ k sin 2 π λ 1 y 1 × sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 I h k ℓ ℓ vol d V y ̲ = Δ k ∑ h = 1 3 ∫ R 3 G 11hh x ̲ − y ̲ sin 2 π λ 1 y 1 × sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 d V y ̲ . Insertion of Eq. 58 into Eq. 91, while considering the identity resulting from derivation of Eq. 63 with respect to x i and x j , reading as ∂ 4 ∂ x i ∂ x j ∂ x 1 2 | x ̲ − y ̲ | + ∂ 4 ∂ x i ∂ x j ∂ x 2 2 | x ̲ − y ̲ | + ∂ 4 ∂ x i ∂ x j ∂ x 3 2 | x ̲ − y ̲ | = 2 ∂ 2 ∂ x i ∂ x j 1 | x ̲ − y ̲ | , yields A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 ∫ R 3 ∂ 2 ∂ x i ∂ x j 1 | x ̲ − y ̲ | × sin 2 π λ 1 y 1 sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 d V y ̲ . Introducing the variable z ̲ = y ̲ − x ̲ , using the relation of Eq. 67, and considering ∂ ∂ z j ∂ ∂ z i ≠ j 1 | − z ̲ | = 3 − z i − z j | − z ̲ | 5 , and the corresponding consequence of even and odd functions ∫ − a a ∂ ∂ z i 1 | − z ̲ | cos 2 π λ i z i d z i = 0 , ∀ a ∈ R , Eq. 93 can be transformed to A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i cos 2 π λ j x j × sin 2 π λ k x k ∫ R 3 3 z i z j | − z ̲ | 5 sin 2 π λ i z i × sin 2 π λ j z j cos 2 π λ k z k d V z ̲ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Introducing the microstructure-related variable change z ̲ ′ = 2 π λ 1 z 1 , 2 π λ 2 z 2 , 2 π λ 3 z 3 into Eq. 96 yields A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i × cos 2 π λ j x j sin 2 π λ k x k × ∫ R 3 3 z i ′ z j ′ | − z ′ ̲ | 5 sin z i ′ sin z j ′ cos z k ′ d V z ′ ̲ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Considering ∫ − ∞ ∞ z i ′ − z i ′ 2 + − z j ′ 2 + − z k ′ 2 5 / 2 sin z i ′ d z i ′ = 2 K 1 − z j ′ 2 + − z k ′ 2 3 − z j ′ 2 + − z k ′ 2 , with K 1 being the modified Bessel function of the second kind, Eq. 97 can be re-written as A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k × ∫ − ∞ ∞ ∫ − ∞ ∞ 2 z j ′ − z j ′ 2 + − z k ′ 2 × K 1 − z j ′ 2 + − z k ′ 2 sin z j ′ cos z k ′ d z j ′ d z k ′ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Applying a transformation towards polar coordinates, z j ′ = ρ ⁡ cos ⁡ θ and z k ′ = ρ ⁡ sin ⁡ θ , yields A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k × ∫ 0 2 π ∫ 0 ∞ 2 ρ ⁡ cos ⁡ θ K 1 ρ sin ρ ⁡ cos ⁡ θ cos ρ ⁡ sin ⁡ θ d ρ d θ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Thanks to the trigonometric identity sin ( a ) cos ( b ) = 1 2 [ sin ( a + b ) + sin ( a − b ) ] , Eq. 100 can be reformulated as A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k × ∫ 0 2 π ⁡ cos ⁡ θ ∫ 0 ∞ ρ K 1 ρ sin cos ⁡ θ + sin ⁡ θ ρ + sin cos ⁡ θ − sin ⁡ θ ρ d ρ d θ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Thus, integrating over ρ yields A 1 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k × ∫ 0 2 π ⁡ cos ⁡ θ cos ⁡ θ + sin ⁡ θ 2 + 2 ⁡ cos ⁡ θ ⁡ sin ⁡ θ + cos ⁡ θ − sin ⁡ θ 2 − 2 ⁡ cos ⁡ θ ⁡ sin ⁡ θ + sinh − 1 cos ⁡ θ + sin ⁡ θ 2 + 2 ⁡ cos ⁡ θ ⁡ sin ⁡ θ 3 / 2 + sinh − 1 cos ⁡ θ − sin ⁡ θ 2 − 2 ⁡ cos ⁡ θ ⁡ sin ⁡ θ 3 / 2 d θ , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. The integral remaining in Eq. 102 amounts to 4 π 3 . Thus, the first element of the series is A 1 , i j ℓ ℓ 0 x ̲ = 4 π 3 − Δ k 1 + ν 0 12 π k 0 1 − ν 0 cos 2 π λ i x i × cos 2 π λ j x j sin 2 π λ k x k , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j.

The following element of the series A 2 , i j l l 0 ( x ̲ ) reads, after the introduction of the corresponding Green’s functions from Eq. 58, as A 2 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 2 ∫ R 3 ∂ 2 ∂ x i ∂ x j 1 | x ̲ − y ̲ | × sin 2 π λ 1 y 1 sin 2 π λ 2 y 2 sin 2 π λ 3 y 3 × ∇ 2 ∫ R 3 1 | y ̲ − y ̲ ′ | sin 2 π λ 1 y 1 ′ × sin 2 π λ 2 y 2 ′ sin 2 π λ 3 y 3 ′ d V y ̲ ′ d V y ̲ .

After application of Poisson’s Eq. 74 and Eq. 104 can be expressed as follows, when using sin 2 ( a ) = 1 2 [ 1 − cos ( 2 a ) ] A 2 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 2 1 8 ∫ R 3 ∂ 2 ∂ x i ∂ x j 1 | x ̲ − y ̲ | × 1 − cos 2 2 π λ 1 y 1 1 − cos 2 2 π λ 2 y 2 × 1 − cos 2 2 π λ 3 y 3 d V y ̲ .

Proceeding with an analogous process to the one carried out to obtain A 1 , i j l l 0 ( x ̲ ) , i.e., variable change to z ̲ = y ̲ − x ̲ , partial derivation with respect to z i and z j , change of the variable z ̲ ′ = 2 π λ 1 z 1 , 2 π λ 2 z 2 , 2 π λ 3 z 3 , integration over z i ′ , change towards polar coordinate system ( z j ′ = ρ ⁡ cos ⁡ θ and z k ′ = ρ ⁡ sin ⁡ θ ), integration over r and, lastly, integration over θ yield A 2 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 2 ⁡ sin 2 2 π λ i x i sin 2 2 π λ j x j × π 4 − π 6 cos 2 2 π λ k x k , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Similarly, the third element was computed as A 3 , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 3 9 π 16 cos 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k − 27 π 176 cos 3 2 π λ i x i cos 2 π λ j x j sin 2 π λ k x k − 27 π 176 cos 2 π λ i x i cos 3 2 π λ j x j sin 2 π λ k x k − 9 π 176 cos 2 π λ i x i cos 2 π λ j x j sin 3 2 π λ k x k + 27 π 304 cos 3 2 π λ i x i cos 3 2 π λ j x j sin 2 π λ k x k + 9 π 304 cos 3 2 π λ i x i cos 2 π λ j x j sin 3 2 π λ k x k + 9 π 304 cos 2 π λ i x i cos 3 2 π λ j x j sin 3 2 π λ k x k − π 48 cos 3 2 π λ i x i cos 3 2 π λ j x j sin 3 2 π λ k x k , with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j.

For the following terms A n , i j l l 0 ( x ̲ ) , corresponding substitution of Green’s function tensor components of Eq. 58 and reiterated application of Poisson’s Eq. 74 yield A n , i j ℓ ℓ 0 x ̲ = − Δ k 1 + ν 0 12 π k 0 1 − ν 0 n ∫ R 3 ∂ 2 ∂ x i ∂ x j 1 | x ̲ − y ̲ | × sin n 2 π λ 1 y 1 sin n 2 π λ 2 y 2 sin n 2 π λ 3 y 3 d V y ̲ . These terms can be computed in the same manner as the previous ones, converging rapidly due to the factor − Δ k ( 1 + ν 0 ) 12 π k 0 ( 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 ̲ ∈ V RVE = n λ 3 λ 1 λ 2 λ 3 ⇔ n λ − λ 1 / 2 ≤ x 1 ≤ n λ λ 1 / 2 , n λ − λ 2 / 2 ≤ x 2 ≤ n λ λ 2 / 2 , n λ − λ 3 / 2 ≤ x 3 ≤ n λ λ 3 / 2 , 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, A 0 ( 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 , M − 1 , will be calculated first, in order to obtain the RVE-to-auxiliary strain conversion tensor as M = M − 1 − 1 .

Like in the previous sections, the components of tensor M − 1 will be obtained individually. The first components studied are M − 1 iiii . Insertion of Eq. 86 into the inverse of Eq. 12 yields ( M − 1 ) iiii = 2 3 1 V RVE ∫ V RVE d V x ̲ + 1 V RVE ∫ V RVE 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 d V x ̲ . Clearly, the term in square brackets is equal to 1, while we must focus on the other integral expression I = 1 V RVE ∫ n λ − λ 1 / 2 n λ λ 1 / 2 ∫ n λ − λ 2 / 2 n λ λ 2 / 2 ∫ n λ − λ 3 / 2 n λ λ 3 / 2 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 × d x 1 d x 2 d x 3 = 1 λ 1 λ 2 λ 3 ∫ − λ 1 / 2 λ 1 / 2 ∫ − λ 2 / 2 λ 2 / 2 ∫ − λ 3 / 2 λ 3 / 2 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 × d x 1 d x 2 d x 3 . Proceeding with a change of variable x ̲ ′ = 2 π λ 1 x 1 , 2 π λ 2 x 2 , 2 π λ 3 x 3 yields I = 1 2 π 3 ∫ − π π ∫ − π π ∫ − π π 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin x 1 ′ sin x 2 ′ sin x 3 ′ × d x 1 ′ d x 2 ′ d x 3 ′ . Solving the integral for x 1 ′ yields I = 1 2 π 3 ∫ − π π ∫ − π π 2 ⁡ tan − 1 ψ x 2 ′ , x 3 ′ + tan π / 2 1 − ψ 2 x 2 ′ , x 3 ′ 1 − ψ 2 x 2 ′ , x 3 ′ d x 2 ′ d x 3 ′ − ∫ − π π ∫ − π π 2 ⁡ tan − 1 ψ x 2 ′ , x 3 ′ + tan − π / 2 1 − ψ 2 x 2 ′ , x 3 ′ 1 − ψ 2 x 2 ′ , x 3 ′ d x 2 ′ d x 3 ′ = 1 2 π 2 ∫ − π π ∫ − π π 1 1 − ψ 2 x 2 ′ , x 3 ′ d x 2 ′ d x 3 ′ , where ψ ( x 2 ′ , x 3 ′ ) = Δ k ( 1 + ν 0 ) 3 k 0 ( 1 − ν 0 ) sin ( x 2 ′ ) sin ( x 3 ′ ) . Integrating I with respect to x 2 ′ yields I = 2 2 π 2 ∫ − π π K χ 2 ⁡ sin 2 x 3 ′ + 1 1 − χ 2 ⁡ sin 2 x 3 ′ K χ 2 ⁡ sin 2 x 3 ′ 1 − χ 2 ⁡ sin 2 x 3 ′ d x 3 ′ , where K(a) is the complete elliptic integral of the first kind with parameter a, and χ = Δ k 1 + ν 0 3 k 0 1 − ν 0 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,

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

M − 1 iiii = 2 3 + I χ .

The next components to be considered are M iijj − 1 , with i ≠ j. They read, from Eq. 87 and Eq. 115, as M − 1 iijj = − 1 3 + I χ , ∀ i ≠ j . The shear components read as M − 1 ijij = 1 , ∀ i ≠ j , and M − 1 ijkℓ = 0 , ∀ i ≠ j , i ≠ k . Thus, the auxiliary-to-RVE tensor reads as M = 1 3 I 1 / 3 + 2 I 1 / 3 − I 1 / 3 − I 0 0 0 1 / 3 − I 1 / 3 + 2 I 1 / 3 − I 0 0 0 1 / 3 − I 1 / 3 − I 1 / 3 + 2 I 0 0 0 0 0 0 3 I 0 0 0 0 0 0 3 I 0 0 0 0 0 0 3 I .

Lastly, the real strain concentration tensor field is computed according to Eq. 14 as A x ̲ = A iiii A iijj A iijj 0 0 0 A iijj A iiii A iijj 0 0 0 A iijj A iijj A iiii 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 , where A iiii x ̲ = 2 3 + 1 3 I χ 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 , A iijj x ̲ = − 1 3 + 1 3 I χ 1 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 , 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 C 0 reads as Δ C hom = C hom − C 0 = 1 V RVE ∫ V RVE c x ̲ − C 0 : A x ̲ d V x ̲ . One more time, the components of this tensor will be obtained individually. The only non-vanishing components ΔC hom,ijkℓ are those with i = j and k = ℓ, due to Eq. 122 and c ijkℓ x ̲ − C ijkℓ 0 = Δ k sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 sin 2 π λ 3 x 3 δ i j δ k ℓ . Therefore, the components ΔC hom,iikk read as Δ C hom , i i k k = 1 V RVE ∫ V RVE Δ k sin 2 π λ 1 x 1 sin 2 π λ 2 x 2 × sin 2 π λ 3 x 3 ∑ p = 1 3 A ppkk x ̲ d V x ̲ . Thus, inserting Eq. 123 and Eq. 124 into Eq. 127 and applying x ̲ ′ = 2 π λ x ̲ yields Δ C hom , i i k k = Δ k I 2 π 3 ∫ V RVE ′ sin x 1 ′ sin x 2 ′ sin x 3 ′ 1 + Δ k 1 + ν 0 3 k 0 1 − ν 0 sin x 1 ′ sin x 2 ′ sin x 3 ′ d V x ̲ ′ . Integrating Eq. 128 with respect to x 1 ′ from −π to π yields Δ C hom , i i k ℓ = 3 k 0 1 − ν 0 I 1 + ν 0 2 π 2 ∫ − π π ∫ − π π d x 2 ′ d x 3 ′ − ∫ − π π ∫ − π π 1 1 − ψ 2 x 2 ′ , x 3 ′ d x 2 ′ d x 3 ′ = − 3 k 0 1 − ν 0 1 + ν 0 I χ − 1 I χ ,

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. 41–43, 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 ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | f z ̲ d V z ̲ = ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | f z ̲ − f 0 d V z ̲ + ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | d V z ̲ f 0 , whereby f z ̲ = cos 2 π λ z 1 cos 2 π λ z 2 cos 2 π λ z 3 , 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 ϵ → 0 ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | f z ̲ − f 0 d V z ̲ = 0 , so that lim ϵ → 0 ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | f z ̲ d V z ̲ = lim ϵ → 0 ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | d V z ̲ f 0 = lim ϵ → 0 ∫ S ϵ ∂ ∂ z 1 1 | z ̲ | n 1 d S z ̲ f 0 , 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 ̲ | = z 1 2 + z 2 2 + z 3 2 1 / 2 , the integrand in the surface integral of Eq. 133 can be transformed to ∂ ∂ z 1 z 1 2 + z 2 2 + z 3 2 − 1 / 2 = − z 1 z 1 2 + z 2 2 + z 3 2 − 3 / 2 = − z 1 ϵ 3 . Then, the surface integral in Eq. 133 is preferably evaluated in spherical coordinates, where. z 1 = ϵ cos ⁡ ϕ sin ⁡ θ , n 1 = cos ⁡ ϕ sin ⁡ θ , d S = ϵ 2 sin ⁡ θ d ϕ d θ , so that use of Eqs. 135–138 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 ϵ → 0 ∫ S ϵ ∂ ∂ z 1 1 | z ̲ | n 1 d S z ̲ f 0 = − ∫ θ = 0 π ⁡ sin 3 ⁡ θ d θ ∫ ϕ = 0 2 π ⁡ cos 2 ⁡ ϕ d ϕ f 0 = − 4 π 3 f 0 , and hence, since f (0) = 1, ∫ V ϵ ∂ 2 ∂ z 1 2 1 | z ̲ | f z ̲ d V z ̲ = − 4 π 3 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. 71–77. 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. 94–103.

 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).

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

Mechanical properties
average bulk modulus k 0 = 100 [MPa]
bulk modulus fluctuation Δk = 75 [MPa]
Poisson ratio ν = 0.3
Associated quantities
Eq. 116 χ = 0.46429
Eq. 115, Figure 2 I(χ) = 1.02975

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 ̲ E 0 − G ̂ k ̲ : τ ̂ k ̲ , with the wave vector k ̲ , and the hat-symbol indicating the Fourier transform, according to ε ̂ k ̲ = 1 2 π 3 ∫ R 3 ε x ̲ exp − i k ̲ ⋅ x ̲ d V x ̲ . 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 A iiii , 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 ̲ i ⊗ e ̲ 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 E ii = 1. Analogously, component A iijj , with i, j = 1, 2, 3 and i ≠ j, is the microstrain field component ɛ ii , with i = 1, 2, 3, arising from a macroscopic strain tensor of the format E = e ̲ j ⊗ e ̲ 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.

Strain concentration field A iiii ( x ̲ ) , i = 1, 2, 3, arising from the sinusoidal microelasticity distribution of Eq. 50, at the plane x i = λ 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.

Strain concentration field A iijj ( x ̲ ) , i, j, = 1, 2, 3, at the plane x i = λ 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.

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.

Stiffness component Analytical model [MPa] FFT [MPa]
323 voxels 523 voxels 1003 voxels
C 1111 0 161.54 –
C hom,1111 156.87 157.01 156.92 156.89
C 1122 0 69.23 –
C hom,1122 64.56 64.70 64.62 64.58
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 grad x grad x G x ̲ − y ̲ = − grad x S grad y G x ̲ − y ̲ . 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 E 0: 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. 41–43, 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 E 0. 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 ̲ , w i t h ‖ σ ‖ ‖ ∂ σ ∂ x ̲ ‖ = d . 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 (x 1, x 2, x 3) yields c x ̲ = C 0 + 3 ∑ k ∞ ∑ l ∞ ∑ m ∞ c klm ⁡ exp 2 π i k x 1 λ 1 + l x 2 λ 2 + m x 3 λ 3 , with i now standing for the imaginary unit, and with the Fourier coefficients c klm being obtained from c klm = ∫ − a a ∫ − a a ∫ − a a ⁡ exp − 2 π i k x 1 λ 1 + l x 2 λ 2 + m x 3 λ 3 × Δ k x 1 , x 2 , x 3 d x 1 d x 2 d x 3 . 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 cos a + b = cos a cos b − sin a sin b , 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. 122–124. 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, C hom , is smaller than or equal to the average stiffness, C 0 . 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 : ⟨ A T : c : A ⟩ − C hom : E ≥ 0 , ∀ E .

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. 10.1016/0167-6636(87)90005-6 Benveniste Y. Dvorak G. J. Chen T. (1991). On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media. J. Mech. Phys. Solids 39, 927–946. 10.1016/0022-5096(91)90012-d Bernard O. Ulm F.-J. Lemarchand E. (2003). A multiscale micromechanics-hydration model for the early-age elastic properties of cement-based materials. Cem. Concr. Res. 33, 1293–1309. 10.1016/s0008-8846(03)00039-5 Bertrand E. Hellmich C. (2009). Multiscale elasticity of tissue engineering scaffolds with tissue-engineered bone: A continuum micromechanics approach. J. Eng. Mech. 135, 395–412. 10.1061/(asce)0733-9399(2009)135:5(395) Brisard S. Dormieux L. (2010). FFT-based methods for the mechanics of composites: A general variational framework. Comput. Mater. Sci. 49, 663–671. 10.1016/j.commatsci.2010.06.009 Buchner T. Königsberger M. Jäger A. Füssl J. (2022). A validated multiscale model linking microstructural features of fired clay brick to its macroscopic multiaxial strength. Mech. Mater. 170, 104334. 10.1016/j.mechmat.2022.104334 Buryachenko V. (2007). Micromechanics of heterogeneous materials. Springer Science and Business Media. Cai X. Brenner R. Peralta L. Olivier C. Gouttenoire P.-J. Chappard C. (2019). Homogenization of cortical bone reveals that the organization and shape of pores marginally affect elasticity. J. R. Soc. Interface 16, 20180911. 10.1098/rsif.2018.0911 Chung M. K. Dalton K. M. Shen L. Evans A. C. 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. 10.1109/tmi.2007.892519 Cochran W. Cooley J. Favin D. Helms H. Kaenel R. Lang W. (1967). What is the fast Fourier transform? Proc. IEEE 55, 1664–1674. 10.1109/PROC.1967.5957 Dvorak G. J. (2012). Micromechanics of composite materials. Springer Science and Business Media. 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. Feng T. Jia M. Xu W. Wang F. Li P. Wang X. (2022). Three-dimensional mesoscopic investigation of the compression mechanical properties of ultra-high performance concrete containing coarse aggregates. Cem. Concr. Compos. 133, 104678. 10.1016/j.cemconcomp.2022.104678 Fredholm I. (1900). Sur les équations de l’équilibre d’un corps solide élastique. Acta Math. 23, 1–42. 10.1007/bf02418668 Fritsch A. Dormieux L. Hellmich C. (2006). Porous polycrystals built up by uniformly and axisymmetrically oriented needles: Homogenization of elastic properties. Comptes Rendus Mécanique 334, 151–157. 10.1016/j.crme.2006.01.008 Fritsch A. Dormieux L. Hellmich C. 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. 10.1002/jbm.a.31727 Fritsch A. Hellmich C. 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. 10.1016/j.jtbi.2009.05.021 Fritsch A. Hellmich C. Young P. (2013). Micromechanics-derived scaling relations for poroelasticity and strength of brittle porous polycrystals. J. Appl. Mech. 80, 020905. 10.1115/1.4007922 Fritsch A. 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. 10.1016/j.jtbi.2006.09.013 Grimal Q. Raum K. Gerisch A. 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. 10.1007/s10237-010-0284-9 Guo W. Han F. Jiang J. Xu W. (2022). A micromechanical framework for thermo-elastic properties of multiphase cementitious composites with different saturation. Int. J. Mech. Sci. 224, 107313. 10.1016/j.ijmecsci.2022.107313 Hashin Z. (1983). Analysis of composite materials—a survey. J. Appl. Mech. 50, 481–505. 10.1115/1.3167081 Hashin Z. (1963). Theory of mechanical behavior of heterogeneous media. Philidelphia, Pa: Towne School of Civil and Mechanical Engineering, University of Pennsylvania. Hashin Z. (1965). Viscoelastic behavior of heterogeneous media. J. Appl. Mech. 32, 630–636. 10.1115/1.3627270 Hellmich C. Barthélémy J.-F. Dormieux L. (2004). Mineral–collagen interactions in elasticity of bone ultrastructure – A continuum micromechanics approach. Eur. J. Mech. - A/Solids 23, 783–810. 10.1016/j.euromechsol.2004.05.004 Hellmich C. Mang H. (2005). Shotcrete elasticity revisited in the framework of continuum micromechanics: From submicron to meter level. J. Mater. Civ. Eng. 17, 246–256. 10.1061/(asce)0899-1561(2005)17:3(246) Hervé E. Zaoui A. (1993). inclusion-based micromechanical modelling. Int. J. Eng. Sci. 31, 1–10. 10.1016/0020-7225(93)90059-4 Hill R. (1963). Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids 11, 357–372. 10.1016/0022-5096(63)90036-x Hofstetter K. Hellmich C. 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. 10.1016/j.euromechsol.2005.05.006 Horwitz A. (2001). A version of Simpson’s rule for multiple integrals. J. Comput. Appl. Math. 134, 1–11. 10.1016/s0377-0427(00)00444-1 Jiménez Segura N. Pichler B. L. Hellmich C. (2023). Concentration tensors preserving elastic symmetry of multiphase composites. Mech. Mater. 178, 104555. 10.1016/j.mechmat.2023.104555 Jiménez Segura N. Pichler B. L. 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. 10.1002/zamm.202200091 Kneer G. (1965). Über die Berechnung der Elastizitätsmoduln vielkristalliner Aggregate mit Textur. Phys. Status Solidi (b) 9, 825–838. 10.1002/pssb.19650090319 Königsberger M. Hlobil M. Delsaute B. Staquet S. Hellmich C. Pichler B. (2018). Hydrate failure in itz governs concrete strength: A micro-to-macro validated engineering mechanics model. Cem. Concr. Res. 103, 77–94. 10.1016/j.cemconres.2017.10.002 Korringa J. (1973). Theory of elastic constants of heterogeneous media. J. Math. Phys. 14, 509–513. 10.1063/1.1666346 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. Kröner E. (1977). Bounds for effective elastic moduli of disordered materials. J. Mech. Phys. Solids 25, 137–155. 10.1016/0022-5096(77)90009-6 Levin V. M. (1967). Thermal expansion coefficient of heterogeneous materials. Mekhanika Tverd. Tela 2, 83–94. Lipinski P. Barhdadi E. H. Cherkaoui M. (2006). Micromechanical modelling of an arbitrary ellipsoidal multi-coated inclusion. Philos. Mag. 86, 1305–1326. 10.1080/14786430500343868 Lippmann B. A. Schwinger J. (1950). Variational principles for scattering processes. I. Phys. Rev. 79, 469–480. 10.1103/physrev.79.469 Lucarini S. Upadhyay M. Segurado J. (2021). Fft based approaches in micromechanics: Fundamentals, methods and applications. Model. Simul. Mater. Sci. Eng. 30, 023002. 10.1088/1361-651x/ac34e1 Moës N. Cloirec M. Cartraud P. Remacle J. F. (2003). A computational approach to handle complex microstructure geometries. Comput. Methods Appl. Mech. Eng. 192, 3163–3177. 10.1016/s0045-7825(03)00346-3 Molinari A. El Mouden M. (1996). The problem of elastic inclusions at finite concentration. Int. J. Solids Struct. 33, 3131–3150. 10.1016/0020-7683(95)00275-8 Mori T. Tanaka K. (1973). Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 21, 571–574. 10.1016/0001-6160(73)90064-3 Moulinec H. 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. Moulinec H. 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. 10.1016/s0045-7825(97)00218-1 Pahr D. H. Zysset P. K. (2008). Influence of boundary conditions on computed apparent elastic properties of cancellous bone. Biomechanics Model. Mechanobiol. 7, 463–476. 10.1007/s10237-007-0109-7 Pichler B. Hellmich C. Eberhardsteiner J. Wasserbauer J. Termkhajornkit P. Barbarulo R. (2013). Effect of gel–space ratio and microstructure on strength of hydrating cementitious materials: An engineering micromechanics approach. Cem. Concr. Res. 45, 55–68. 10.1016/j.cemconres.2012.10.019 Pichler B. Hellmich C. (2011). Upscaling quasi-brittle strength of cement paste and mortar: A multi-scale engineering mechanics model. Cem. Concr. Res. 41, 467–476. 10.1016/j.cemconres.2011.01.010 Rosen B. W. Hashin Z. (1970). Effective thermal expansion coefficients and specific heats of composite materials. Int. J. Eng. Sci. 8, 157–173. 10.1016/0020-7225(70)90066-2 Sanahuja J. Dormieux L. Meille S. Hellmich C. 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. 10.1061/(asce)em.1943-7889.0000072 Scheiner S. Sinibaldi R. Pichler B. Komlev V. Renghini C. Vitale-Brovarone C. (2009). Micromechanics of bone tissue-engineering scaffolds, based on resolution error-cleared computer tomography. Biomaterials 30, 2411–2419. 10.1016/j.biomaterials.2008.12.048 Sevostianov I. Kachanov M. (2014). On some controversial issues in effective field approaches to the problem of the overall elastic properties. Mech. Mater. 69, 93–105. 10.1016/j.mechmat.2013.09.010 Ting T. C. T. 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. 10.1093/qjmam/50.3.407 Tonon F. Pan E. Amadei B. (2001). Green’s functions and boundary element method formulation for 3D anisotropic media. Comput. Struct. 79, 469–482. 10.1016/s0045-7949(00)00163-2 Torquato S. (1997). Effective stiffness tensor of composite media—I. Exact series expansions. J. Mech. Phys. Solids 45, 1421–1448. 10.1016/s0022-5096(97)00019-7 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. 10.1002/andp.18892741206 Wang H. Hellmich C. Yuan Y. Mang H. 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. 10.1016/j.cemconres.2018.05.008 Whittaker E. 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. Willis J. R. (1977). Bounds and self-consistent estimates for the overall properties of anisotropic composites. J. Mech. Phys. Solids 25, 185–202. 10.1016/0022-5096(77)90022-9 Wolfram U. Peña Fernández M. McPhee S. Smith E. Beck R. J. Shephard J. D. (2022). Multiscale mechanical consequences of ocean acidification for cold-water corals. Sci. Rep. 12, 8052. 10.1038/s41598-022-11266-w Xie L. Zhang C. Sladek J. 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. 10.1098/rspa.2015.0272 Xu W. Wu F. Jiao Y. 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. 10.1016/j.matdes.2017.04.075 Xu W. Wu Y. 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. 10.1016/j.cma.2019.03.021 Xu W. Wu Y. 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. 10.1016/j.compstruct.2018.07.009 Zaoui A. (2002). Continuum micromechanics: Survey. J. Eng. Mech. 128, 808–816. 10.1061/(asce)0733-9399(2002)128:8(808) Zeller R. Dederichs P. (1973). Elastic constants of polycrystals. Phys. Status Solidi (b) 55, 831–842. 10.1002/pssb.2220550241 Zienkiewicz O. C. Taylor R. L. Zhu J. Z. (2005). The finite element method: Its basis and fundamentals. Elsevier.