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

Abstract

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

1 Introduction–Motivation and scope

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 ; 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 , somewhat neglecting the concentration tensor field . 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 , from a given microstiffness distribution 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⁰. The solution of Eshelby (1957) then relates the auxiliary strains to the homogeneous strains ɛ^I in the inhomogeneity, according towith as the concentration tensor associated with inhomogeneity I embedded in an infinite matrix subjected to E⁰. 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 formwith as a concentration tensor field associated with an arbitrary inhomogeneous microstructure represented by a stiffness distribution , spreading throughout an infinite domain, see Figure 1C. The stiffness distribution is considered as fluctuation around a homogeneous stiffness . The strains then consist of two portions: (i) auxiliary strains E⁰ prevailing in the homogeneous infinite domain of stiffness , see Figure 1A, and (ii) fluctuations around E⁰ which arise from the fluctuations in the stiffness field, , 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 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 thatwith as the RVE-to-auxiliary strain conversion tensor. Multiplication from the left, of Eq. 12 with , and insertion of the corresponding result into Eq. 11, yieldsand comparison of Eq. 13 with Eq. 7 yields the strain concentration tensor asand when considering, in addition, Eq. 9, the homogenized stiffness tensor is eventually retrieved assee Figures 1E, F. Eq. 15 can be reformulated by setting , which yieldsWe are left with the determination of . 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

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 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)Equating Eq. 17 with Eq. 4 yieldssee 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

with homogeneous boundary conditionsThe homogeneous solution reads asThe particular solution can be given in the formwhere 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)with boundary conditionsIn Eq. 27, 1 denotes the second-order identity tensor and δ denotes the Dirac function, with the following properties.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, yieldingInsertion of Eq. 32 into Eq. 26, adding the respective result to Eq. 25, and using the obtained expression in Eq. 3 yield the microscopic strain field throughout the RVE asIn Eq. 33, the following gradients of the Green’s functions were explicitly introducedConsidering 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 yieldsEq. 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 in the integral of Eq. 36. Accordingly, the term in the integral of Eq. 36 is numbered in a way reflecting these insertion processes, namely, by , i = 1, 2, …, ∞.

As a starting point, Eq. 36 is specified for in the integral of Eq. 36, yieldingIn order to come up with an expression for to be inserted into the integral in Eq. 37, we specify Eq. 37 for and for , yieldingInsertion of Eq. 38 into Eq. 37 yields

where 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, . For this purpose, we specify Eq. 37 for and . Repeating this process over and over again gives access to a relation involving an auxiliary macrostrain concentration tensor field and a residual term comprising the implicit strain field , reading asIn more detail, the concentration tensor field reads aswhereand, for n > 1,

The residual term after N iterations reads as

Note that for N → ∞, , provided the following requirement is met: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 , 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 by , which is integrated over the volume of the RVE, V_RVE, (while every other is integrated over the entire space, ), allows for the following expressionThus, a comprehenive integral format for the strain concentration tensor field is obtained from inserting Eqs. 41–43 and Eq. 46 into Eq. 14, yielding

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 by , yielding

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₀ and Poisson’s ratio ν₀; it reads as (Dvorak, 2012)where, with respect to the actual formula (4.5.12) given on page 107 of (Dvorak, 2012), we consider , with μ₀ being the shear modulus. Moreover, we consider a sinusoidal stiffness distribution across this infinitely extended body, according towhere x₁, x₂, and x₃ are the components of location vector with respect to an orthonormal base frame , , , such that , and stands for the volumetric part of the symmetric fourth-order unity tensor . The latter has the components I_ijrs = 1/2 (δ_irδ_js + δ_isδ_jr), while the former reads as , with components , 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

The auxiliary strain downscaling tensor is the result of the sum of an infinite series, see Eq. 41. The first term of this series, , is defined through Eq. 42, so that consideration of Eq. 50 yields the component 1111 of asIn order to retrieve the components G₁₁₁₁, G₁₁₂₂, and G₁₁₃₃, we start with the general expression for the components of the Green’s function of Eq. 49, reading asThe fourth-order Green’s function gradient according to Eq. 35 exhibits the following components

