# Contribution of microscale stochastic truss models to investigate the macroscale elasticity constants of porous ceramics

Thierry Canet, Gilles Dusserre, Thierry Cutard

► To cite this version:

Thierry Canet, Gilles Dusserre, Thierry Cutard. Contribution of microscale stochastic truss models to investigate the macroscale elasticity constants of porous ceramics. European Journal of Mechanics - A/Solids, 2025, 111, pp.105561. 10.1016/j.euromechsol.2024.105561. hal-04882712

HAL Id: hal-04882712

https://imt-mines-albi.hal.science/hal-04882712v1

Submitted on 13 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Full Length Article

# Contribution of microscale stochastic truss models to investigate the macroscale elasticity constants of porous ceramics

Thierry Canet, Gilles Dusserre* and Thierry Cutard

Institut Clément Ader (ICA), Université de Toulouse, CNRS, IMT Mines Albi, INSA, ISAE-SUPAERO, UPS, Campus Jarlard, F-81013, Albi, France

ARTICLE INFO

Keywords:

Porous ceramic
Truss
Ceramic matrix composite
Homogenization
Connector element
Modal analysis
Lattice

ABSTRACT

This paper deals with modelling the mechanical behaviour of silica-alumina open-cell porous ceramics obtained by viscous flow sintering. The modelling approach is based on the similarity of the material microstructure to a truss of sintering bridges connecting alumina particles. This makes it possible to use two-node elements, leading to a low computational cost. The method includes the building of a random packing of spheres (alumina particles) and the setup of connections between their centres (silica sintering bridges). An equivalent stiffness is then assigned to each bridge, based on the material parameters. The macroscale elasticity constants have been derived from the natural frequencies of such micrometric cylindrical volume elements made of thousands of particles. The reliability of the underlying assumptions is discussed and the dependence to the material parameters is emphasized. The method is suitable to handle more complex behaviours, which opens the door to fracture modelling.

1. Introduction

Porous ceramics can be classified in open-cell or closed-cell (Chen et al., 2021) and can be mainly manufactured by four conventional techniques (partial sintering, replica template, sacrificial template and direct foaming) and by additive manufacturing (Ahmad et al., 2014). The matrix of oxide/oxide ceramic matrix composites (CMCs) is a porous ceramic usually obtained by partial sintering of ceramic powder (Zok and Levi, 2001), often infiltrated into the fibrous reinforcement by slurry infiltration. The mechanical properties of the matrix, which are very sensitive to the porosity ratio, cf. (Marchelli, 2011) and references within, are key parameters of the composite behaviour and properties. Indeed, these materials belong to the family of weak matrix composites, whose non-brittle behaviour is based on matrix damaging to avoid or to delay fibre failure. The range of porosity ratio leading to a non-brittle behaviour without a prohibitive drop of mechanical properties, is narrow, between 30 and 40 % according to (Zok and Levi, 2001) and references within.

Matrix cracking and fibre-matrix debonding are the first damages in the scenarios proposed for such materials (Ben Ramdane et al., 2017; Guel, 2018; Koch et al., 2008; Marin, 2021) and play a major role in the initiation of delamination. Moreover, part of these damages already exists in the initial state of the composite due to sintering shrinkage and

mismatch of thermal expansion coefficients. It is therefore mandatory to investigate the mechanical behaviour of the matrix, which is a random assembly of particles of diameter typically about 0.2 µm, connected by narrow sintering necks. At the submicronic scale, matrix cracking results from the correlated failure of a large number of sintering necks of typically several tens of nanometres in diameter. This length-scale is about eight orders of magnitude smaller than the CMC part, which implies a broad range of experimental and numerical methods to investigate these phenomena. The scope of the present study is related to a new numerical method to investigate the microscale behaviour of the matrix.

The present paper focuses on a family of porous alumina-silica ceramics (cf. Fig. 1) formulated to act as the matrix of an oxide-oxide ceramic composite (Boutenel, 2020). The composite is obtained by slurry infiltration of the fibrous reinforcement with an aqueous suspension of dispersed alumina and silica particles. The suspension is then dried and the composite sintered to obtain the desired matrix porosity ratio. The sintering process mainly involves the viscous flow of the amorphous silica. The silica acts as a binder, and its amount must be optimized by considering its effect on the suspension rheology, on the densification kinetic and on the properties of the resulting matrix. This paper aims to contribute to more deeply understand the relationship between the latter and the material composition.

To this end, the work reported in this paper is based on the

<footnote>* Corresponding author. IMT Mines Albi, Institut Clément Ader, Campus Jarlard, route de Teillet, CT cedex 9, Albi, 81013, France.</footnote>

<footnote>E-mail address: Gilles.Dusserre@mines-albi.fr (G. Dusserre).</footnote>
consideration that the microstructure of such a material is an open-cell porous ceramic, which may be described as a truss structure. A large number of modelling tools are therefore available, allowing to consider a large number of particles at a limited computational cost.

A truss structure can be described as a set of members connected together at nodes, able to withstand high loads with an amount of material as low as possible and therefore widespread in construction and building. There has been renewed interest in such structures in recent years, related to the improvement of many additive manufacturing techniques (Qi et al., 2019, 2022; Aghajani et al., 2024; Khoda et al., 2021; Pei et al., 2022), making it possible to build lattice truss meso-structures, resulting in materials with tailorable or unusual properties such as anisotropy (Xu et al., 2016) or Poisson's ratio (Shaat and Wagih, 2020; Wang et al., 2022), or mimicking natural materials (Bhate et al., 2019; Ramirez-Chavez et al., 2022; Mohammadi et al., 2024; Li et al., 2022). Various strategies for optimizing the local non-uniform properties may be carried out, either focusing on the members section (Wang et al., 2018, 2020; Thillaithevan et al., 2021), or on the architecture (Garnier et al., 2022).

Several classifications of so-called lattice structures have been proposed (Bhate et al., 2019; Pei et al., 2022; Ramirez-Chavez et al., 2022), based on the periodicity of a repeating unit cell and its nature (strut- or shell-based or both) and the arrangement of the members within it. Indeed, these are the three main features that drive the emerging properties of the whole structure.

Most of the recent studies investigate periodic patterns, but a few are dedicated to aperiodic networks (Imediegwu et al., 2023). In this paper, a stochastic distribution of particles is considered, which is also the case in foams (Ashby, 2006; Belardi et al., 2024), granular (Li et al., 2022) or natural materials (Donaldson et al., 2008).

The arrangement of the members in the unit cell defines the connectivity of the lattice, which, together with the behaviour of the connections (pinned or rigid), drives the overall behaviour of the assembly, as early theorized by Maxwell (1864).

In the recent works, the members are generally slender struts of constant cross-section, which allows assumptions of bending- or stretching-dominated behaviours (Ashby, 2006). In the present work, the members are short beams of varying cross-section and made of two materials. Their behaviour is not likely to be dominated by a single solicitation.

If elastic behaviour is often the first studied property (Mohammadi

et al., 2024; Shaat and Wagih, 2020; Wang et al., 2022), many other behaviours are considered in such materials, as strength (Qi et al., 2022), anisotropy (Wang et al., 2022), plasticity, buckling or fracture (Ashby, 2006).

The basic principles of truss models were historically developed for stress analysis in structures (Ibrahimbegovic and Mejia-Nava, 2023; Maxwell, 1864). Recent works on periodic lattice meso-structures more often rely to the assessment of emerging macroscale properties by homogenization, using asymptotic (Cheng et al., 2019; Mohammadi et al., 2024; Qi et al., 2019; Wang et al., 2022) or numerical methods. A wide range of analytical (Ashby, 2006; Shaat and Wagih, 2020; Wang et al., 2022) or numerical models, mainly based on finite element methods, are described in the literature. Among the numerical approaches, slender struts are simulated by beam elements with constant cross-section (Maurizi et al., 2022; Mohammadi et al., 2024) or not (Belardi et al., 2024), continuum solid elements (Belardi et al., 2024; Donaldson et al., 2008; Meza et al., 2017; Qi et al., 2022), or a combination of both (Portela et al., 2018). The modelling strategy depends on the number of members required in the simulation and the richness of the expected results.

Other methods developed for granular materials may also be of interest. For instance, Discrete Element Method is based on explicit integration of interactions between particles in many-body problems and may be applied to model the studied material (Asadi et al., 2023; Li et al., 2022). However, this type of modelling requires to address transient dynamics.

Another family of models uses two-node elements, initially to simulate bolted joints (Guidault and Marie-France Soulé de Lafont, 2015) and also to account for the elasticity and fracture behaviour of granular materials (Asahina et al., 2017). In the case of dense granular materials, these methods succeed in modelling the macroscale elastic behaviour only at the cost of a tedious calibration. However, this kind of elements turns out to be well suited to model truss microstructures with non-slender members, and their low computational cost make it possible to consider a large number of members.

