椭球夹杂体复合材料有效性能预测：有效场方法与有限元模型的理论与应用对比研究

本综述系统比较了采用椭球夹杂体的解析模型与数值模型在预测复合材料有效性能方面的表现，重点评估了有效场方法（EFMs，含Mori-Tanaka、Maxwell等方法）与伪晶粒分解法（PGDM）的精度，揭示了EFMs在弹性、热膨胀和传导问题中与有限元（FE）模型高度一致，而PGDM仅适用于特定高纵横比夹杂体的弹性预测，为计算密集型FE模型提供了优选解析替代方案。

　　

Highlight

本研究比较了采用椭球夹杂体的解析模型与数值模型在预测复合材料有效性能方面的表现。重点评估了三类影响因素：均匀化问题类型（弹性、膨胀和传导）、材料相特性（各向同性、夹杂取向和纵横比）以及建模方法。通过参数化研究，将有效场方法（EFMs）——包括非相互作用、Mori-Tanaka和Maxwell方法——与伪晶粒分解法（PGDM）和数值有限元（FE）模型进行对比。研究涵盖纵横比0.2至5的椭球夹杂体，以及从完全随机到完全定向的取向分布。结果表明，均匀化问题类型对EFMs的预测精度影响不显著，Mori-Tanaka和Maxwell方法与FE模型在所有性能预测上表现出极好的一致性。PGDM仅对纵横比大于1的夹杂体复合材料的某些弹性性能（如E₁₁）能给出可靠结果。因此，Mori-Tanaka和Maxwell方法可作为计算密集型FE模型的最优解析替代方案。

Introduction

复合材料有效（均匀化）性能对于理解其宏观行为至关重要。通过确定这些性能，具有特定微结构的复合材料可被等效均匀材料替代，从而展现相同的整体响应。一旦有效性能得以估算，复杂几何可被简化，降低计算需求并提高模拟效率（Heller等，2016，2017）。这些有效性能可通过基于微力学的解析（Kachanov和Sevostianov，2018；Mityushev，2022；Mityushev等，2025）、数值（Gao和Chen，2021；Oskay，2015；Polyzos，Ravindranath等，2021）和半解析（Polyzos等，2022a；2023b；2023c）方法进行评估。这些方法广泛应用于估算有效刚度（Sevostianov和Kachanov，2002b）、热膨胀系数（Sevostianov，2012）以及热/电导率（Maxwell，1873；Sevostianov和Kachanov，2007a；Sevostianov等，2019）等多种复合材料。

Among the three, numerical and semi-analytical methods are generally considered more accurate than analytical techniques due to the higher geometric detail in modeling inhomogeneities (fibers, particles, etc.). These inhomogeneities can be obtained directly from imaging methods such as Scanning Electron Microscopy or micro-Computed Tomography, thus allowing for exact representation of the microstructure. However, such high-fidelity models are computationally intensive due to the fine meshes required (Nikolaou et al., 2024, Polyzos et al., 2022b).

To mitigate computational cost, idealized shapes – most commonly ellipsoids (and their limiting cases such as cylinders, spheres, and penny-shaped disks) – are often used in all three methodologies. The accuracy of such geometry simplifications depends on how well the ellipsoidal shapes approximate the real microstructure. Nonetheless, for many engineering applications, these approximations yield sufficiently accurate results and are widely adopted in industry (Sevostianov & Kachanov, 2013). It is important to note that ellipsoidal geometries are typically referred to as “inclusions”, whereas “inhomogeneities” more broadly refer to real-shape geometries.

Despite the simplification through ellipsoids, homogenization remains complex due to numerous influencing factors. These can be categorized into three main groups: the type of homogenization problem, characteristics of the material phases, and the modeling methodologies employed. The following section outlines these groups along with the contributions and novelties presented in this study for each. It is important to highlight that this study primarily aims to compare analytical and numerical models employing ellipsoidal inclusions and the ensuing discussion focuses on elements affecting such models.

Homogenization of a given property involves determining the volume-averaged response of a material by solving a partial differential equation (PDE), which varies depending on the physical problem under consideration (e.g., elasticity, conductivity). Accordingly, the governing PDE and the associated tensorial quantities differ in order: second-order PDEs arise in thermal, electrical, diffusion, and permeability problems (e.g., Fourier’s, Ohm’s, Fick’s, and Darcy’s laws), leading to second-rank tensors. In contrast, elasticity problems are governed by fourth-order PDEs and involve fourth-rank tensors, such as stiffness and compliance (Kachanov & Sevostianov, 2018).