In order to evaluate Eq. 53, it is useful to recall the following properties of the spatial derivatives of the norm . The first-order derivative of the aforementioned norm reads asrevealing the interesting propertyThe first derivative of the inverse of the norm  reads asrevealing the interesting propertyEq. 55 and Eq. 57 allow for re-writing Eq. 53 as

so that the components occurring in Eq. 51 read asThe 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,Namely, Eq. 62 entails the following identitytwofold derivation of which yields an identity comprising the derivatives of the norm  occurring in Eqs. 59–61, reading asInsertion of Eqs. 59–61 into Eq. 51, while considering Eq. 64, yieldsA change of variable to results in

Transforming Eq. 66 by means ofand considering the second derivative of the inverse of the norm of asas well as that the even and odd functions appearing as factors in the integral expression of Eq. 66 imply vanishing integralswe arrive at

Noting that and , the integral in Eq. 70 can be expressed by means of an auxiliary function in , according towhere the auxiliary function stands forEq. 72 exhibits two remarkable properties. Firstly, its formatis the solution of the Poisson’s equationSecondly, Eq. 72 is symmetric in the sense of

Use of the latter relations in Eq. 74, while accounting for the structure of f according to Eq. 72 and Eq. 73, yieldswhere we made use of Eq. 75. Solving the equation for and inserting the corresponding result into Eq. 70 yieldThe 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

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

Considering the last integral in Eq. 79 as the Green’s function solving the Poisson’s equationyieldsRecalling from Eq. 65 and Eq. 77 thatfor any with the symmetry properties of Eq. 75, the integral in Eq. 81 can be solved, yieldingRepeating this derivation for the 1111-component of any other member of the series, i.e., for with n > 2, one notices that

Moreover, because of the symmetry of the considered micro-stiffness distribution, 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 .

The components of the auxiliary strain downscaling tensor are calculated as an infinite sum, see Eq. 41. Therefore, the explicit expression for the sum of an infinite geometric series, reading asis applied to . This yields the normal strain concentration tensor components as

5.3 Shear strain-related components of

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 , 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 thatCombining Eq. 88 with Eq. 41 implies thatand

Let us now turn to the remaining non-vanishing shear-related components of the concentration tensor, i.e., to , 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 asInsertion 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 asyieldsIntroducing the variable , using the relation of Eq. 67, and consideringand the corresponding consequence of even and odd functionsEq. 93 can be transformed towith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Introducing the microstructure-related variable change into Eq. 96 yieldswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Consideringwith K₁ being the modified Bessel function of the second kind, Eq. 97 can be re-written aswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Applying a transformation towards polar coordinates, and , yields with i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Thanks to the trigonometric identity , Eq. 100 can be reformulated aswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Thus, integrating over ρ yieldswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. The integral remaining in Eq. 102 amounts to . Thus, the first element of the series iswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j.

The following element of the series reads, after the introduction of the corresponding Green’s functions from Eq. 58, as

After application of Poisson’s Eq. 74 and Eq. 104 can be expressed as follows, when using

Proceeding with an analogous process to the one carried out to obtain , i.e., variable change to , partial derivation with respect to z_i and z_j, change of the variable , integration over , change towards polar coordinate system ( and ), integration over r and, lastly, integration over θ yieldwith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j. Similarly, the third element was computed aswith i,j,l = 1,2,3, whereby i ≠ j and k ≠ i,j.

For the following terms , corresponding substitution of Green’s function tensor components of Eq. 58 and reiterated application of Poisson’s Eq. 74 yieldThese terms can be computed in the same manner as the previous ones, converging rapidly due to the factor .

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., 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, , is related to the auxiliary concentration tensor calculated in the previous sections, , according to Eq. 14. Thus, this section is devoted to the derivation of the tensorial link , see Eq. 12. For the sake of simplicity, the inverse of , , will be calculated first, in order to obtain the RVE-to-auxiliary strain conversion tensor as

Like in the previous sections, the components of tensor will be obtained individually. The first components studied are . Insertion of Eq. 86 into the inverse of Eq. 12 yieldsClearly, the term in square brackets is equal to 1, while we must focus on the other integral expressionProceeding with a change of variable yieldsSolving the integral for yieldswhere . Integrating I with respect to yieldswhere K(a) is the complete elliptic integral of the first kind with parameter a, andis 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