This work proposes to assimilate the random open-cell porous microstructure to a truss of short members that correspond to the sintering silica necks between two alumina particles. One half of both particles are considered in the members behaviour, too. The use of two-node elements to mimic the behaviour of the members allows to build a microscale model involving several thousands of particles. The aim of

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Silica and alumina particles after drying of the suspension before sintering (a.) and resulting microstructure after sintering at 1250 °C for 65 min (b.) (scanning electron microscopy). The solid fraction contains 59 vol% Al₂O₃ and 41 vol% SiO₂, the porosity ratio is 29.5 ± 0.5%. In c., d. and e. (same magnification) the dotted lines highlight sintering silica bridges, hourglass-shaped, between alumina particles.</figcaption></figure>
the paper is to assess the reliability of the model and the underlying assumptions to simulate the elastic behaviour of the material at the macroscale, as a function of the material parameters. It is a first step toward the optimisation of the porous material to act as the matrix of a ceramic matrix composite.

## 2. Material and methods

The studied material results from the partial sintering (from a few minutes to several hours between 1150 and 1300 °C) of a green body consisting of nanometric amorphous silica particles (density $2.2 \text{ g cm}^{-3}$) and submicronic α-alumina particles (density $3.96 \text{ g cm}^{-3}$). The green body is obtained by drying the aqueous suspension prepared with a nanometric silica suspension (Ludox® AS-40, Grace Davison, USA), an α-alumina powder (AKP-50, Sumimoto Chemical, Japan) and an ammonium polymethacrylate solution (Darvan® C-N, Vanderbilt Minerals LLC, USA) acting as dispersant for alumina particles. The total solid fraction is fitted by adding distilled water and the suspension is deagglomerated by ultrasonic treatment. The detailed procedure is available in (Boutenel et al., 2021). Initially, the silica particles are almost perfect spheres of typical diameter 24 nm and the alumina particles are more likely ellipsoids of median diameter 206 nm, Fig. 1a.

During the sintering step, the silica particles enter a viscous state which allows them to coalesce and build viscous bridges that connect alumina particles. Under the effect of surface tension, an attractive force arises between the alumina particles connected to such a bridge, leading to densification during the sintering (Margarida et al., 2008). During this process, the silica remains amorphous, as checked by X-ray diffraction, and the sintering mechanism is assumed to involve only viscous flow. The duration of the process is adjusted to match the required porosity ratio. The resulting material exhibits an open-cell microstructure with random arrangement of alumina particles, linked by solid silica bridges of various length and cross-section (cf. Fig. 1b). Due to the high porosity ratio required, the densification process is stopped at its early stage and most of the bridges are hourglass-shaped.

For a porosity ratio ranging between 28 and 39 % and a ratio between 59 and 92 vol% of alumina in the solid fraction, it has been shown in (Boutenel, 2020), that the porosity is the only parameter, at given composition, that drives both microstructure and rupture strength of the material, provided that silica remains amorphous after sintering. Therefore, only two parameters control the microstructure, for instance porosity ratio and another parameter related to the composition (alumina content, silica content or silica to alumina ratio).

The elasticity constants used in this paper for α-alumina and amorphous silica at room temperature are reported in Table 1 (De Faoite et al., 2012; Deschamps et al., 2014).

## 3. Theory and calculation

The model proposed in this paper places at the submicronic scale of the alumina particles and silica bridges. It is based on an analogy of the open-cell porous microstructure (Fig. 1b) with a truss whose nodes are the centres of the alumina particles, and the members coincide with the silica bridges. Each member is modelled thanks to a single two-node element. Abaqus connectors of bushing type, with 6 degrees of freedom, are used in this paper. Only six degrees of freedom are therefore required for each silica bridge, making it possible to consider a large

<caption>**Table 1**</caption>

Young's modulus and Poisson's ratio of α-alumina (De Faoite et al., 2012) and amorphous silica (Deschamps et al., 2014) at room temperature.

<table><thead><tr><th></th><th>Young's modulus (GPa)</th><th>Poisson's ratio</th></tr></thead><tbody><tr><td>α-alumina</td><td>352.8</td><td>0.208</td></tr><tr><td></td><td>71.5</td><td>0.176</td></tr><tr><td>amorphous silica</td><td></td><td></td></tr></tbody></table>

number of particles, typically several thousands, even with an office computer. In this paper, only linear elasticity is considered in order to assess the merit of the approach, but much more complex behaviours can be assigned to those elements (Guidault and Marie-France Soulé de Lafont, 2015).

Let's consider a sample of volume V of the material, containing a volume $V_{SiO_2}$ of silica and a volume $V_{Al_2O_3}$ of alumina. The residual volume is porosity. Regarding the specific microstructure of the studied material, parameterized by its porosity ratio, p, and silica to alumina volume ratio, $V_{SiO_2}/V_{Al_2O_3}$ (with $p + V_{SiO_2}/V + V_{Al_2O_3}/V = 1$), the following hypotheses are made:

- H1: The alumina particles are idealized as spheres of radius R = 103 nm. Particle size distribution spreads between 100 and 400 nm, cf. (Boutenel et al., 2021)

- H2: The silica bridges are assumed to contain the whole amount of silica, $V_{SiO_2}$, available as particles in the precursor suspension, i.e. the sintering process is sufficiently far advanced, so that the initial granular structure is completely lost. However, the sintering process is not long enough to allow the liquid bridges to merge, i.e. the material is in the pendular state (Weigert and Ripperger, 1999).

- H3: The member stiffness is calculated by idealizing the silica bridge geometry as the one of an axisymmetric liquid meniscus, supposed to be formed at high temperature during sintering (cf. Fig. 1b).

- H4: The alumina particles are positioned by applying an isotropic scaling to an initial random packing occupying a volume $V_0$. The scaling factor, $\alpha$, depends on the initial ($V_{Al_2O_3}/V_0$) and target

$$ (V_{Al_2O_3}/V) \text{ alumina volume fraction: } \alpha = \sqrt[3]{V/V_0} = \sqrt[3]{\frac{V_{Al_2O_3}}{V_0}} $$

### 3.1. Bridge geometry

As reported in (Boutenel, 2020), viscous sintering leads to a volumetric shrinkage of a few tens of percent as a result of the capillary forces acting on a pair of alumina particles through the viscous silica bridge, which enables to formulate the hypothesis H2. As a default hypothesis, it may be further assumed that the liquid bridge at the end of sintering is close to equilibrium. Accordingly, the geometry of the bridge in the porous ceramic is similar to a liquid meniscus (H3), with a given contact angle, $\delta$, uniform along all the three-phase contact lines. Finally, the hypothesis H1 implies that the meniscus is axisymmetric, and for the sake of simplicity, a toroidal approximation can be made (Fisher, 1926; Lian et al., 1993).

In such conditions, the geometry of the bridge can be fully described with only two parameters, for instance the separation distance between both particles, $\alpha$, and either the volume of the toroidal meniscus, $V_b$, or the so-called half-filling angle, $\beta$, Fig. 2. Indeed, $V_b$ and $\beta$ are related through Eqs. (1)–(4), which, according to (Weigert and Ripperger, 1999), within such assumptions, give the volume of a liquid bridge as a function of $\alpha$, $R$, $\beta$ and $\delta$.

Building the bridge geometry in section 3.2 is more convenient using $\alpha$ and $\beta$ as parameters, whereas the setup of the truss structure is based on the distance between the centres of the particles, $l = a + 2R$ (cf.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Idealized schematization of a silica bridge between two alumina particles, under the toroidal approximation (Weigert and Ripperger, 1999).</figcaption></figure>
section 3.3), and the volume of the toroidal meniscus, $V_b$ (cf. section 3.4). Both parametrizations are equivalent.

equivalent stiffness properties to the two-node element that acts as the bridge in the modelled truss-like microstructure. Thanks to the revolu-

$$ \frac{V_b}{2\pi} = (R_1^2 + b^2)R_1 \cos(\beta + \delta) - \frac{R_1^3 \cos^2(\beta + \delta)}{3} - b(R_1^2 \cos(\beta + \delta)\sin(\beta + \delta) + R_1^2(\frac{\pi}{2} - \beta - \delta)) - \frac{R^3}{3}(2 + \cos \beta)(1 - \cos \beta)^2 \quad \text{Eq. (1)} $$

$$ b = R_1 + R_2 \qquad \text{Eq. (2)} $$

$$ R_1 = \frac{2R(1 - \cos \beta) + a}{2 \cos(\beta + \delta)} \qquad \text{Eq. (3)} $$

$$ R_2 = R\sin\beta + R_1(\sin(\beta + \delta) - 1) \qquad \text{Eq. (4)} $$

The assumption of a uniform contact angle is very strong since the wetting phenomenon is a complex dynamic process, which depends on many parameters. Nevertheless, it leads to a bridge geometry roughly similar to the observations (Fig. 1b). The value of the contact angle has been set to 15°, based on the data reported in (Comeforo and Hursh, 1952), where the contact angle at 1200 °C between a molten soda-lime-silica glass and a dense 99% alumina substrate has been shown to quickly tend toward a value between 15° and 20° after a few tens of minutes, which is often the order of magnitude of the sintering step duration. The decrease is afterwards much slower, making the previous assumption possible.

Two physical limitations constrain the bridge geometry and ensure that the model keeps a physical meaning. The first one is based on a minimum surface free energy criterion of bridge stability at the molten state. According to (Lian et al., 1993), for a contact angle up to 40°, this criterion is equivalent to a good approximation to a direct relationship between the critical separation distance at break, $a_c$, and the bridge volume, $V_b$, Eq. (5), with $\delta$ expressed in radians. Eq. (5) is used to set a minimum allowable silica volume, $V_{b min}$, for each bridge likely to exist in the structure, Eq. (6).