Although the homogenization procedure is conceptually similar across different physical problems (Kachanov & Sevostianov, 2018), most studies in the literature focus on a single type of problem. In this work, both a second-order and a fourth-order PDE problem are addressed, and an effort is made to present the analytical formulation in a unified framework.

Several material factors influence the prediction of the models. These include (a) the isotropy of the considered material phases (the symmetry of the relevant tensors, e.g., elasticity tensor, conductivity tensor, etc.), (b) the orientation of the inclusions, and (c) the shape of the inclusions.

In most studies, both the matrix and inclusions are typically assumed to be isotropic (Jain et al., 2013, Pierard et al., 2004), primarily to simplify the governing equations (Jain et al., 2013). However, this assumption restricts the broader applicability of the analytical methods.

For analytical models, matrix isotropy is often a necessary condition. Closed-form solutions are available only for isotropic matrices with orientation-scattered inclusions (Tandon & Weng, 1984), or for transversely isotropic matrices with aligned inclusions (Sevostianov et al., 2005, Withers, 1989). Therefore, in cases where inclusions are randomly oriented, the matrix must be isotropic to ensure the analytical tractability of the problem.

From a mathematical standpoint, full isotropy in both phases enables the commutative property of tensor multiplication. That is, for two tensors A and B of any given rank, the operation AB = BA holds, where the inner product (e.g., A_imB_mj = C_ij for second-rank tensors and A_ijmnB_nmjk = C_ijkl for fourth-rank tensors) is implied. This property plays a critical role in the application of analytical methods. Consequently, comparisons involving isotropic inclusions should not be generalized to cases involving transversely isotropic inclusions (Sevostianov & Kachanov, 2014).

In this work, both isotropic and transversely isotropic inclusions are considered to systematically investigate the influence of phase isotropy.

Most studies consider only simplified orientation cases, such as fully aligned or completely random inclusions (Mishurova, Cabeza, Bruno, & Sevostianov, 2016). A notable exception is the work of Jain et al. (2013), which incorporates planar orientation scatter. However, comprehensive studies addressing three-dimensional (3D) orientation scatter are lacking in the literature.

The suitability of planar orientation scatter depends heavily on the geometry of the inclusions and the overall structure of the composite. For example, in thin composites reinforced with long fibers, where out-of-plane deviation is minimal, planar modeling may be sufficient (Bauer & B?hlke, 2022). In contrast, 3D-printed composites with short fibers typically exhibit significant orientation scatter in all directions (Yu, Hwang, Hwang, & Hong, 2019), making full 3D orientation modeling necessary. Accordingly, this work investigates various 3D orientation scatter scenarios to evaluate their influence on model predictions.

Analytically, orientation scatter is treated through the integration of a given tensor over all possible orientations, a process referred to as “orientation averaging” (Advani & Tucker, 1987). This integral can be approximated by discretizing the orientation space and summing the tensor, rotated into a finite set of directions (Desrumaux, Meraghni, & Benzeggagh, 2001).

However, this discretization introduces a fundamental physical inconsistency (B?hlke & Bertram, 2001): The orientation averaging process in dual problems – such as when the driving field¹ in an elasticity problem is either a displacement or a traction – can lead to different model predictions.

This inconsistency becomes apparent in the discretized form of the problem. Consider, for example, the stiffness tensor C and the compliance tensor S in their local orientations. When applying the Voigt and Reuss averages – also referred to as the arithmetic and harmonic means, respectively (B?hlke & Bertram, 2001) – the following relation is observed: C^eff = Σ_i ?ⁱ Cⁱ ≠ Σ_i ?ⁱ (Sⁱ)^-1 = (S^eff)^-1 where the superscripts “eff” and “i” denote the effective tensor and the tensor of phase i, respectively, and ? represents the volume fraction. The notation (Sⁱ)^-1 refers to the inverse of the compliance tensor. Here, the inverse of a tensor A is defined by the relation AA^-1 = I, with I being the unit tensor of the same rank as A. Although, physically, C^eff = (S^eff)^-1, the two averaging approaches yield different results (see also the discussion in B?hlke and Bertram (2001)).

In the literature, orientation averaging is typically performed for only one of the dual tensors—for example, the stiffness tensor in elasticity problems (Desrumaux et al., 2001, Jain et al., 2013). However, it is of interest to examine whether both dual formulations can be consistently applied in orientation averaging, and this is investigated in the present work.