The next components to be considered are , with i ≠ j. They read, from Eq. 87 and Eq. 115, asThe shear components read asandThus, the auxiliary-to-RVE tensor reads as

Lastly, the real strain concentration tensor field is computed according to Eq. 14 aswheresee 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 reads asOne 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 andTherefore, the components ΔC_hom,iikk read asThus, inserting Eq. 123 and Eq. 124 into Eq. 127 and applying yieldsIntegrating Eq. 128 with respect to from −π to π yields

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 , 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 , i.e., at . 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 . 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,wherebyis fully in line with the developments of Eq. 72 and Eq. 73. As stated before, we are interested in the limit case of ϵ → 0 whereso thatwhereby we made use of the divergence theorem, with 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 thatthe integrand in the surface integral of Eq. 133 can be transformed toThen, the surface integral in Eq. 133 is preferably evaluated in spherical coordinates, where.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 haveand hence, since f (0) = 1,The result of Eq. 140 is fully equivalent with evaluating Eq. 76 for 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 , 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).

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 , reading aswith the wave vector , and the hat-symbol indicating the Fourier transform, according toNamely, 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 , with 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 .

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 100³ voxels last by a factor of over 6,500 longer, namely, for 26.7 s.

FIGURE 3

FIGURE 4

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 – this differs from our definition (35). Mathematically speaking, the aforementioned sign change is due toHowever, 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 ɛ⁰ plays exactly the role of our auxiliary strain E⁰: 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 ɛ⁰, the role of which is comparable to our auxiliary strain E⁰. 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 , rather than , 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 ɛ⁰.

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 asDue 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₁, x₂, x₃) yieldswith i now standing for the imaginary unit, and with the Fourier coefficients c_klm being obtained fromIn 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 andthe 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, , is smaller than or equal to the average stiffness, . 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)

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.

References