$$ a_c = (1 + 0.5\delta)^{\frac{3}{2}} \sqrt{V_b} \qquad \text{Eq. (5)} $$

$$ V_{b min} = a^3 (1 + 0.5\delta)^{-3} \qquad \text{Eq. (6)} $$

The second limitation is related to the hourglass-shape of the bridges. It implies that $\beta < \beta_{max}$, with $\beta_{max} = \frac{\pi}{2} - \delta$, to avoid barrel-shape. At a given separation distance, the maximal volume of a bridge, of cylinder shape, $V_{b max}$, is thus obtained for $\beta = \beta_{max}$ and is given by Eq. (7).

$$ V_{b max} = \pi R^3 \left( \cos^2 \delta \left( 2(1 - \sin \delta) + \frac{a}{R} \right) - \frac{2}{3}(1 - \sin \delta)^2 (2 + \sin \delta) \right) \qquad \text{Eq. (7)} $$

For a given contact angle and particle radius, Eq. (6) and Eq. (7) provide an upper bound for the separation distance, $a_{max}$, that arises when the minimum volume given by Eq. (6) equals the maximum volume given by Eq. (7). With $R = 103$ nm and $\delta = 15°$, the higher separation distance likely to be built in the microstructure is $a_{max} = 238$ nm. As many as 1285 silica particles would be required to build such a bridge, longer than the average alumina particle diameter.

### 3.2. Bridge equivalent stiffness

The approximated bridge geometry described in the previous section is controlled with only two geometrical parameters, $a$ and $\beta$. A quasi-static finite element model has been used to compute the stiffness of a bridge as a function of those parameters, in order to assign some

tion symmetry of the bridge, the shear and bending properties are the same in all transverse directions. Therefore, only four loadings per couple $(a, \beta)$ are required, each related to one degree of freedom of the two-node element: tension, shear, torsion and bending.

The finite element model (Fig. 3a) includes two half-spheres that model the pair of alumina particles, and one silica bridge as described previously. The convenient elasticity properties given in Table 1 are assigned to each material, assumed to be homogeneous, isotropic and linear elastic. A typical mesh size of 10 nm is defined (cf. convergence study in appendix A and a comparison of several element types in appendix B). Depending on the parameters $a$ and $\beta$, 4163–14637 C3D10 continuum elements (respectively 6465 to 21417 nodes) are required. These tetrahedral elements allow an efficient automatic meshing of the complex geometry for a large number of $a$ and $\beta$ values, while second order interpolation is required for curved surfaces and stress concentration issues.

For each alumina particle, the three degrees of freedom of the nodes in the median plane perpendicular to the bridge axis are kinematically coupled to a reference point at the centre of the sphere, on which boundary conditions and loads are applied (Fig. 3b and Table 2). The total CPU time required to perform the four computations for a single bridge is typically of 50 s (laptop, CPU @1.6 GHz, 16 GB memory). The equivalent bridge stiffnesses to be assigned to the two-node elements are given in Table 2.

The set of possible couples $(a, \beta)$ is restricted by the aforementioned assumptions: $a_{max} < 238$ nm to ensure that the minimum volume of the bridge is lower than its maximum value, $\beta \le \frac{\pi}{2} - \delta$ to avoid barrel-shaped bridges. Moreover, the inequality $V_b > V_{b min}$, cf. Eq. (6), based on the stability of the viscous bridge, together with Eq. (1), provides a lower bound of $\beta$ for a given $a$. To avoid meshing issues, a 10% margin is applied to these criteria, and the minimum separation distance $a$ is set to 10 nm, which is lower than the radius of a silica particle.

Within the aforementioned bounds, the $(a, \beta)$ plane is mapped with 100 couples, with $a$ between 10 and 238 nm, and $\beta$ between the lower value allowed by $V_b > V_{b min}$ and $\frac{\pi}{2} - \delta$. For each couple, the four equivalent stiffnesses are computed using the previous finite element model and their values are fitted using Bezier surfaces of degree (6,6), cf. Fig. 4.

### 3.3. Setup of nodes location

The material is considered as a random open-cell porous microstructure. According to the hypothesis H4, the node location is obtained by applying an isotropic scaling to an initial random packing of 206 nm-diameter spheres mimicking the alumina particles. The scaling factor, $a$, controls the amount of alumina in the material.

The starting random packing consists in a set of $N$ nodes included in a volume of the desired shape collected into a large densely packed cluster of nodes, Fig. 5a. The cluster, that contains typically 70000 nodes, is built using a method similar to the so-called “3D clew approach” proposed by (Falco et al., 2017) to model particle reinforced composites. This method is based on building a cluster of spherical particles one by one, starting from a central one and including sequentially the next particle by maximizing the tangency between spheres. Building the initial cluster has a high computational cost, but once implemented, many volumes of any shape are rapidly collected (typically 5 s for 13000
<figure><img src="image_4.png" /><figcaption>Fig. 3. Finite element model of a bridge: a. geometry (a = 111.3 nm and β = 53.4°) and mesh, b. Mises stress fields for each elementary solicitation.</figcaption></figure>

<caption>**Table 2**</caption>

Boundary conditions and loads applied to the reference points in the bridge finite element model (cf. Fig. 3). *Ui*, *URi*, *Fi* and *Mi* refer respectively to the displacement, rotation, force and moment in direction *i*.

<table><thead><tr><th rowspan="2"></th><th colspan="2">Tension</th><th colspan="2">Torsion</th><th colspan="2">Bending</th><th colspan="2">Shear</th></tr><tr><th>RP1</th><th>RP2</th><th>RP1</th><th>RP2</th><th>RP1</th><th>RP2</th><th>RP1</th><th>RP2</th></tr></thead><tbody><tr><td>U1</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>0</td></tr><tr><td>U2</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td></tr><tr><td>U3</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>0</td></tr><tr><td>UR1</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>0</td></tr><tr><td>UR2</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>0</td></tr><tr><td>UR3</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>0</td></tr><tr><td>F1</td><td>-</td><td>F1</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>-</td></tr><tr><td>F2</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>F2</td></tr><tr><td>M1</td><td>-</td><td>0</td><td>-</td><td>M1</td><td>-</td><td>0</td><td>-</td><td>-</td></tr><tr><td>M3</td><td>-</td><td>0</td><td>-</td><td>0</td><td>-</td><td>M3</td><td>-</td><td>-</td></tr><tr><td>Stiffness</td><td colspan="2">K1 = F1/U1</td><td colspan="2">KR1 = M1/<br>UR1</td><td colspan="2">KR2 = KR3 = M3/<br>UR3</td><td colspan="2">K2 = K3 = F2/<br>U2</td></tr></tbody></table>

particles). According to the radial distribution function centred at the initial particle and plotted in Fig. 5b, it is checked that the resulting cluster doesn't exhibit any long-range order.

### 3.4. Setup of connections between nodes

Starting from the set of nodes collected in a large cluster and scaled to the desired amount of alumina (section 3.3), the separation distance between each pair of alumina particles is fixed and connections must be setup between nodes to act as the members of the truss-like structure, intended to simulate the silica bridges. Regarding the previous assumptions, some restrictions arise. The total volume of the bridges should be consistent with the amount of silica in the volume element, $V_{SiO_2}$. Moreover, separation distances higher than 238 nm are prohibited, cf. section 3.1, but it turns out that this condition is not actually restrictive. Finally, the bridge volume must lay between the allowed upper and lower bounds. Connections are setup using the following process:

- i) The first step consists in distributing $V_{SiO_2}$ among the N alumina particles. In this paper a uniform distribution has been used as a default assumption.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Equivalent tensile stiffness as a function of a and β, computed with the finite element model (dots) and interpolated by a Bezier surface (R²: coefficient of determination). The shape of some bridges is illustrated for instance (red: alumina particles; blue: silica bridge). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

- ii) Then, a first oversized set of connections is setup by pairing the nodes spaced at a distance lower than a given value. In order to account for the material parameters, this value was chosen empirically as a function of the critical separation distance of a bridge containing all the silica assigned to the particle. With a uniform distribution of the silica among the alumina particles and a uniform contact angle, the critical separation distance is the same for each particle and Eq.

(5) reduces to $a_c = (1 + 0.5\delta) \left(\frac{V_{SiO_2}}{N}\right)^{1/3}$. The results reported in the paper have been obtained with the following criterion: $l < 2 a_c$.

- iii) In a third step, the silica volume assigned to each particle, namely $\frac{V_{SiO_2}}{N}$, is evenly distributed between all the bridges connected to the considered particle.

- iv) The volume of each bridge is then compared to the maximum allowable value given by Eq. (7). If required, the volume in excess in a given bridge is uniformly distributed to the bridges linked to both concerned particles, unless their volume has already reached the maximum value. This step is resumed if necessary until all the bridges match the maximum volume criterion.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Selection of a 12 µm long cylindrical volume element of 13000 particles from a large (12 × 8 × 8 µm³) cluster of 70186 densely packed particles (a.) and its radial distribution function (b.).</figcaption></figure>