The shape of the inclusions is typically characterized by the aspect ratios ξ₁ = r₁/r₂ and ξ₂ = r₁/r₃ of the ellipsoid, where r_i denotes the semi-axes of the ellipsoid. In literature ellipsoids of equal aspect ratios (r₂ = r₃) are examined (Kachanov & Sevostianov, 2018) and are mentioned as spheroids. Spheroids with aspect ratios greater than one are termed prolate spheroids, while those with aspect ratios less than one are oblate. A spheroid with an aspect ratio of one corresponds to a sphere.

The aspect ratio of a spheroid plays a significant role in influencing model predictions. In Jain et al. (2013), an aspect ratio of 3 was used due to challenges in meshing at higher ratios. However, as shown in Sevostianov et al. (2005), aspect ratios between 0.1 and 10 account for approximately 90% of the variability caused by changes in the spheroid’s shape. This sensitivity follows a logarithmic scale, with around 80% of the variability captured within the narrower range of 0.2 to 5. In the present study, aspect ratios from 0.2 to 5 were selected and represent the upper boundary of feasible meshing while allowing for a meaningful assessment of variability attributable to aspect ratio. Nevertheless, the results are limited to these aspect ratios and studies outside these bounds should be also conducted.

The modeling methodologies examined in this work include both analytical and numerical approaches. A brief overview of each is provided below.

Among the various analytical methods available for homogenizing material properties, two are particularly relevant when spheroidal inclusions with orientation averaging are considered: the effective field methods (EFMs) and the pseudo grain decomposition method (PGDM). Note that EFMs are also commonly referred to as mean-field methods.

EFMs account for inclusion interactions by treating each inclusion as isolated, subjected to an effective field that represents the influence of the surrounding medium. The specific form of the effective field depends on the chosen method. Since each inclusion is modeled independently, it can be arbitrarily rotated to account for different orientations.

Among the many EFMs available – such as the Non-Interaction, Mori–Tanaka, Maxwell, Self-Consistent, Kanaun–Levin, and differential schemes (Sevostianov & Kachanov, 2014) – only the Non-Interaction, Mori–Tanaka, and Maxwell methods possess explicit analytical solutions (Sevostianov, 2012). Among these, the Mori–Tanaka method is the most commonly employed. However, in certain material configurations, it can produce physically inadmissible results, such as non-symmetric homogenized tensors (Sevostianov & Kachanov, 2014). To remedy this, a symmetrized version of the Mori–Tanaka method was proposed in Sevostianov and Kachanov (2014), though it has yet to be validated against numerical models.

Pierard et al. (2004) proposed a different way to circumvent the mathematical problems of the Mori–Tanaka method. In their work, they consider a number of “pseudo-grains”, defined as a two-phase composite consisting of inclusions having the same orientation and aspect ratio. The pseudo-grains are, thus, unidirectional (UD) with respect to their local orientation system and can be homogenized as a UD composite. The Mori–Tanaka method can be then used for the homogenization of the UD composite and afterwards the resulting tensor can be rotated. Since for UD materials, the Mori–Tanaka method leads to symmetric tensors, the PGDM can be used without physical inconsistencies. Due to its applicability, the PGDM has found increased application and has been used in several works (see Jain et al., 2013) and has been included in commercial software such as ABAQUS (Smith, 2009) and DIGIMAT (HEXAGON, 2023).

Despite the fact that the PGDM eliminated the issues of the Mori–Tanaka method, it introduced additional approximations with regard to the interactions between the inclusions. In particular, the interactions are limited to the ones inside each grain which implies that the interactions refer to these of aligned inhomogeneities of the same kind. Physically this issue is non-trivial but in Doghri and Tinel (2005) it has been demonstrated that in general the predictions of the PGDM are in good agreement with numerical models. However, the models considered in Doghri and Tinel (2005) referred to composites with isotropic inclusions and 2D cases which raises questions related to the applicability of the PGDM for 3D cases.

Considering the above discussion, in this work three EFMs are considered (the Non-Interaction, the Mori–Tanaka, the Maxwell methods) and are compared to the PGDM and to 3D numerical models containing both isotropic and transversely isotropic inclusions to better understand the influence of the above assumptions.

Numerical models are influenced by several interrelated factors, including the spatial arrangement of inclusions within the representative volume element (RVE), the choice of boundary conditions, the RVE size, and the driving field.