• 1

  BenvenisteY. (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

  • CrossRef
  • Google Scholar

• 2

  BenvenisteY.DvorakG. J.ChenT. (1991). On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media. J. Mech. Phys. Solids39, 927–946. 10.1016/0022-5096(91)90012-d

  • CrossRef
  • Google Scholar

• 3

  BernardO.UlmF.-J.LemarchandE. (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

  • CrossRef
  • Google Scholar

• 4

  BertrandE.HellmichC. (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)

  • CrossRef
  • Google Scholar

• 5

  BrisardS.DormieuxL. (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

  • CrossRef
  • Google Scholar

• 6

  BuchnerT.KönigsbergerM.JägerA.FüsslJ. (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

  • CrossRef
  • Google Scholar

• 7

  BuryachenkoV. (2007). Micromechanics of heterogeneous materials. Springer Science and Business Media.

  • Google Scholar

• 8

  CaiX.BrennerR.PeraltaL.OlivierC.GouttenoireP.-J.ChappardC.et al (2019). Homogenization of cortical bone reveals that the organization and shape of pores marginally affect elasticity. J. R. Soc. Interface16, 20180911. 10.1098/rsif.2018.0911

  • CrossRef
  • Google Scholar

• 9

  ChungM. K.DaltonK. M.ShenL.EvansA. C.DavidsonR. J. (2007). Weighted Fourier series representation and its application to quantifying the amount of gray matter. IEEE Trans. Med. Imaging26, 566–581. 10.1109/tmi.2007.892519

  • CrossRef
  • Google Scholar

• 10

  CochranW.CooleyJ.FavinD.HelmsH.KaenelR.LangW.et al (1967). What is the fast Fourier transform?Proc. IEEE55, 1664–1674. 10.1109/PROC.1967.5957

  • CrossRef
  • Google Scholar

• 11

  DvorakG. J. (2012). Micromechanics of composite materials. Springer Science and Business Media.

  • Google Scholar

• 12

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

• 13

  FengT.JiaM.XuW.WangF.LiP.WangX.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. 10.1016/j.cemconcomp.2022.104678

  • CrossRef
  • Google Scholar

• 14

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

  • CrossRef
  • Google Scholar

• 15

  FritschA.DormieuxL.HellmichC. (2006). Porous polycrystals built up by uniformly and axisymmetrically oriented needles: Homogenization of elastic properties. Comptes Rendus Mécanique334, 151–157. 10.1016/j.crme.2006.01.008

  • CrossRef
  • Google Scholar

• 16

  FritschA.DormieuxL.HellmichC.SanahujaJ. (2009a). Mechanical behavior of hydroxyapatite biomaterials: An experimentally validated micromechanical model for elasticity and strength. J. Biomed. Mater. Res. Part A88A, 149–161. 10.1002/jbm.a.31727

  • CrossRef
  • Google Scholar

• 17

  FritschA.HellmichC.DormieuxL. (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

  • CrossRef
  • Google Scholar

• 18

  FritschA.HellmichC.YoungP. (2013). Micromechanics-derived scaling relations for poroelasticity and strength of brittle porous polycrystals. J. Appl. Mech.80, 020905. 10.1115/1.4007922

  • CrossRef
  • Google Scholar

• 19

  FritschA.HellmichC. (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

  • CrossRef
  • Google Scholar

• 20

  GrimalQ.RaumK.GerischA.LaugierP. (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

  • CrossRef
  • Google Scholar

• 21

  GuoW.HanF.JiangJ.XuW. (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

  • CrossRef
  • Google Scholar

• 22

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

  • CrossRef
  • Google Scholar

• 23

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

  • Google Scholar

• 24

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

  • CrossRef
  • Google Scholar

• 25

  HellmichC.BarthélémyJ.-F.DormieuxL. (2004). Mineral–collagen interactions in elasticity of bone ultrastructure – A continuum micromechanics approach. Eur. J. Mech. - A/Solids23, 783–810. 10.1016/j.euromechsol.2004.05.004

  • CrossRef
  • Google Scholar

• 26

  HellmichC.MangH. (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)

  • CrossRef
  • Google Scholar

• 27

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

  • CrossRef
  • Google Scholar

• 28

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

  • CrossRef
  • Google Scholar

• 29

  HofstetterK.HellmichC.EberhardsteinerJ. (2005). Development and experimental validation of a continuum micromechanics model for the elasticity of wood. Eur. J. Mech. - A/Solids24, 1030–1053. 10.1016/j.euromechsol.2005.05.006

  • CrossRef
  • Google Scholar

• 30

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

  • CrossRef
  • Google Scholar

• 31

  Jiménez SeguraN.PichlerB. L.HellmichC. (2023). Concentration tensors preserving elastic symmetry of multiphase composites. Mech. Mater.178, 104555. 10.1016/j.mechmat.2023.104555

  • CrossRef
  • Google Scholar

• 32

  Jiménez SeguraN.PichlerB. L.HellmichC. (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

  • CrossRef
  • Google Scholar

• 33

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

  • CrossRef
  • Google Scholar

• 34

  KönigsbergerM.HlobilM.DelsauteB.StaquetS.HellmichC.PichlerB. (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

  • CrossRef
  • Google Scholar

• 35

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

  • CrossRef
  • Google Scholar

• 36

  KrönerE. (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

• 37

  KrönerE. (1977). Bounds for effective elastic moduli of disordered materials. J. Mech. Phys. Solids25, 137–155. 10.1016/0022-5096(77)90009-6

  • CrossRef
  • Google Scholar

• 38

  LevinV. M. (1967). Thermal expansion coefficient of heterogeneous materials. Mekhanika Tverd. Tela2, 83–94.

  • Google Scholar

• 39

  LipinskiP.BarhdadiE. H.CherkaouiM. (2006). Micromechanical modelling of an arbitrary ellipsoidal multi-coated inclusion. Philos. Mag.86, 1305–1326. 10.1080/14786430500343868

  • CrossRef
  • Google Scholar

• 40

  LippmannB. A.SchwingerJ. (1950). Variational principles for scattering processes. I. Phys. Rev.79, 469–480. 10.1103/physrev.79.469

  • CrossRef
  • Google Scholar

• 41

  LucariniS.UpadhyayM.SeguradoJ. (2021). Fft based approaches in micromechanics: Fundamentals, methods and applications. Model. Simul. Mater. Sci. Eng.30, 023002. 10.1088/1361-651x/ac34e1

  • CrossRef
  • Google Scholar

• 42

  MoësN.CloirecM.CartraudP.RemacleJ. 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

  • CrossRef
  • Google Scholar

• 43

  MolinariA.El MoudenM. (1996). The problem of elastic inclusions at finite concentration. Int. J. Solids Struct.33, 3131–3150. 10.1016/0020-7683(95)00275-8

  • CrossRef
  • Google Scholar

• 44

  MoriT.TanakaK. (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

  • CrossRef
  • Google Scholar

• 45

  MoulinecH.SuquetP. (1994). A fast numerical method for computing the linear and nonlinear properties of composites. Comptes-Rendus l’Académie Sci. Série II318, 1417–1423.

  • Google Scholar

• 46

  MoulinecH.SuquetP. (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

  • CrossRef
  • Google Scholar

• 47

  PahrD. H.ZyssetP. 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

  • CrossRef
  • Google Scholar

• 48

  PichlerB.HellmichC.EberhardsteinerJ.WasserbauerJ.TermkhajornkitP.BarbaruloR.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. 10.1016/j.cemconres.2012.10.019

  • CrossRef
  • Google Scholar

• 49

  PichlerB.HellmichC. (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

  • CrossRef
  • Google Scholar

• 50

  RosenB. W.HashinZ. (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

  • CrossRef
  • Google Scholar

• 51

  SanahujaJ.DormieuxL.MeilleS.HellmichC.FritschA. (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

  • CrossRef
  • Google Scholar

• 52

  ScheinerS.SinibaldiR.PichlerB.KomlevV.RenghiniC.Vitale-BrovaroneC.et al (2009). Micromechanics of bone tissue-engineering scaffolds, based on resolution error-cleared computer tomography. Biomaterials30, 2411–2419. 10.1016/j.biomaterials.2008.12.048

  • CrossRef
  • Google Scholar

• 53

  SevostianovI.KachanovM. (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

  • CrossRef
  • Google Scholar

• 54

  TingT. C. T.LeeV.-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

  • CrossRef
  • Google Scholar

• 55

  TononF.PanE.AmadeiB. (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

  • CrossRef
  • Google Scholar

• 56

  TorquatoS. (1997). Effective stiffness tensor of composite media—I. Exact series expansions. J. Mech. Phys. Solids45, 1421–1448. 10.1016/s0022-5096(97)00019-7

  • CrossRef
  • Google Scholar

• 57

  VoigtW. (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

  • CrossRef
  • Google Scholar

• 58

  WangH.HellmichC.YuanY.MangH.PichlerB. (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

  • CrossRef
  • Google Scholar

• 59

  WhittakerE.RobinsonG. (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

• 60

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

  • CrossRef
  • Google Scholar

• 61

  WolframU.Peña FernándezM.McPheeS.SmithE.BeckR. J.ShephardJ. D.et al (2022). Multiscale mechanical consequences of ocean acidification for cold-water corals. Sci. Rep.12, 8052. 10.1038/s41598-022-11266-w

  • CrossRef
  • Google Scholar

• 62

  XieL.ZhangC.SladekJ.SladekV. (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

  • CrossRef
  • Google Scholar

• 63

  XuW.WuF.JiaoY.LiuM. (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

  • CrossRef
  • Google Scholar

• 64

  XuW.WuY.GouX. (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

  • CrossRef
  • Google Scholar

• 65

  XuW.WuY.JiaM. (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

  • CrossRef
  • Google Scholar

• 66

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

  • CrossRef
  • Google Scholar

• 67

  ZellerR.DederichsP. (1973). Elastic constants of polycrystals. Phys. Status Solidi (b)55, 831–842. 10.1002/pssb.2220550241

  • CrossRef
  • Google Scholar

• 68

  ZienkiewiczO. C.TaylorR. L.ZhuJ. Z. (2005). The finite element method: Its basis and fundamentals. Elsevier.

  • Google Scholar