- v) Afterwards, the minimum volume criterion, Eq. (6), is checked. The volume of any bridge that doesn't match the criterion is temporarily set to zero.

- vi) The whole amount of silica removed in the previous step is once again evenly distributed over all the particles, and the process resumes at the third step until the residual volume reaches a given value, 1% in the results reported in this paper.

- vii) A last iteration is performed from the third step and the bridges that doesn't match the minimum volume criterion are definitively removed.

Finally, it is checked that all the particles are connected to a single network, which is the case for all the results reported in this paper. It is noteworthy that issues related to unconnected networks in the volume element may arise only from 70% of porosity, which is out of scope here.

The resulting truss-like network of bridges is fully characterized by the length, *l*, of each member and the volume of silica assigned to each member, $V_b$. For each member, a couple of parameters $(a, \beta)$ is deduced ($a = l - 2R$ and afterwards $\beta$ from Eqs. (1)–(4)), and the corresponding interpolated values of the four stiffnesses are assigned to the two-node CONN3D2 element that models the member (bushing type connector, linear elastic behaviour, no coupling). This procedure, including interpolation and assignment, is very fast, typically 1 µs, whereas 50 s are required to compute the 4 stiffnesses, which allows to consider sets of a large number of particles, typically several thousands in this paper.

### 3.5. Homogenization procedure

The truss-like structure, based on two-node connector elements, built in the previous section may be used in many types of procedure in Abaqus/Standard and Abaqus/Explicit, including but not limited to quasi-static, dynamic and linear perturbation (for instance buckling or natural frequencies extraction) analyses.

The purpose of this paper is to assess the effect of material parameters (porosity ratio, volume fractions of alumina and silica) on the macroscale linear elasticity constants of the porous ceramic, assumed to be isotropic thanks to the randomness of the particle arrangement. In order to avoid a tedious implementation of periodic boundary conditions on a random arrangement of nodes, it is proposed here to assess the elasticity constants from the natural frequencies of a truss-like structure, encased in a cylinder, acting as a representative volume element (cf. Fig. 5a).

#### 3.5.1. Natural frequencies extraction

The natural frequencies of the representative volume element are extracted using the Lanczos eigenvalue solver in a linear frequency

analysis, which requires nodal masses to be assigned to the nodes. Accounting for the mass and rotary inertia of alumina particles, assumed to be spherical, is straightforward since their centre of gravity coincides with a node of the truss-like structure, which supports the mass in the frequency analysis. It is trickier regarding the silica mass, since their centre of mass are by definition between two nodes.

As a basic assumption, if the members connected at a given node are similar and oriented almost evenly in all directions, the centre of mass of all of these members should be close to the node. In accordance, the half mass of each bridge is assigned to both particles to which it is connected.

Four assumptions are compared in section 4.3, regarding the rotary inertia of silica bridges.

- I1: The rotary inertia of both alumina and silica is neglected.

- I2: The mass of alumina particle and the half-mass of the bridges are supposed uniformly distributed in the volume of the particle, leading to an underestimation of the rotary inertia.

- I3: The half-mass of the bridges is supposed to be uniformly spread around the particle, leading to a weaker underestimation of the rotary inertia than I2.

- I4: The half silica bridges are modelled as cylinders of the same mass, length and axis, and their contributions to the rotary inertia are calculated at the centre of the particle, leading to a very slight overestimation of the rotary inertia.

The mass of an alumina particle is of the order of $10^{-23}$ kg and the stiffness of a bridge, for instance in tension, is of the order of $10^3$ N.m$^{-1}$, leading to natural frequency of the order $10^{13}$ Hz. In order to avoid the solver to handle such small or large numbers, possibly leading to rounding errors, the mass unit chosen is mg. Since N and mm are used as force and length units respectively, the time unit is affected accordingly (for instance $10^{-\frac{9}{2}}$ s). As a consequence, the whole process does not affect the final results, but only some intermediate computed values.

The CPU time required to process the whole method, including the truss setup and the extraction of the first ten natural frequencies, for a truss-like volume element of 1000–13500 particles is respectively between 5 min and 2 h.

#### 3.5.2. Macroscale elasticity constants

The natural vibration modes of a cylindrical representative volume element are extracted as described in the previous sections. Among them, the first order torsion mode and both first order bending modes are selected for homogenization purpose. Their eigenfrequencies are respectively $f_t$, $f_{b1}$ and $f_{b2}$, with $f_{b1} < f_{b2}$. Fig. 6 shows that the shapes of these modes are very close to the expected kinematics for pure first order bending and torsion modes.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Bending (a.) and torsion (b.) mode shapes of the connector-based model (cf. section 4.2 for a description of the volume element).</figcaption></figure>

The macroscale shear modulus, G, is deduced from $f_t$ using Eq. (8). Then three couples of macroscale Young's modulus and Poisson's ratio, ($E_i, \nu_i$), with $i = 1, 2, m$, are deduced from $f_{bi}$ using Eqs. 9–12 (ASTM C215, 2020; ASTM E1875, 2021; ASTM E1876, 20222). These equations are originally involved in the experimental measurement of elasticity constants from cylindrical samples. $E_m$ and $\nu_m$ are respectively considered as the Young's modulus and Poisson's ratio of the equivalent homogeneous isotropic linear elastic material. L, D and m are respectively the length, diameter and mass of the cylindrical sample. In order to limit the surface effect, short cylinders of aspect ratio $L/D = 3$ are used in this paper.

$$G = \frac{16}{\pi} m \frac{L}{D^2} f_t^2 \qquad \text{Eq. (8)}$$

$$E_i = 1.6067m \frac{L^3}{D^4} f_{bi}^2 T_i \qquad \text{Eq. (9)}$$

$$T_i = 1 + 4.939(1 + 0.0752\nu_i + 0.8109\nu_i^2)\left(\frac{D}{L}\right)^2 - 0.4883\left(\frac{D}{L}\right)^4 - \left[ \frac{4.6911(1 + 0.2023\nu_i + 2.173\nu_i^2)\left(\frac{D}{L}\right)^4}{1.000 + 4.7541(1 + 0.1408\nu_i + 1.536\nu_i^2)\left(\frac{D}{L}\right)^2} \right] \qquad \text{Eq. (10)}$$

$$\nu_i = \frac{E_i}{2G} - 1 \qquad \text{Eq. (11)}$$

$$f_{bm} = \frac{1}{2}(f_{b1} + f_{b2}) \qquad \text{Eq. (12)}$$

## 4. Results and discussion

In this section, the homogenized elasticity constants of the porous ceramic are computed with the method described in section 3, as a function of the material composition, namely porosity ratio, p, and silica to alumina volume ratio, $V_{SiO_2}/V_{Al_2O_3}$. Prior to that, the representativeness of a volume element is assessed. Then, the relevance of the method is checked i) by comparing the frequencies of a connector-based model with those of a continuum with corresponding homogeneous properties, and ii) by studying the effect of two assumptions. The first one is related to the assignment of silica rotary inertia (cf. section 3.5.1) and the second one to the computation of bridge equivalent stiffness (cf. sections 3.1 and 3.2).

### 4.1. Representativeness of the volume element

In order to assess the representativeness of a volume element, homogenizations of several cylindrical truss-like microstructures of the same aspect ratio $L/D = 3$ have been performed, with a length ranging from 5 µm to 12 µm before scaling (respectively 1007 to 13445 particles). The porosity ratio, p, and the silica to alumina volume ratio, $V_{SiO_2}/V_{Al_2O_3}$, are respectively set to 0.4 and 0.695. The corresponding homogenized shear and Young's moduli are reported as a function of the number of particles in Fig. 7a, and Poisson's ratio in Fig. 7b.

Fig. 7 shows that, above 5687 particles ($L = 9$ µm), the elasticity moduli are within a range of ±2% around the value obtained with 13445 particles, and the Poisson's ratio within a range of ±0.01. The Young's moduli, $E_1$ and $E_2$, and Poisson's ratios, $\nu_1$ and $\nu_2$, corresponding to each of both bending eigenfrequencies, are within a range of respectively ±4% and ±0.03. Between 10360 ($L = 11$ µm) and 13450 particles ($L = 12$ µm), the variation is only of 0.5% (less than 0.2 GPa) for both shear and Young's moduli, and of 0.01 for the Poisson's ratio. In these cases, $E_1$ and $E_2$ are both within the range of ±2%, which confirms

<figure><img src="image_8.png" /><figcaption>Fig. 7. Macroscale elasticity constants of cylindrical volume elements with aspect ratio $L/D = 3$ as a function of the number of particles: a. shear and Young's moduli and b. Poisson's ratio. Dot lines display a range of ±2% and ±4% (a.), and ±0.01 and ±0.03 (b.) around the value for the larger number of particles (13445).</figcaption></figure>
that the volume element is quasi isotropic.

As a result, depending on the required accuracy, a volume element of 5687 (L = 9 µm) or 10360 particles (L = 11 µm) may be used. The corresponding computation time is respectively below 1000 s and 2000 s (laptop, CPU @1.6 GHz, 16 GB memory). In the following sections, volume elements of length of 10, 11 or 12 µm are considered.

## 4.2. Comparison of connector network with homogenized continuum