The positions of inclusions inside the RVE can be either random or periodic, often referred to as random and periodic packings. It is important to distinguish between the positioning of inclusions within the volume of the RVE and their relation to the RVE boundaries. Random packings, where inclusions are randomly distributed throughout the RVE volume, may be implemented with either periodic or non-periodic boundary conditions (see discussions in Polyzos et al. (2022b)). In contrast, periodic packings are constructed so that inclusions repeat periodically both within the RVE volume and at its boundaries. Fig. 1 illustrates the difference between a random packing with periodic boundaries and a fully periodic packing.

Positioning inclusions periodically within the RVE volume is generally straightforward to implement, as the inclusion locations are predetermined. However, periodic packings have several drawbacks. First, they can introduce errors since inclusions in real materials are rarely arranged in a perfectly periodic manner (Nguyen, Béchet, Geuzaine, & Noels, 2012). Second, periodic packings cannot capture orientation scattering and are typically limited to symmetric shapes such as cylinders, spheres, or aligned spheroids.

For these reasons, most studies in the literature prefer random packings (Jain et al., 2013), which better reflect realistic microstructures. Nevertheless, generating random packings is more challenging because inclusion positions are not predefined, often requiring specialized algorithms (Drach, Kuksenko and Sevostianov, 2016).

In this work, both periodic and random packings are employed for cases involving aligned geometries to facilitate comparison. For scenarios involving orientation scattering, only random packings are used.

The position of the inclusions in the RVEs is also related to the boundary conditions. Periodic boundary conditions require the nodes on opposite surfaces of the RVE to be also periodically positioned (Xia et al., 2003, Xia et al., 2006). Therefore, periodic boundary conditions also necessitate the inclusions to be periodically positioned. In practice, however, there are several techniques that can be used for the imposition of periodic boundary conditions to RVEs with boundary nodes that are not periodic (see for example (Nguyen et al., 2012, Ouchetto et al., 2018) for a comparison between techniques on 2D structures).

Another type of boundary conditions that do not have these periodicity requirements is the uniform boundary conditions. However, such conditions can lead to “over-constrained” RVEs and to high errors in the predictions of the effective properties (Hollister and Kikuchi, 1992, Nguyen et al., 2012) and their imposition should be carefully examined.

The size of the RVEs must also be regulated accordingly to avoid the influence of the edges. For elastic problems, it has been shown that RVEs with 20–30 inclusions lead to convergence (Polyzos, Polyzos, Van Hemelrijck and Pyl, 2023) (over several realizations ²) for RVEs with periodic boundary conditions. For the case of uniform boundary conditions the number that leads to convergence is higher, but to the best of the authors’ knowledge no parametric studies have been conducted to quantify this number for the above mentioned 3D microstructures with spheroids. It is, therefore, of importance to quantify such numbers so that uniform boundary conditions can be imposed.

In this work, both periodic and uniform boundary conditions are imposed on RVEs with variable sizes to investigate the above issues.

Finally, the choice of the driving field can also influence the predictions. This driving field is, for example, a strain or traction in the framework of elasticity, and a temperature gradient or a heat flux in the framework of thermal conductivity (Sevostianov & Kachanov, 2009). In general, different driving fields can lead to differences in the predictions, which stands both for numerical (Polyzos, Van Hemelrijck and Pyl, 2021) and for analytical methodologies (Sevostianov & Kachanov, 2009). All possible driving fields are examined in this work to investigate their influence on the predictions of the effective properties.

Considering all the above, the primary goal of this study is to compare predictions from analytical and numerical models that utilize ellipsoidal inclusions to determine effective composite properties. The numerical models are considered here to represent accurate solutions to the homogenization problem, given their extensive validation in literature. In doing so, this work also aims to investigate various factors (detailed in Section 1.2) that influence the models and, consequently, the prediction accuracy of effective properties. To the best of the authors’ knowledge, such a comprehensive comparison has not yet been published in the open literature.

The structure of the study is as follows. Section 2 presents the analytical modeling framework, which provides background information and afterwards introduces a novel implementation of effective field methods (EFMs) for composites with orientation scattering using property contribution tensors. Section 3 details the development of the numerical models. Section 4 compares predictions from the analytical and numerical approaches. Finally, Section 5 summarizes the key conclusions drawn from the results.

Conclusions

The goal of this study was to compare the predictions between analytical and numerical models that use ellipsoids for effective properties and, simultaneously, to investigate factors (see Section 1.2) influencing the models and thus the prediction of the effective properties.

The main conclusions coming from the results (Section 4) regarding these factors can be summarized as follows.

The accuracy of predictions of the EFMs is not affected by the choice of homogenization problem. In general, both