In order to check the relevance of the homogenization process, the natural vibration modes of a connector-based volume element, as described in section 4, are compared to the vibration modes of a cylindrical continuum of the same dimensions and exhibiting the corresponding homogenized elasticity constants and the same apparent density. The same material characteristics as in the previous section are considered, i.e. $p = 0.4$ and $V_{SiO_2}/V_{Al_2O_3} = 0.695$.

The homogenization procedure is applied to a 12 µm length cylindrical volume element (13445 particles) with $L/D = 3$, providing a Young's modulus of 36.2 GPa, a shear modulus of 15.5 GPa and a Poisson's ratio of 0.17. The natural frequencies of the corresponding continuum have been extracted by a finite element analysis (the same unit system as in section 3.5.1. is applied). The cylinder is meshed with 13600 C3D8 elements (15249 nodes) about 0.25 µm in size. The resulting relative deviation between the squared natural frequencies of each of the nine first modes, taking the continuum as reference, are reported in Fig. 8.

As shown is section 3.5.2, the shape of the three first modes used for homogenization purpose are very close to the ones of the continuum, and obviously, the deviation between the squared natural frequencies is low, less than 1.2%. The difference in the mode shape is higher for the

six next modes, especially for the first order tension and second order torsion, where the expected axial symmetry is broken with the connector-based model. Indeed, high order modes are much more sensitive to local heterogeneities, that may result in mixed modes.

Actually, comparing the frequencies of the so-called tension mode is not really meaningful, because they are very close to the frequencies of the second-order bending mode. It turns out that the connector-based model is unable to exhibit a pure tension mode and instead displays a mixed mode involving tension and second order bending, as supported by the frequencies of each mode, respectively of 167, 171, 174 and 175 MHz for tension with continuum elements, pseudo-tension with two-node elements and second-order bending with two-node and continuum elements.

Nevertheless, the deviation is lower than 3.3% for the six first modes and lower than 6% for all modes, except the second order torsion (12%). These results show that the connector-based model, with a volume element as small as a ten of micrometres, is able to account for the macroscale behaviour of the equivalent continuum, while providing an insight into local effects related to heterogeneity.

## 4.3. Sensitivity of homogenized properties to mass distribution

The four assumptions given in section 3.5.1, regarding the rotary inertia of silica bridges, are compared in this section, based on the resulting macroscale elasticity constants of a cylindrical volume element (length $L = 10$ µm, aspect ratio $L/D = 3$, and silica to alumina volume ratio $V_{SiO_2}/V_{Al_2O_3} = 0.695$). Taking the assumption 14 (bridges are modelled as cylinders) as a reference, the relative deviation on elasticity moduli and the deviation on the Poisson's ratio are reported in Fig. 9, for porosity ratios between 0.35 and 0.55. These results show that the rotary

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of the natural vibration modes of a connector-based model with its equivalent continuum and relative deviation between the squared natural frequencies. (a. relative deviation between the squared natural frequencies, b. natural modes shape obtained with connector elements and c. natural modes shape obtained with continuum elements).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Relative (a. and b.) or absolute (c.) deviation on the macroscale elasticity constants: (a. Young's modulus, b. shear modulus and c. Poisson's ratio, cf. notations in section 3.5.2) computed with several assumptions for rotary inertia (I1, I2 and I3, cf. section 3.5.1) taking assumption I4 as a reference.</figcaption></figure>

inertia has no significant influence on the elasticity moduli, with a relative deviation lower than 0.3%, except for assumption I1 and porosity ratio of 0.35, and affects slightly more the shear modulus than the Young's modulus. The deviation on Poisson's ratio remains lower than 0.003 in absolute value for assumptions I1, I2 and I3. Assumption I2 provides the lower deviations, with relative deviations on elasticity moduli below 0.07% and deviations on Poisson's ratio below $6 \cdot 10^{-3}$. Moreover, computing the rotary inertia is slightly faster with assumption I2 than with assumption I4 (about 5 s less for 10000 particles).

As a result, the assumption on the rotary inertia of the silica bridges,

which is a limitation of using two-nodes elements to model the microstructure, doesn't affect significantly the computed macroscale elasticity constants. In the following sections, assumption I2 has been used by default.

### 4.4. Reliability of member stiffness computation method

Among the main assumptions on which is based the model proposed in this paper, the way of computing the stiffness of the members is likely to affect quantitatively the computed macroscale elasticity constants. Indeed, the boundary conditions used in the quasi-static finite element model presented in section 3.2 assume that the stiffness of a member is not affected by the potential presence of other members linked to the same particle. Moreover, the kinematic coupling constraints are likely to slightly overestimate the stiffness of the alumina particles.

The purpose of this section is to assess the effect of these assumptions by comparing the natural frequencies of two simple arrangements of two and eight particles (respectively one and twelve bridges), extracted either with a connector-based model or with a continuum element model (respectively 11886 and 140192 C3D10 elements of typical mesh size 15 nm). In these models, the mass unit is fg. Figs. 10 and 11 display both particles arrangements. The first one (two particles) allows to highlight the effect of the coupling constraints only, whereas the second one (eight particles) also includes the effect of several bridges linked to a given particle. In order to avoid differences related to silica bridge inertia, only the mass of the alumina particles is considered in both models, either as a nodal mass with convenient rotary inertia in the connector-based models, or as a distributed mass in the continuum models.

Fig. 10a and 11 show the relative deviation between the squared natural frequencies of similar eigenmodes extracted with both connector-based and continuum models using the latter as the baseline, respectively for two- and eight-particle models.

It is noteworthy that the strain energy of each mode of the connector-based two-particle model is mainly related to a single degree of freedom (cf. appendix C). In Fig. 10, the modes are labelled according to the connector degree of freedom that embeds the strain energy of the mode.

Fig. 10 shows that the torsion and bending modes of the two-particle models are barely affected by the coupling constraints used to compute the connector stiffness. Indeed, the relative deviation between both connector-based and continuum models is lower than 0.8%. This value rises up to 2.4% for the tension mode and reaches almost 8% for the shear mode. As discussed more in-depth in appendix C, these relative deviations are mainly related to the coupling constraints used to compute the connector stiffnesses, which constraint the distortion of the alumina particles (cf. Fig. 10b) and thus affect the strain energy.

In the case of the eight-particle arrangement, Fig. 11, the strain energy is distributed over several degrees of freedom of the connectors, depending on the mode. However, a global correlation may be drawn between the amount of strain energy embedded in a given degree of freedom and the relative deviation of the mode, through the relative deviation reported in Fig. 10 for each degree of freedom.

Indeed, the modes labelled 2 and 4 in Fig. 11 only involve the shear degree of freedom, and they correspond to the higher relative deviations, respectively 7 and 9%, i.e. the same order of magnitude than the relative deviation for the shear mode in Fig. 10 (8%). The lower relative deviation, 0.5%, is reported in Fig. 11 for the mode labelled 3, which only involves bending strain energy, which is the degree of freedom that exhibits the lower relative deviation in Fig. 10, less than 0.2%. The mode 9 only involves tension degree of freedom and its relative deviation is about 1%, whereas the relative deviation for the tension mode in Fig. 10 is 2.4%. The strain energy is shared between two degrees of freedom for modes 1 (85% shear and 15% torsion) and 5 (70% bending and 30% torsion) and their relative deviations are respectively of 5% and 1%, to be compared to the relative deviation of shear, torsion and bending modes in Fig. 10, respectively 8, 0.8 and
<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison of the natural vibration modes of a two-particle arrangement extracted with the connector-based model and the continuum model (the latter is the reference): relative deviation between the squared natural frequencies (a.) and cross-section view to highlight the deformed shape of the alumina particles (b.). The magnitude of the displacement fields computed with the continuum model are plotted. a = 103 nm, β = 45°, the other parameters are given in section 3.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Relative deviation between the squared natural frequencies of an eight-particle arrangement extracted with the connector-based model and the continuum model for similar vibration modes (the latter is the reference). The displacement fields computed with the continuum model are plotted. a = 103 nm, β = 45°, the other parameters are given in section 3.</figcaption></figure>

0.2%. The modes 6, 7 and 8 involve mainly the tension degree of freedom, which corresponds respectively to 72, 76 and 76% of the strain energy, and their relative deviations are respectively about 3, 2 and 4% (2.4% for the tensile mode in Fig. 10).

As a result, the relative deviations reported in Fig. 10 explain the major part of those reported in Fig. 11, which may indicate that the presence of three bridges connected to each particle in the eight-particle model has a weak influence on the results. Most of the deviation between both models is related to the kinematic coupling constraints used to compute the connector stiffness.

A tentative extrapolation of this analysis to the cylindrical volume element used in section 4.2 (cf. Fig. 8), would provide the relative deviation expected between the squared natural frequencies (and thus roughly the elasticity moduli, cf. Eqs. (8) and (9)) computed with the connector-based model and a hypothetical continuum model, at least the part related to the shear mode. All the modes reported in Fig. 8 exhibit a similar distribution of strain energy over the four degrees of freedom of the connectors, for instance 58–62% for tension, 34–38% for shear, 3% for bending and 1% for torsion. It would result in a relative

deviation typically about 5%. Such a level of deviation is not expected to affect the tendencies reported in the next section.

### 4.5. Homogenized properties versus material composition

The purpose of the truss-like model proposed in this paper is to provide homogenized elasticity constants as a function of material parameters, namely porosity ratio, p, and silica to alumina volume ratio, $V_{SiO_2}/V_{Al_2O_3}$. In this section, the porosity ratio ranges between 0.33 and 0.51. The method employed to build the truss-like model limits the alumina content to the initial value available in the unscaled random packing, $V_{Al_2O_3}/V_0$ (cf. section 3), which implies a lower bound to the silica to alumina volume ratio, for a given porosity ratio, Eq. (13), which is equivalent to a scaling factor $\alpha \ge 1$. The elasticity constants are computed for each porosity ratio up to $V_{SiO_2}/V_{Al_2O_3} = 2$. Such a high silica content is not likely to be suitable for the matrix of a ceramic matrix composite. A cylindrical volume element of 11 µm–15 µm in length after scaling, with $L/D = 3$, is considered. The number of particles is 10362 in these volume elements, and all of them are considered as representative of the material at macroscale according to the results of section 4.1.

$$ \frac{V_{SiO_2}}{V_{Al_2O_3}} \ge \frac{1-p}{V_{Al_2O_3}/V_0} - 1 \qquad \text{Eq. (13)} $$

As a whole, the natural frequencies of 77 representative volume elements have been extracted. The total computation time was 17 h (laptop, CPU @1.6 GHz, 16 GB memory). The resulting macroscale elasticity constants are reported in Fig. 12 as a function of porosity ratio and silica to alumina volume ratio.

Fig. 12 shows that the elasticity constants follow two distinct trends when increasing silica to alumina volume ratio. Up to $V_{SiO_2}/V_{Al_2O_3} \approx 0.5$, both Young’s and shear moduli increase while the Poisson’s ratio drops from 0.45 to 0.2. Then, above $V_{SiO_2}/V_{Al_2O_3} \approx 0.5$, the moduli slowly decrease and the Poisson’s ratio increases almost linearly. This behaviour seems to apply to all considered porosity ratios, even if it is concealed by the lower limit of silica to alumina ratio, Eq. (13), for $p \le 0.39$. Indeed, at a given porosity ratio, while increasing the silica to alumina volume ratio, the amount of silica asymptotically tends toward $1-p$. Concomitantly, the scale factor continuously increases and tends toward infinity, which means that the particles are increasingly spaced out. It results, in a first step, in an increase in the number of bridges thanks to a higher silica content, but the faster increase of inter-particle distances finally tends to decrease the number of bridges, cf. Fig. 13a. The mean nodal connectivity turns out to be a key characteristic that drives the
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Homogenized elasticity constants as a function of the material parameters: porosity ratio, *p*, and alumina to silica volume ratio. The grey filled areas correspond to alumina particles intersecting themselves ($\alpha < 1$) and are excluded.</figcaption></figure>

macroscale elasticity constants of the material.

The mean half-filling angle $\beta$ is plotted, in Fig. 13b, as a function of porosity ratio and silica to alumina volume ratio. It shows that this parameter is almost independent on porosity ratio. This is a direct consequence of the geometrical assumptions used in section 3.1, which lead to a high sensitivity of the inter-particle distance *a* to the angle $\beta$. However, the half-filling angle is widely spread around its mean value,

as shown in Fig. 13b. Indeed, the standard deviation ranges between 3.8° and 7.3° for the whole data reported in this graph.

As expected, Fig. 13c displays an increasing trend of the elasticity moduli as a function of the mean nodal connectivity, which seems to mainly drive the macroscale stiffness. The role of mean nodal connectivity on Poisson's ratio is of major influence. Indeed, Fig. 13d shows that the porosity ratio has a second order influence compared to connectivity.

These results show that a large range of porosity ratio and Young's modulus may be achieved by this family of porous ceramics, but only high porosity materials are suitable to reach a high value of Poisson's ratio. The model provides a useful insight into ways of designing porous materials with desired properties for a specific purpose.

## 5. Conclusion

This paper proposes to simulate the elastic behaviour of an open-cell porous ceramics, obtained by partial viscous sintering of submicronic alumina and nanometric silica particles, by a stochastic truss structure in which each silica bridge is modelled with a two-node connector element. Such a model is able to handle several thousands of alumina particles at a limited computational cost. This method is suitable for quasi-static analyses but also for transient dynamic procedures or natural frequencies extraction. The latter was used in the paper to assess the macroscale elasticity constants of the homogenized material, without requiring to apply any boundary conditions. The underlying assumptions related to inertia distribution in the material have been proved to have a minor influence on the results.

The behaviour of the members of the truss are computed, assuming an idealized geometry inspired by the hourglass-shape of the silica bridges and the theory of liquid menisci, assumed to drive the viscous sintering process.

The model is able to compute macroscale elasticity constants as a function of material parameters in a large range of porosity ratio and silica to alumina volume ratio, respectively between 33% and 51% and between 0.225 and 2. A large set of reasonable values of Young's modulus and Poisson's ratio were computed, spreading between 10 and 45 GPa for the former and between 0.15 and 0.45 for the latter. The elasticity moduli exhibit a maximum for a silica to alumina volume ratio of 0.5, irrespective of the porosity ratio. These results open the door to a more insightful tailored choice of the material parameters to match with the requirements of a given application.

Moreover, the limitations of the underlying assumptions have been clearly identified. The most restrictive one is due to the assumptions made on the boundary conditions, involving coupling constraints, in the computation of the connector stiffnesses, which lead to a poor description of the shear mode (cf. appendix C). This issue should be solved by calibrating the connector stiffness based on two-particle natural frequencies. This is the major route to improvement, together with a more accurate modelling of the silica mass distribution, especially in dynamic analyses. This latter issue is likely to be fixed by adding nodes to accurately support the silica mass, but the computational cost would increase consequently.

Before considering the fracture of the bridges and an application to composite issues, further works include a comparison with experimental data to improve the reliability of the model to the actual material. Among the numerous possibilities, the model may be improved to deal with higher porosity ratios, and a variability in the shape and size of alumina particles may be introduced. Such a variability would require to increase the number of parameters that describe the shape of the bridges, leading to a higher computational cost for member stiffness calculation only.

### CRediT authorship contribution statement

Thierry Canet: Writing – original draft, Software, Methodology,
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Mean nodal connectivity (a.) and mean half-filling angle $\beta$ (b.) as a function of the material parameters: porosity ratio, $p$, and alumina to silica volume ratio (the grey filled area marks the maximum standard deviation above and below the mean value of $\beta$); Young's (E) and shear (G) moduli (c.) and Poisson's ratio (d.) as a function of mean nodal connectivity for each porosity ratio.</figcaption></figure>

**Investigation, Data curation, Conceptualization.** **Gilles Dusserre:**
Writing – review & editing, Writing – original draft, Validation, Methodology, Funding acquisition, Data curation, Conceptualization.
**Thierry Cutard:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Funding acquisition.

## Funding sources

This research received a grant from the french Région Occitanie (N°00090338 / 21012975).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This research received a grant from the french Région Occitanie.

## Appendices.

### A. Sensitivity of single bridges stiffness calculation to elements size.

Four bridge geometries have been selected to study the sensitivity of the calculation of single bridge stiffness to elements size. The three first ones correspond to extreme values of separation distance and half-filling angle (references a., b. and c. in Fig. A.1). The last geometry corresponds to the average point of all the bridges of the cylindrical volume element described in section 4.2 (reference d. in Fig. A.1). For each of these four bridges, the stiffness under tension, torsion, bending and shear loadings have been calculated, as reported in section 3.2, with several nominal size of the C3D10 elements ranging between 7.5 nm and 25 nm (the mesh has been refined up to 5 nm for the smallest bridge, reference a. in Fig. A.1). Fig. A.2 reports the relative deviation from the results obtained with the finer mesh for each of the other meshes.

For a nominal element size of 10 nm, the relative deviation compared to a mesh size of 7.5 nm is below $1.6 \cdot 10^{-4}$ in all cases (Fig. A.2.b, Fig. A.2.c and Fig. A.2.d), except for the smallest bridge (up to 1.5% compared to a mesh size of 5 nm in bending, Fig. A.2.a, and below 1% for the other loading conditions). However, such short and thin bridges are expected to be exceptional in the volume elements (cf. Fig. A.1). As a result, the stiffness used in this paper, computed with an element size of 10 nm, can be considered to be accurate regarding the sensitivity to element size.

### B. Sensitivity of single bridges stiffness calculation to element type.

The influence of the element type on the member stiffness has also been studied based on the average bridge geometry ($a = 41$ nm, $\beta = 47.43°$, reference d. in Fig. A.1). C3D10 (tetrahedral, quadratic interpolation), C3D8 (hexahedral, linear interpolation) and C3D20 (hexahedral, quadratic interpolation) elements with nominal size ranging from 7.5 nm to 25 nm have been used. Fig B.1 shows that both quadratic elements C3D10 and C3D20 provide converged stiffness values, irrespective of the loading mode, within $\pm 0.055$ % around the stiffness computed with C3D20 elements of size 7.5 nm, assumed to be the most accurate. The results obtained with C3D10 elements of size 10 nm exhibit a relative deviation from this baseline below $10^{-4}$, which allows to conclude that the results reported in section 3.2 are accurate regarding the sensitivity to element type.

### C. Single connector natural frequencies.
Eq. (C.1-C.4) provide the analytical expressions of the natural angular frequencies, $\omega_{dof} = 2\pi f_{dof}$, of a two-mass system (mass $m$ and moment of inertia $I$) linked by a linear elastic connector of length $l$ of independent stiffness $K_{dof}$ for each degree of freedom. On the one hand, Eq. (C.1-C.3) are straightforward. Indeed, assuming a harmonic motion, the Newton's second law of motion, applied to the isolated system, implies that each of the three first modes (namely torsion, Eq. (C.1), bending, Eq. (C.2) and tension, Eq. (C.3)) involves out-of-phase displacements of a single degree of freedom. On the other hand, the so-called shear mode, Eq. (C.4), involves coupled shear and bending displacements. Indeed, the conservation of angular momentum requires that $U_{R3} = \frac{m}{I} \frac{L}{2} U_2$ and the Newton's second law of motion leads to an out-of-phase shear displacement and an in-phase bending rotation. As a result, the bending moment vanishes, together with the relative rotation. Therefore, Eq. (C.4) only involves shear stiffness, $K_2$, but accounts for linear and angular momenta. Moreover, the displacement, used to compute the force transmitted by the CONN3D2 elements, considers the relative displacement in the element local coordinate system. The relative shear displacement, work conjugate to the shear force, is therefore a function of both translations and rotations of the nodes. Considering infinitesimal displacements, all these considerations lead to Eq. (C.4).

$$ \omega_{R1}^2 = \frac{2K_{R1}}{I} \qquad \text{Eq. (C.1)} $$

$$ \omega_{R2}^2 = \omega_{R3}^2 = \frac{2K_{R2}}{I} \qquad \text{Eq. (C.2)} $$

$$ \omega_1^2 = \frac{2K_1}{m} \qquad \text{Eq. (C.3)} $$

$$ \omega_2^2 = \omega_3^2 = \frac{2K_2}{m} \left( 1 + \frac{m}{I} \left( \frac{L}{2} \right)^2 \right) \qquad \text{Eq. (C.4)} $$

The relative deviations between the squared natural frequencies extracted using a single two-node element or using a continuum model (cf. section 4.4 and Fig. 10) are respectively of – 0.742%, – 0.153%, – 1.97% and – 7.70% in torsion, bending, tension and shear (Table C.1). The relative deviations are below $4 \cdot 10^{-5}$ between connector-based and analytical models.

Moreover, the natural frequencies of another continuum model have been extracted (cf. Table C.1), assuming the same configuration used in section 3.2 to calibrate the connector stiffnesses (cf. Fig. 3). As previously, the silica mass is not considered. The alumina mass and rotary inertia are assumed to be punctual at the reference node of the rigid surface of both half-particles. The relative deviations are below 0.5% with the connector-based model. This comparison shows that the deviations between the connector-based model and the continuum model (Fig. 10) are mainly related to the way the stiffnesses are calibrated from quasi-static analyses with simplified boundary conditions. This configuration affects the displacement field (Fig. 10b) and tends to underestimate the strain energy.

Another calibration strategy based on Eq. (C.1-C.4) would consist in calibrating the connector stiffnesses using the natural frequencies of the two-particle continuum model used in section 4.4. This would reduce the maximum relative deviations reported in Fig. 11 for an eight-particle system from 9.3% to 2.1%. This strategy was not investigated yet for a large number of particles.

<table><caption>**Table C.1**

Natural frequencies (GHz) of a two-particle system extracted using several models.</caption><thead><tr><th></th><th>f<sub>R1</sub></th><th>f<sub>R2</sub></th><th>f<sub>1</sub></th><th>f<sub>2</sub></th></tr></thead><tbody><tr><td>Continuum model (Fig. 10)</td><td>1.43325</td><td>1.5463</td><td>3.7122</td><td>4.6704</td></tr><tr><td>Single connector model</td><td>1.42792</td><td>1.5451</td><td>3.6755</td><td>4.4870</td></tr><tr><td>Analytical model (Eq. (C.1-C.4))</td><td>1.427924</td><td>1.54511</td><td>3.67558</td><td>4.48689</td></tr><tr><td>Continuum model with punctual mass (Fig. 3)</td><td>1.42799</td><td>1.54514</td><td>3.6673</td><td>4.4851</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. A.1.** Separation distance and half-filling angle of each bridge of a cylindrical volume element of 13445 particles (black dots, cf. section 4.2 for a description of the volume element), and particular bridges used for sensitivity study to elements size (red dots: $a = 10$ nm, $\beta = 11.22^{\circ}$ (a.), $a = 10$ nm, $\beta = 67.89^{\circ}$ (b.), $a = 41$ nm, $\beta = 47.43^{\circ}$ (c.) and $a = 237$ nm, $\beta = 74.8^{\circ}$ (d.).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. A.2.** Relative deviation of the stiffnesses in tension, torsion, bending and shear computed with different element sizes, compared to the finer mesh (a.: 5 nm; b., c. and d.: 7.5 nm), for four different bridge geometries: *a* = 10 nm, *β* = 11.22° (a.), *a* = 10 nm, *β* = 67.89° (b.), *a* = 237 nm, *β* = 74.8° (c.) and *a* = 41 nm, *β* = 47.43° (d.).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. B.1.** Stiffness of an average bridge (a = 41 nm, β = 47.43°, reference d. in Fig. A.1) computed using C3D10 (tetrahedral, quadratic interpolation), C3D8 (hexahedral, linear interpolation) and C3D20 (hexahedral, quadratic interpolation) with different nominal element sizes (a. Tension, b. shear, c. Torsion, d. bending). The grey filled areas mark a range of ±0.1% around the a priori most accurate value (C3D20 elements of size 7.5 nm).</figcaption></figure>

Data availability

Data will be made available on request.
## References

Aghajani, S., Wu, C., Li, Q., Fang, J., 2024. Additively manufactured composite lattices: a state-of-the-art review on fabrications, architectures, constituent materials, mechanical properties, and future directions. Thin-Walled Struct. 197, 111539. https://doi.org/10.1016/j.tws.2023.111539.

Ahmad, R., Ha, J.-H., Song, I.-H., 2014. Processing methods for the Preparation of porous ceramics. J. Korean Powder Metall. Inst. 21, 389-398. https://doi.org/10.4150/KPMI.2014.21.5.389.

Asadi, F., André, D., Emam, S., Doumalin, P., Huger, M., 2023. Advances in micro-mechanical modeling using a bonded-particle material and periodic homogenization within discrete element framework applied to heterogeneous ceramics. J. Eur. Ceram. Soc. 43, 1205-1219. https://doi.org/10.1016/j.jeurceramsoc.2022.10.073.

Asahina, D., Aoyagi, K., Kim, K., Birkholzer, J.T., Bolander, J.E., 2017. Elastically-homogeneous lattice models of damage in geomaterials. Comput. Geotech. 81, 195-206. https://doi.org/10.1016/j.compgeo.2016.08.015.

Ashby, M.F., 2006. The properties of foams and lattices. Phil. Trans. R. Soc. A. 364, 15-30. https://doi.org/10.1098/rsta.2005.1678.

ASTM C215, 2020. Test method for fundamental transverse, longitudinal, and torsional resonant frequencies of concrete specimens. https://doi.org/10.1520/C0215-19.

ASTM E1875, 2021. Test method for dynamic youngs modulus, shear modulus, and Poisson's ratio by. Sonic Resonance. https://doi.org/10.1520/E1875-20A.

ASTM E1876, 2022. Test method for dynamic youngs modulus, shear modulus, and. Poisson's Ratio by Impulse Excitation of Vibration. https://doi.org/10.1520/E1876-22.

Belardi, V.G., Trupiano, S., Fanelli, P., Vivio, F., 2024. Overall elastic characterization of equivalent FE models for aluminum foams through computational homogenization approach and genetic algorithm optimization. Eur. J. Mech. Solid. 103, 105189. https://doi.org/10.1016/j.euromechsol.2023.105189.

Ben Ramdane, C., Julian-Jankowiak, A., Valle, R., Renollet, Y., Parlier, M., Martin, E., Diss, P., 2017. Microstructure and mechanical behaviour of a Nextel™610/alumina weak matrix composite subjected to tensile and compressive loadings. J. Eur. Ceram. Soc. 37, 2919-2932. https://doi.org/10.1016/j.jeurceramsoc.2017.02.042.

Bhate, D., Penick, C.A., Ferry, L.A., Lee, C., 2019. Classification and selection of cellular materials in mechanical design: engineering and biomimetic approaches. Design 3, 19. https://doi.org/10.3390/designs3010019.

Boutenel, F., 2020. Comportement d'une matrice alumine-silice au cours de l'élaboration de composites à matrice céramique oxyde/oxyde. Ecole des Mines d'Albi-Carmes.

Boutenel, F., Dusserre, G., Aimable, A., Chartier, T., Cutard, T., 2021. Rheophysical study of dispersed alumina suspensions. Powder Technol. 393, 630-638. https://doi.org/10.1016/j.powtec.2021.08.016.

Chen, Y., Wang, N., Ola, O., Xia, Y., Zhu, Y., 2021. Porous ceramics: light in weight but heavy in energy and environment technologies. Mater. Sci. Eng. R Rep. 143, 100589. https://doi.org/10.1016/j.mser.2020.100589.

Cheng, L., Bai, J., To, A.C., 2019. Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput. Methods Appl. Mech. Eng. 344, 334-359. https://doi.org/10.1016/j.cma.2018.10.010.

Comeforo, J.E., Hursh, R.K., 1952. Wetting of Al2O3-SiO2 refractories by molten glass: I, measurement of wetting. J. Am. Ceram. Soc. 35, 130-134. https://doi.org/10.1111/j.1151-2916.1952.tb13085.x.

De Faoite, D., Browne, D.J., Chang-Díaz, F.R., Stanton, K.T., 2012. A review of the processing, composition, and temperature-dependent mechanical and thermal properties of dielectric technical ceramics. J. Mater. Sci. 47, 4211-4235. https://doi.org/10.1007/s10853-011-6140-1.

Deschamps, T., Margueritat, J., Martinet, C., Mermet, A., Champagnon, B., 2014. Elastic moduli of permanently densified silica glasses. Sci. Rep. 4, 7193. https://doi.org/10.1038/srep07193.

Donaldson, F.E., Pankaj, P., Law, A.H., Simpson, A.H., 2008. Virtual trabecular bone models and their mechanical response. Proc. Inst. Mech. Eng. H 222, 1185-1195. https://doi.org/10.1243/09544119JEM408.

Falco, S., De Cola, F., Petrinic, N., 2017. A method for the generation of 3D representative models of granular based materials. Int. J. Numer. Methods Eng. 112, 338-359. https://doi.org/10.1002/nme.5521.

Fisher, R.A., 1926. On the capillary forces in an ideal soil; correction of formulae given by W. B. Haines. J. Agric. Sci. 16, 492-505. https://doi.org/10.1017/S0021859600007838.

Garnier, D.-H., Schmidt, M.-P., Rohmer, D., 2022. Growth of oriented orthotropic structures with reaction/diffusion. Struct. Multidiscip. Optim. 65, 327. https://doi.org/10.1007/s00158-022-03395-7.

Guel, N., 2018. Comportement mécanique de composites oxydes: Relations procédé-microstructure-propriétés. Université de Lyon.

Guidault, P.-A., Marie-France Soulé de Lafont, Pierre-Alain Boucard, 2015. Modélisation et identification d'un connecteur élément fini non-linéaire pour la simulation d'assemblages boulonnés. 12e Colloque national en calcul des structures.

Ibrahimbegovic, A., Mejia-Nava, R.-A., 2023. Truss model: general theorems and methods of force, displacement and finite elements. In: Ibrahimbegovic, A., Mejia-Nava, R.-A. (Eds.), Structural Engineering: Models and Methods for Statics, Instability and Inelasticity. Springer International Publishing, Cham, pp. 15-58. https://doi.org/10.1007/978-3-031-23592-4_2.

Imediegwu, C., Clarke, D., Carter, F., Grimm, U., Jowers, I., Moat, R., 2023. Mechanical characterisation of novel aperiodic lattice structures. Mater. Des. 229, 111922. https://doi.org/10.1016/j.matdes.2023.111922.

Khoda, B., Ahsan, A.M.M.N., Shovon, A.N., Alam, A.I., 2021. 3D metal lattice structure manufacturing with continuous rods. Sci. Rep. 11, 434. https://doi.org/10.1038/s41598-020-79826-6.

Koch, D., Tushtev, K., Grathwohl, G., 2008. Ceramic composites: experimental analysis and modeling of mechanical properties. Compos. Sci. Technol. 68, 1165-1172. https://doi.org/10.1016/j.compscitech.2007.06.029.

Li, T., Li, L., Liu, J., Zhou, H., 2022. Influence of hydrate participation on the mechanical behaviour of fine-grained sediments under one-dimensional compression: a DEM study. Granul. Matter 24, 32. https://doi.org/10.1007/s10035-021-01195-2.

Lian, G., Thornton, C., Adam, J.M., 1993. A theoretical study of the liquid bridge forces between 2 rigid spheric bodies. J. Colloid Interface Sci. 161, 138-147.

Marchelli, G.L.S., 2011. Highly porous 3D printed stoneware ceramic. J. Ceram. Sci. Technol. https://doi.org/10.4416/JCST2011-00042.

Margarida, M., Lima, R.A., Monteiro, R. da C.C., 2008. Viscous sintering in a glass-alumina system. MSF 587-588, 143-147. https://doi.org/10.4028/www.scientific.net/MSF.587-588.143.

Marin, E., 2021. Comportement mécanique de composites oxyde/oxyde pour structures aéronautiques. Université de Bordeaux.

Maurizi, M., Edwards, B.W., Gao, C., Greer, J.R., Berto, F., 2022. Fracture resistance of 3D nano-architected lattice materials. Extre. Mech. Lett. 56, 101883. https://doi.org/10.1016/j.eml.2022.101883.

Maxwell, J.C., 1864. L. On the calculation of the equilibrium and stiffness of frames. London, Edinburgh Dublin Phil. Mag. J. Sci. 27, 294-299. https://doi.org/10.1080/1478644608463668.

Meza, L.R., Phlipot, G.P., Portela, C.M., Maggi, A., Montemayor, L.C., Comella, A., Kochmann, D.M., Greer, J.R., 2017. Reexamining the mechanical property space of three-dimensional lattice architectures. Acta Mater. 140, 424-432. https://doi.org/10.1016/j.actamat.2017.08.052.

Mohammadi, A., Sadighi, M., Hedayati, R., Ganghoffer, J.-F., 2024. Asymptotic homogenization of tesseract lattice structures. Eur. J. Mech. Solid. 107, 105368. https://doi.org/10.1016/j.euromechsol.2024.105368.

Pei, E., Kabir, I., Breski, T., Godec, D., Nordin, A., 2022. A review of geometric dimensioning and tolerancing (GD&T) of additive manufacturing and powder bed fusion lattices. Prog Addit Manuf 7, 1297-1305. https://doi.org/10.1007/s40964-022-00304-8.

Portela, C.M., Greer, J.R., Kochmann, D.M., 2018. Impact of node geometry on the effective stiffness of non-slender three-dimensional truss lattice architectures. Extre. Mech. Lett. 22, 138-148. https://doi.org/10.1016/j.eml.2018.06.004.

Qi, D., Yu, H., Liu, M., Huang, H., Xu, S., Xia, Y., Qian, G., Wu, W., 2019. Mechanical behaviors of SLM additive manufactured octet-truss and truncated-octahedron lattice structures with uniform and taper beams. Int. J. Mech. Solid. 163, 105091. https://doi.org/10.1016/j.ijmecsci.2019.105091.

Qi, G., Zeng, Y., Chen, J., 2022. Preparation of porous SnO2-based ceramics with lattice structure by DLP. Ceram. Int. 48, 14568-14577. https://doi.org/10.1016/j.ceramint.2022.02.0350.

Ramirez-Chavez, I.E., Anderson, D., Sharma, R., Lee, C., Bhate, D., 2022. A classification of aperiodic architected cellular materials. Design 6, 63. https://doi.org/10.3390/designs6040063.

Shaat, M., Wagih, A., 2020. Hinged-3D metamaterials with giant and strain-independent Poisson's ratios. Sci. Rep. 10, 2228. https://doi.org/10.1038/s41598-020-59205-x.

Thillaithevan, D., Bruce, P., Santer, M., 2021. Stress-constrained optimization using graded lattice microstructures. Struct. Multidiscip. Optim. 63, 721-740. https://doi.org/10.1007/s00158-020-02723-z.

Wang, C., Gu, X., Zhu, J., Zhou, H., Li, S., Zhang, W., 2020. Concurrent design of hierarchical structures with three-dimensional parameterized lattice microstructures for additive manufacturing. Struct. Multidiscip. Optim. 61, 869-894. https://doi.org/10.1007/s00158-019-02408-2.

Wang, C., Zhu, J.H., Zhang, W.H., Li, S.Y., Kong, J., 2018. Concurrent topology optimization design of structures and non-uniform parameterized lattice microstructures. Struct. Multidiscip. Optim. 58, 35-50. https://doi.org/10.1007/s00158-018-2009-0.

Wang, S., Ma, Y., Deng, Z., 2022. Stretching-dominated truss lattice materials: elastic anisotropy evaluation, control, and design. Compos. Struct. 298, 116004. https://doi.org/10.1016/j.compstruct.2022.116004.

Weigert, T., Ripperger, S., 1999. Calculation of the liquid bridge volume and bulk saturation from the half-filling angle. Part. Part. Syst. Char. 16, 200-244.

Xu, S., Shen, J., Zhou, S., Huang, X., Xie, Y.M., 2016. Design of lattice structures with controlled anisotropy. Mater. Des. 93, 443-447. https://doi.org/10.1016/j.matdes.2016.01.007.

Zok, Frank W., Levi, Carlos G., 2001. Mechanical properties of porous-matrix ceramic composites. Adv. Eng. Mater. 3, 15-23. https://doi.org/10.1002/1527-2648 (200101)3:1/2<15::AID-ADEM15>3.0.CO;2-A.