# Micromechanical explanation of elasticity and strength of gypsum: from elongated anisotropic crystals to isotropic porous polycrystals

Julien Sanahuja, Luc Dormieux, Sylvain Meille, Christian Hellmich, Andreas Fritsch

► To cite this version:

Julien Sanahuja, Luc Dormieux, Sylvain Meille, Christian Hellmich, Andreas Fritsch. Micromechanical explanation of elasticity and strength of gypsum: from elongated anisotropic crystals to isotropic porous polycrystals. Journal of Engineering Mechanics - ASCE, 2010, 136 (2), pp.239-253. 10.1061/(ASCE)EM.1943-7889.0000072. hal-00587082

HAL Id: hal-00587082

https://hal.science/hal-00587082v1

Submitted on 26 Feb 2021

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 CC BY 4.0 - Attribution - International License
# Micromechanical Explanation of Elasticity and Strength of Gypsum: From Elongated Anisotropic Crystals to Isotropic Porous Polycrystals

Julien Sanahuja, Ph.D.<fnref n="1" />; Luc Dormieux<fnref n="2" />; Sylvain Meille, Ph.D.<fnref n="3" />; Christian Hellmich<fnref n="4" />; and Andreas Fritsch, Ph.D.<fnref n="5" />

**Abstract:** Gypsum is made up of interlocked and elongated crystals. The random nature of its morphology suggests to resort to homogenization of random media to investigate its mechanical properties from the scale of the single crystals upwards. Unfortunately, the usual homogenization schemes fail to quantitatively predict the influence of the porosity on the effective Young's modulus of gypsum. This is clearly due to the inability of such approaches to take into account the elongated nature of the crystals. A modification of the classical self-consistent scheme is proposed. It is validated against elastic characteristics computed by finite element analyses, and also against experiments on real dried gypsum samples (with empty pores). Finally, a strength model based on brittle failure is presented. The whole strength domain in the space of macroscopic principal stresses is derived. The comparison to experimental data in both simple tension and simple compression is remarkably good.

**Keywords:** Homogenization; Gypsum; Microstructures; Elasticity; Strength.

## Introduction

As shown on the scanning electron microscope (SEM) picture on Fig. 1, gypsum is made up of interlocking elongated crystals with random orientation. This suggests that gypsum belongs to the family of polycrystal-type materials. In the framework of classical Eshelby-based homogenization tools for random media, such materials are typically modeled by the self-consistent scheme (Kröner 1977). Unfortunately, it is well known that a vanishing effective stiffness is obtained for a solid volume fraction of 0.5, while experimental data clearly reveal that gypsum exhibits nonzero stiffness properties at least for solid volume fractions as low as 0.3 (Sanahuja et al. 2008).

On the other hand, it is not surprising that the other classical schemes (differential, Mori Tanaka, and Christensen and Lo) fail to predict correctly the effective stiffness as a function of the solid

volume fraction. Indeed, all these schemes predict a nonvanishing stiffness as long as the solid volume fraction is nonzero. This is related to the fact that they all attempt to take into account the connectedness of a matrix phase which obviously does not exist in gypsum at the scale revealing the crystals.

In conclusion, keeping the idea of dealing with a polycrystal-line material, there is a need for extending the domain of applicability of the classical self-consistent scheme, so as to be able to predict a critical solid volume fraction much lower than 0.5. This in short is the topic of the first part of the paper, presenting a microelastic model for gypsum and its experimental validation. After having discussed the elastic properties of gypsum, we then address its strength in the framework of brittle failure. The purpose is to define the domain of admissible stress states. This analysis is validated by the comparison of theoretical and experimental results from uniaxial compression and tension tests.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** SEM picture of a gypsum sample</figcaption></figure>

<footnote n="1">¹Lafarge Centre de Recherche, 95 rue du Montmurier, BP15, 38291 St. Quentin Fallavier cedex, France; presently, EDF R&D, Route de Sens, Ecuelles, 77250 Moret sur Loing, France (corresponding author). E-mail: julien.sanahuja@edf.fr</footnote>

<footnote n="2">²Professor, École Nationale des Ponts et Chaussées, 6-8 Ave. Blaise Pascal, 77455 Champs-sur-Marne, France.</footnote>

<footnote n="3">³INSA Lyon, Mateis Laboratory, Université de Lyon, UMR CNRS 5510, 20 Ave. Albert Einstein, bâtiment Blaise Pascal, 69621 Villeurbanne, France.</footnote>

<footnote n="4">⁴Professor, Institut für Mechanik der Werkstoffe und Strukturen, Technische Universität Wien, Karlsplatz 13/202, A-1040 Wien, Austria.</footnote>

<footnote n="5">⁵Institut für Mechanik der Werkstoffe und Strukturen, Technische Universität Wien, Karlsplatz 13/202, A-1040 Wien, Austria.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Orientation of a given crystal characterized by the three Euler angles $\theta$, $\phi$, $\psi$; the crystal frame $(\underline{u}_1, \underline{u}_2, \underline{u}_3)$ is defined with respect to the reference frame $(\underline{\epsilon}_1, \underline{\epsilon}_2, \underline{\epsilon}_3)$ [the spherical frame $(\underline{\epsilon}_r, \underline{\epsilon}_\theta, \underline{\epsilon}_\phi)$ appears as dotted lines]</figcaption></figure>

## Modified Self-Consistent Scheme for Elongated Particle Phases

As pointed out in the introduction, classical homogenization schemes fail to properly predict the effective stiffness of gypsum as a function of porosity. In order to overcome this difficulty, a modified self-consistent scheme is developed in the present section. The aim is to take into account the elongated shape of the solid particles. This new scheme is validated against numerical computations (numerical homogenization performed on clearly defined microstructures) available in the literature.

## Theoretical Developments

We propose to modify the classical self-consistent scheme through replacing the spherical inclusion used to estimate the average strain in the solid phase by a set of infinitely many elongated spheroids (prolates) representing single crystals characterized by anisotropic elasticity. All these spheroids are identical with respect to shape and material behavior, while being isotropically oriented in all directions in space, i.e., the spheroid's Euler angles $\theta$, $\phi$, and $\psi$ (see Fig. 2) being uniformly distributed between zero and $\pi$, $2\pi$, and $2\pi$, respectively. This has two implications: (1) the stiffness tensor related to a single crystal is a function of the Euler angles $C_s(\theta, \phi, \psi)$, while its components are orientation-independent in the frame $(\underline{u}_1, \underline{u}_2, \underline{u}_3)$ (see Fig. 2), and (2) the (overall) effective stiffness tensor of polycrystal “gypsum” is isotropic; it is denoted by $C^{\text{SCS}}$.

This micromechanical representation implies that two morphological parameters characterize the microstructure: (1) the porosity and (2) the aspect ratio of the solid spheroids, denoted by $\varphi$ and $r_s$, respectively.

The average strain in the pore space $\langle \boldsymbol{\epsilon} \rangle_p$ is classically estimated by the uniform strain arising in a spherical cavity embedded into an infinite medium whose stiffness is $C^{\text{SCS}}$, with uniform strain (tensor $E_0$) boundary conditions at infinity (see left part of Fig. 3)

$$ \langle \boldsymbol{\epsilon} \rangle_p = (I - S_{\text{sph}}^{\text{SCS}})^{-1} : E_0 \quad (1) $$

with $S_{\text{sph}}^{\text{SCS}}$ the Eshelby tensor (Eshelby 1957) of a sphere in the medium of stiffness $C^{\text{SCS}}$.

The average strain in the set of crystals whose orientation is parametrized by $(\theta, \phi, \psi)$ is estimated by the uniform strain arising in a prolate spheroid whose axis of revolution is $\underline{\epsilon}_r$ (characterized by $\theta, \phi$), made up of the same material as these crystals, embedded in the infinite medium of stiffness $C^{\text{SCS}}$, with uniform strain (tensor $E_0$) boundary conditions at infinity (see right part of Fig. 3); being a specific case of the Eshelby inhomogeneity problem (Eshelby 1957)

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Auxiliary problems of elasticity to solve to implement the self-consistent scheme proposed (sphere to represent the pore space $p$ and spheroids to represent the solid $s$)</figcaption></figure>

$$ \langle \boldsymbol{\epsilon} \rangle_{\theta, \phi, \psi} = \{I + P_{\text{pro}}^{\text{SCS}}(\theta, \phi): [C_s(\theta, \phi, \psi) - C^{\text{SCS}}]\}^{-1} : E_0 \quad (2) $$

with $P_{\text{pro}}^{\text{SCS}}(\theta, \phi)$ the Hill tensor of a prolate spheroid, whose axis of revolution is characterized by $\theta, \phi$, in the medium of stiffness $C^{\text{SCS}}$.

The average strain in the whole solid domain is then obtained from an angular average, since the orientation distribution of the crystals was introduced as isotropic

$$ \langle \boldsymbol{\epsilon} \rangle_s = \int_{\psi=0}^{2\pi} \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} \langle \boldsymbol{\epsilon} \rangle_{\theta, \phi, \psi} \frac{\sin \theta}{8\pi^2} d\theta d\phi d\psi \quad (3) $$

A similar expression holds for the stresses averaged over the entire solid domain

$$ \langle \boldsymbol{\sigma} \rangle_s = \int_{\psi=0}^{2\pi} \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} C_s(\theta, \phi, \psi) \langle \boldsymbol{\epsilon} \rangle_{\theta, \phi, \psi} \frac{\sin \theta}{8\pi^2} d\theta d\phi d\psi \quad (4) $$

Finally, the strain tensor $E$ applied on the boundary of the representative element volume (REV) and the average stress $\Sigma$ over the REV then read

$$ E = \varphi \langle \boldsymbol{\epsilon} \rangle_p + (1 - \varphi) \langle \boldsymbol{\epsilon} \rangle_s \quad (5a) $$

and

$$ \Sigma = (1 - \varphi) \langle \boldsymbol{\sigma} \rangle_s \quad (5b) $$

Eq. (5a) allows for relating the REV-related macroscopic strains $E$ to the (auxiliary) strains $E_0$ imposed to the infinite matrices of Fig. 3. More specifically, insertion of Eqs. (1)–(3) into Eq. (5a) yields a relation between $E_0$ and $E$. Use of this relation in Eqs. (2), (4), and (5b) yields a relation between the macroscopic stresses $\Sigma$ and the macroscopic strains $E$, expressed through the (homogenized and effective) self-consistent stiffness tensor of the polycrystal

$$ \begin{aligned} C^{\text{SCS}} = (1 - \varphi) & \int_{\psi=0}^{2\pi} \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} C_s(\theta, \phi, \psi): \{I + P_{\text{pro}}^{\text{SCS}}(\theta, \phi): \\ & [C_s(\theta, \phi, \psi) - C^{\text{SCS}}]\}^{-1} \frac{\sin \theta}{8\pi^2} d\theta d\phi d\psi: \\ & \left[ (1 - \varphi) \int_{\psi=0}^{2\pi} \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} \{I + P_{\text{pro}}^{\text{SCS}}(\theta, \phi): \right. \\ & \left. [C_s(\theta, \phi, \psi) - C^{\text{SCS}}]\}^{-1} \frac{\sin \theta}{8\pi^2} d\theta d\phi d\psi + \varphi(I - S_{\text{sph}}^{\text{SCS}})^{-1} \right]^{-1} \end{aligned} \quad (6) $$
The self-consistent effective stiffness being isotropic, it is written as $C^{\text{SCS}}=3k^{\text{SCS}}J+2g^{\text{SCS}}K$. From Eq. (6), two nonlinear scalar equations can be extracted, by projection of $C^{\text{SCS}}$ onto $J$ and $K$. The self-consistent estimates $k^{\text{SCS}}$ and $g^{\text{SCS}}$ of the effective bulk and shear moduli are the positive solutions to these equations. Unfortunately, the expressions of these equations are rather lengthy, except for the classical case of the sphere ($r_s \to 1$) and for the limit case of the needle ($r_s \to \infty$) [see Eqs. (7) and (8)].

## Numerical Validation

The proposed self-consistent scheme considers the material's morphology in an implicit, rather than a direct way. This concerns two aspects: (1) rather than modeling the porous space between the elongated needles in a geometrically exact way, the isotropic distribution of different pore geometries between the crystals is simply represented by a spherical pore geometry and (2) all crystals and pores are set into a matrix made up of the homogenized material, instead of considering the precise load transfer mechanisms between individual crystals. This situation calls for careful validation of the proposed scheme, first with respect to computations which do consider both the precise pore geometries between crystal particles and the load carrying mechanisms between individual particles. Fortunately, such FEM homogenization simulations have been performed by Meille at NIST (Meille and Garboczi 2001), on microstructures which were digitally generated, by means of introduction of elongated solid particles into an initially empty cube. Simulations performed on these clearly defined microstructures allow for checking of whether the proposed self-consistent scheme is appropriate to model porous polycrystals made up of elongated solid particles or not.

In (Meille and Garboczi 2001), solid bars (rectangular parallelepipeds) with aspect ratio $7:1:1$ were put into a cube (unit cell) with a size 5 to 10 times larger than the longest dimension of the bars. Three ways of placing the solid bars were realized:

* All randomly placed bars are parallel to one of the axes of the cube, and they may overlap (in this case, all bars consist of $21 \times 3 \times 3$ voxels);

* All randomly placed bars are parallel to one of the axes of the cube, but their hard cores (measuring $19 \times 1 \times 1$ voxels) are not allowed to overlap (also in this case, all bars consist of $21 \times 3 \times 3$ voxels); and

* The bars are randomly oriented in space and they may overlap.

Thus, the first two configurations of bars (elongated crystals) lead to cubically symmetric material behavior, while the third configuration implies isotropic material behavior. The three configurations (morphologies) were realized in various FEM discretizations, as to avoid mesh dependencies. Interestingly, the (mesh-independent) results of all three configurations yield virtually the same relationship between the effective Young's modulus $E^{\text{eff}}$ and the porosity $\phi$. Thus, data computed from the three types of microstructural configurations are represented by one and the same cross symbol in Fig. 4. More specifically, the FEM computations presented on the left part of Fig. 4 were performed with solid particles exhibiting Young's modulus $E_s=46$ GPa and Poisson's ratio $\nu_s=0.33$. These values correspond to the isotropic stiffness tensor which is obtained by angular averaging of the anisotropic gypsum crystal stiffness tensor (see Haussühl 1965 for its components).

The left part of Fig. 4 compares the FEM computations to estimates from the proposed modified self-consistent scheme. The latter have been determined by using prolate spheroids with three different values for the aspect ratio, namely, 7, 10, and 15.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Effective Young's modulus and Poisson's ratio estimated by the self-consistent scheme with spherical pores and elongated solid particles (with different aspect ratios $r_s$), compared to values computed by FEM analyses (Meille and Garboczi 2001); the classical schemes indicated by dotted lines are (when ordered by increasing stiffness): self-consistent scheme with spherical solid inclusions, differential scheme, Christensen and Lo scheme, and Mori-Tanaka scheme</figcaption></figure>

Note that a critical porosity appears, above which the effective Young's modulus vanishes. This critical porosity is denoted by $\phi^c$. It is found to depend on the aspect ratio $r_s$ of the prolates.

An excellent agreement between the estimates and the computations is observed. At lower porosities ($\phi \le 0.4$), and in the range (7,15) of aspect ratios explored here, the self-consistent estimates hardly depend on the aspect ratio. At higher porosities, an aspect ratio of about 15 (to be compared to the aspect ratio 7 of the parallelepipeds) yields a very good agreement between the self-consistent estimates and the FEM computations.

The agreement between self-consistent estimations and FEM computations is confirmed by the prediction of the effective Poisson's ratio $\nu^{\text{eff}}$. The latter is plotted on the right part of Fig. 4, as function of the porosity, for three different values of the Poisson's ratio $\nu_s$ of the solid spheroids, namely, 0, 0.2 and 0.33 (see $\nu^{\text{eff}}$ at $\phi=0$). These different values have been chosen in (Meille and Garboczi 2001), in order to investigate the influence of $\nu_s$ on the curve $\nu^{\text{eff}}(\phi)$. We here confirm a fact already found by means of FEM computations: the effective Poisson's ratio tends to a unique value (independent of the Poisson's ratio of the solid) when the porosity reaches the critical porosity. This fact has also been found in plane stress elasticity from two-dimensional numerical computations (Day et al. 1992) and from a theoretical point of view (Cherkaev et al. 1992). On Fig. 4, we also plotted, in terms of dotted lines, the porosity-effective characteristics relationships defined through the most frequently used classical homogenization schemes. Going from lower to higher effective moduli, they are: self-consistent scheme with spherical shapes, differential scheme, Christensen and Lo scheme, and Mori Tanaka scheme. The elastic characteristics of the solid phase are still $E_s=46$ GPa and $\nu_s=0.33$. Clearly enough, these classical homogenization schemes are not able to properly predict the effective elastic characteristics of polycrystals made up of lengthened solid particles, as soon as $\phi \ge 0.3$.

This comparison of the effective elastic characteristics estimated from the self-consistent scheme to the ones computed by FEM on clearly defined digital microstructures allows us to highlight a family of morphologies whose effective stiffness can be reasonably estimated by means of the proposed self-consistent scheme. Thus, a self-consistent scheme with prolate spheroids provides remarkably good estimates of the effective elastic characteristics of a porous polycrystal made up of lengthened solid grains. The next subsection pushes forward the investigation, comparing the critical porosities $\phi^c$, already uncovered on the left part of Fig. 4, to results from geometrical computations.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Critical solid volume fraction as a function of the aspect ratio of the solid spheroids, estimated by the self-consistent scheme, and percolation threshold numerically computed (Garboczi et al. 1995) (logarithmic scale to reveal the evolutions for very low and very high aspect ratios)</figcaption></figure>

## Critical Porosity

The classical self-consistent scheme is based on spherical phase shapes of both the solid and the pore phase. This scheme implies a critical porosity of $\varphi^c=1/2$, above which the effective Young's modulus vanishes. However, when considering spheroidal solid phases, this critical porosity is found to depend on the aspect ratio $r_s$ of the solid inhomogeneities. While the dependence of $\varphi^c$ with respect to $r_s$ has already been observed for prolate spheroids in (Fritsch et al. 2006), we here show this dependence for oblates. In fact, $\varphi^c$ even exclusively depends on $r_s$; the critical porosity is independent of the elastic characteristics of the solid phase. Consistently, when $r_s \to 1$, i.e., when the limit case of the spherical shape is reached, the classical value of $\varphi^c=1/2$ is retrieved. On Fig. 5, the critical solid volume fraction $f_s^c=1-\varphi^c$ is represented as a function of the aspect ratio $r_s$ of the spheroids, it is indicated

<table><caption>**Table 1.** Critical Porosity and Effective Poisson's Ratio Associated, for Some Values of the Aspect Ratio of the Spheroids</caption><thead><tr><th>r<sub>s</sub></th><th>0</th><th>1</th><th>&infin;</th></tr></thead><tbody><tr><th>&phi;<sup>c</sup></th><td>1</td><td>1/2</td><td>1</td></tr><tr><th>&nu;<sup>c</sup></th><td>&asymp;0.173</td><td>1/5</td><td>(17&minus;&radic;79)/35 &asymp; 0.232</td></tr></tbody></table>

<table><caption>**Table 2.** Polynomials $a_{i_k, i_g}(m_s, \varphi)$ Occurring in Eq. (7)</caption><thead><tr><th></th><th colspan="3" style="text-align: center;">i<sub>k</sub></th></tr><tr><th style="text-align: center;">i<sub>g</sub></th><th style="text-align: center;">0</th><th style="text-align: center;">1</th><th style="text-align: center;">2</th></tr></thead><tbody><tr><th style="text-align: center;">0</th><td style="text-align: center;">0</td><td style="text-align: center;">0</td><td style="text-align: center;">(m<sub>s</sub>+4)&phi;</td></tr><tr><th style="text-align: center;">1</th><td style="text-align: center;">m<sub>s</sub><sup>2</sup>(&phi;-1)</td><td style="text-align: center;">m<sub>s</sub>(m<sub>s</sub>+8&phi;-4)</td><td style="text-align: center;">m<sub>s</sub>(4-&phi;)</td></tr><tr><th style="text-align: center;">2</th><td style="text-align: center;">3m<sub>s</sub><sup>2</sup>(&phi;-1)</td><td style="text-align: center;">3m<sub>s</sub><sup>2</sup></td><td style="text-align: center;">0</td></tr></tbody></table>

<table><caption>**Table 3.** Polynomials $b_{i_k, i_g}(m_s, \varphi)$ Occurring in Eq. (8)</caption><thead><tr><th></th><th colspan="3" style="text-align: center;">i<sub>k</sub></th></tr><tr><th style="text-align: center;">i<sub>g</sub></th><th style="text-align: center;">0</th><th style="text-align: center;">1</th><th style="text-align: center;">2</th></tr></thead><tbody><tr><th style="text-align: center;">0</th><td style="text-align: center;">0</td><td style="text-align: center;">0</td><td style="text-align: center;">48(1-&phi;)</td></tr><tr><th style="text-align: center;">1</th><td style="text-align: center;">0</td><td style="text-align: center;">116m<sub>s</sub>(1-&phi;)</td><td style="text-align: center;">4m<sub>s</sub>(36-61&phi;)+16(27-52&phi;)</td></tr><tr><th style="text-align: center;">2</th><td style="text-align: center;">56m<sub>s</sub><sup>2</sup>(1-&phi;)</td><td style="text-align: center;">m<sub>s</sub>[m<sub>s</sub>(312-587&phi;)+4(207-482&phi;)]</td><td style="text-align: center;">4m<sub>s</sub>(84-209&phi;)-16(47&phi;+3)</td></tr><tr><th style="text-align: center;">3</th><td style="text-align: center;">m<sub>s</sub><sup>2</sup>[m<sub>s</sub>(144-319&phi;)+20(18-53&phi;)]</td><td style="text-align: center;">m<sub>s</sub>[35m<sub>s</sub>(16-51&phi;)-4(299&phi;+101)]</td><td style="text-align: center;">-4m<sub>s</sub>(139&phi;+36)+16(2&phi;-27)</td></tr><tr><th style="text-align: center;">4</th><td style="text-align: center;">m<sub>s</sub><sup>2</sup>[m<sub>s</sub>(224-949&phi;)-8(69&phi;+31)]</td><td style="text-align: center;">m<sub>s</sub>[-m<sub>s</sub>(869&phi;+456)+20(2&phi;-27)]</td><td style="text-align: center;">12m<sub>s</sub>(3&phi;-28)</td></tr><tr><th style="text-align: center;">5</th><td style="text-align: center;">m<sub>s</sub><sup>2</sup>[-5m<sub>s</sub>(77&phi;+48)+4(17&phi;-42)]</td><td style="text-align: center;">m<sub>s</sub><sup>2</sup>(41&phi;-416)</td><td style="text-align: center;">0</td></tr><tr><th style="text-align: center;">6</th><td style="text-align: center;">m<sub>s</sub><sup>3</sup>(53&phi;-128)</td><td style="text-align: center;">0</td><td style="text-align: center;">0</td></tr></tbody></table>

by a solid line. Starting from a maximum related to spherical solid phase shape, the critical volume fraction decreases when the spheroids are either lengthened or flattened. Indeed, it seems intuitive that upon lengthening or flattening, less solid is required to reach a given stiffness level. As already pointed out in (Fritsch et al. 2006), the critical porosity disappears ($\varphi^c \to 1$) if $r_s \to \infty$ (see Table 1). This is also the case if $r_s \to 0$. Table 1 gives the effective Poisson's ratio showing up at the critical porosity $\nu^c=\lim_{\varphi \to \varphi^c} \nu^{\text{SCS}}$. The latter only depends on the aspect ratio $r_s$. It is independent of the solid Poisson's ratio $\nu_s$.

The fact that the critical porosity depends exclusively on the aspect ratio $r_s$, which is, besides the porosity, the only other morphological parameter of the model, suggests that the critical porosity might reflect a purely geometrical phenomenon. Garboczi et al. (1995) performed numerical simulations of percolation thresholds of intersectable spheroids. These thresholds, obtained from a purely geometrical method, are represented by dots on Fig. 5. The percolation thresholds (from geometry) and the critical porosities (from mechanics) are found to describe remarkably similar tendencies. This strengthens the relevance of the proposed self-consistent scheme, for describing porous polycrystals made up of elongated particles.

## Limit Case of Needles

In the section devoted to strength, we will make use of the self-consistent elastic moduli in the case of needles ($r_s \to \infty$) whose elasticity tensor is isotropic (bulk and shear moduli $k_s$ and $g_s$). In this case, relatively simple scalar equations can be extracted from Eq. (6)

$$ \sum_{i_k=0}^{2} \sum_{i_g=0}^{2} a_{i_k, i_g}(m_s, \varphi) \left( \frac{k^{\text{SCS}}}{k_s} \right)^{i_k} \left( \frac{g^{\text{SCS}}}{g_s} \right)^{i_g} = 0 \quad (7) $$

$$ \sum_{i_k=0}^{2} \sum_{i_g=0}^{6} b_{i_k, i_g}(m_s, \varphi) \left( \frac{k^{\text{SCS}}}{k_s} \right)^{i_k} \left( \frac{g^{\text{SCS}}}{g_s} \right)^{i_g} = 0 \quad (8) $$

with $m_s=4g_s/3/k_s$, and with the polynomials $a_{i_k, i_g}(m_s, \varphi)$ and $b_{i_k, i_g}(m_s, \varphi)$, respectively, being given in Tables 2 and 3, respectively.

## Micromechanics of Gypsum Elasticity

Being validated against numerical computations on microstructures made up of elongated particles, the modified self-consistent scheme proposed in the previous section can confidently be used
<table><caption>**Table 4.** Nonzero Components of the Stiffness Tensor of the Gypsum Crystal (Haussühl 1965)</caption><thead><tr><th>ijkl</th><th>c<sub>ijkl</sub><br>(GPa)</th></tr></thead><tbody><tr><td>1111</td><td>78.59</td></tr><tr><td>2222</td><td>62.74</td></tr><tr><td>3333</td><td>72.59</td></tr><tr><td>1122</td><td>41.01</td></tr><tr><td>1133</td><td>26.85</td></tr><tr><td>2233</td><td>24.24</td></tr><tr><td>2323</td><td>9.10</td></tr><tr><td>1313</td><td>26.41</td></tr><tr><td>1212</td><td>10.44</td></tr><tr><td>1113</td><td>-7.0</td></tr><tr><td>2213</td><td>3.1</td></tr><tr><td>3313</td><td>-17.4</td></tr><tr><td>2312</td><td>-1.55</td></tr></tbody></table>

to explain the stiffness of gypsum. Thereby, the model input parameters are: (1) porosity; (2) aspect ratio of the prolates; and (3) elasticity tensor of single gypsum crystal.

## Stiffness of the Gypsum Crystals

### Elasticity Tensor of a Single Gypsum Crystal

Each single gypsum crystal exhibits an anisotropic elastic behavior, arising from its crystallographic structure. Its stiffness tensor has been measured by ultrasonic means (Haussühl 1965). The nonzero components of the latter are listed in Table 4, which is to be completed by symmetry ($c_{ijkl}=c_{jikl}=c_{iklj}=c_{klij}$).

These components are given in a particular frame, defined with respect to the elementary crystallographic lattice of De Jong and Bouman (1938) (see Fig. 6).

$$ \epsilon_1 = \frac{a}{|a|}, \quad \epsilon_2 = \frac{b}{|b|} \quad \text{and} \quad \epsilon_3 = \epsilon_1 \wedge \epsilon_2 \qquad (9) $$

### Axis of the Crystal in the Elementary Lattice

Next, the position of the single crystal's principal material directions relative of the direction of the longest dimension of the elongated (axisymmetric) crystal is due. In other words, the axes of De Jong and Bouman's lattice need to be properly oriented with respect to the longest dimension of the elongated crystal. In this aspect, two contradictory options are reported in the open literature:

Observing the growth of single gypsum crystals, Klima (1983) noted that $l_{\text{meas}} \gg w_{\text{meas}}$ (see Fig. 7). Accordingly, the direction

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Crystallographic lattice of De Jong and Bouman (1938) and definition of the frame ($\epsilon_1, \epsilon_2, \epsilon_3$) in which the gypsum stiffness tensor is expressed (see Table 4)</figcaption></figure>

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

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

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Schematic representation of growing crystals, as observed by Klima (1983)</figcaption></figure>

into which the crystals develop their largest length is the intersection of the crystallographic planes (010) and (011), that is $a$. For the appropriate representation of this case, we will consider the vector $a$ of the lattice of De Jong and Bouman as the axis of the spheroids representing the solid phase.

On the other hand, according to Simon and Bienfait (1965) [see the detailed review in Meille (2001)], the gypsum crystals grow preferentially in the direction $c$ of the elementary lattice of De Jong and Bouman (see Fig. 8). We will not enter into the discussion on which of these two options should be preferred. Instead, we consider also an appropriate mechanical representation of the second option, by introducing the vector $c$ of the lattice of De Jong and Bouman as the axis of symmetry of the axisymmetric solid spheroids. In the following, we will compare model-predicted effective stiffnesses corresponding to the aforementioned two options.

## Experimental Validation of the Microelastic Model for Gypsum

Herein, we validate the model described and validated in the previous section, with the solid stiffness tensor given in Table 4, by comparing respective self-consistent estimates evaluated for three aspect ratios of the solid spheroids, (namely, $r_3=7, 10, \text{ and } 15$), to corresponding experimental data for gypsum [dried samples, that is with empty (drained) pores], reported in Meille (2001), Colak (2006), Ali and Singh (1975), Phani (1986), and Tazawa (1998)

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Schematic representation of a gypsum crystal with its faces and the direction of largest length (Meille 2001)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Experimental data for Young's modulus of dried gypsum (Meille 2001; Colak 2006; Ali and Singh 1975; Phani 1986; Tazawa 1998), and corresponding self-consistent estimates, evaluated for two lattice orientations with respect to the axes of the single prolate crystals, and for three aspect ratios of the solid spheroids</figcaption></figure>

(see Fig. 9). This figure shows two plots: the left one relates to the hypothesis of Klima (1983) where the direction of the largest dimension of each single crystal coincides with the lattice axis *a*, while the right one relates to the hypothesis of Simon and Bienfait (1965) where the direction of the largest dimension of each single crystal coincides with the lattice axis *c* (see previous subsecession "Axis of the Crystal in the Elementary Lattice").

A satisfying agreement between model predictions and experiments is reached for an aspect ratio of 15, over the whole range of experimentally investigated porosities. Note that the choice of the axis of the prolates, either *a* or *c*, hardly affects the micromechanical prediction of the effective Young's modulus as a function of the porosity (see also the two virtually superimposed curves of Fig. 10). As regards the overall homogenized Young's modulus of the polycrystals, the relation between lattice and crystal orientations is eventually negligible. This allows us, for the sake of economy, to restrict, from now onwards, our considerations to one of the two orientation options; we will consider the crystal's long axis being aligned with the lattice axis *c*.

Conclusively, the self-consistent scheme with prolate-shaped solid inhomogeneities thus provides a morphological representation which is well suited to predict the effective Young's modulus of materials made up of lengthened crystals, such as gypsum. For the sake of completeness, let us state that the model's limit case of needles ($r_s \to \infty$) has already been successfully applied to hydroxyapatite biomaterials (Fritsch et al. 2006, 2009).

## Influence of the Aspect Ratio of the Particles

The conditions for preparing the mix between plaster powder and water can have a great impact on the shape of the crystals ob-

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Young's modulus estimated by self-consistent scheme, comparing two hypotheses on the orientation of the axis of the prolates with respect to the lattice of De Jong and Bouman ($r_s=15$)</figcaption></figure>

<table><caption>**Table 5.** Effect of the Preparation-Specific Average Aspect Ratio of the Single Crystals Found in One Material Sample, on the Effective Young's Modulus of Gypsum (Meille 2001)</caption><thead><tr><th>Mixing condition</th><th>Aspect ratio</th><th>Young's modulus (GPa)</th></tr></thead><tbody><tr><td>Reference</td><td>12.5</td><td>3.7</td></tr><tr><td>Crushed gypsum</td><td>13.3</td><td>4.5</td></tr><tr><td>70°C</td><td>9.8</td><td>4.2</td></tr><tr><td>Citric acid</td><td>0.5</td><td>1.3</td></tr></tbody></table>

tained at the end of hydration. Considering a constant water over plaster ratio and using various mixing conditions, Meille (2001) prepared gypsum samples with crystals of different aspect ratios, all exhibiting the same intercrystalline porosity $\varphi=0.55$. Sample-specific spatial "averages" over the aspect ratios of the individual single crystals were estimated from SEM observations (see Table 5 for the average values and the corresponding Young's moduli). Conclusively, experiments show that the aspect ratio of the single crystals modifies the effective Young's modulus of the overall polycrystal-type material. Subsequently, we evaluate whether the proposed self-consistent scheme is able to predict this tendency.

The standard mixing conditions (mixing with tap water without any additive at room temperature) are taken as a reference (first line of Table 5). Therefore, in Fig. 11, the model-predicted Young's moduli are normalized with respect to the Young's modulus corresponding to the aforementioned reference case [which is $E_{ref}^{eff}=E^{eff}(r_s=12.5)$], and these predictions are studied as function of the aspect ratio (see Fig. 11). The model-predicted stiffness estimates agree reasonably well with the measured data, in particular when considering the experimental scattering of both the Young's modulus data [dispersion of 10% (Meille 2001)] and of the aspect ratio data (see Fig. 11). The (solid) curve representing model-predicted stiffnesses falls into three parts, which are separated by values for the aspect ratio which satisfy $\varphi^c(r_s)=\varphi^c$, whereby $\varphi^c$ is the critical porosity defined in subsection "critical porosity." Visualizing this situation, the dashed line of Fig. 11 represents the curve $\varphi^c(r_s)$, the horizontal dotted line represents $\varphi=0.55$, and the two vertical lines correspond to the two aspect ratios where $\varphi^c(r_s)=\varphi$. Let us shortly discuss these three parts of the solid curve illustrated in Fig. 11:

* For flattened particles (oblate spheroids $r_s<1$), the Young's modulus is the higher the flatter the particles (the smaller $r_s$). Once $r_s$ becomes so large that $\varphi^c(r_s)=\varphi$, the effective stiffness vanishes.

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Normalized Young's modulus as a function of the aspect ratio of the crystals, measured (Meille 2001) and estimated by a self-consistent scheme ($\varphi=0.55$); $\varphi^c(r_s)$ is the critical porosity beyond which the overall effective stiffness vanishes</figcaption></figure>
* Conversely, for elongated particles (prolate spheroids $r_s > 1$), the Young's modulus is the higher the more elongated the particles (the larger $r_s$). For example, for $r_s = 5$, the Young's modulus only reaches 45% of the reference value. Once $r_s$ becomes so small that $\phi^c(r_s) = \phi$, the effective stiffness vanishes.

* Third, there is a region around $r_s = 1$ where the particles are not elongated or flattened enough as to allow for any nonzero stiffness.

The model suggests that lengthened or flattened particles are desirable for optimization of the effective stiffness. However, the effective Young's modulus reaches a plateau at 120% of the reference value ($r_s = 12.5$) when $r_s \to \infty$, which renders the reference case as already being quite well optimized. Practically, it is difficult to produce very elongated crystals. Moreover, bundles of elongated crystals have sometimes been observed; from a viewpoint of precise representation of the microstructure, this could suggest consideration of nonisotropic and nonrandom distribution of crystals, and/or of a morphological building block comprising these bundles of crystals.

It is important to note that to plot the continuous curve on Fig. 11, the aspect ratio $r_s$ has been used as the aspect ratio of the solid inhomogeneities entering the auxiliary elasticity problems of the self-consistent scheme. The abscissa of each experimental point (dot on Fig. 11) corresponds to an average aspect ratio measured on the real microstructure. From a theoretical point of view, the aspect ratio of the solid spheroids used in the auxiliary problems of elasticity need not necessarily to be the average aspect ratio of the particles making up the microstructure. This suggests some care in interpreting the comparison proposed in Fig. 11 (see for example the data point at $r_s = 0.5$ which falls in the domain of vanishing effective modulus predicted by the model). In addition, the random and isotropic distribution of crystals is certainly a simplified representation of the real microstructure of gypsum.

## “Equivalent” Isotropic Stiffness Tensor

All stiffness estimations of the two previous subsections were based on the anisotropic stiffness tensor of the single gypsum crystals (see Table 4). The related complexity and computational expenditure is tractable for purely elastic calculations, but it becomes a serious issue when entering questions of elastic limit and strength. Therefore, the question arises whether this anisotropic single crystal elasticity can be represented in a computationally less expensive way, namely through an *equivalent* isotropic elasticity tensor characterizing single crystal behavior. Such an approach has already been successfully in the framework of finite element based numerical homogenization schemes (Meille and Garboczi 2001). Also for porous (polycrystal) hydroxyapatite biomaterials, anisotropic, and equivalent isotropic hydroxyapatite stiffness tensors yielded virtually the same effective stiffnesses estimated from needle-based self-consistent schemes (Fritsch et al. 2006).

As equivalent isotropic stiffness of the single crystals, we propose to use the homogenized isotropic stiffness of gypsum, evaluated for the porosity tending toward zero (see Fig. 9). This yields $E^{\text{eff}}(\varphi \to 0^+) = 40$ GPa and $\nu^{\text{eff}}(\varphi \to 0^+) = 0.34$. The relevance of such an approximation is evidenced in Fig. 12, which depicts the effective Young's modulus as a function of porosity, using prolates with the anisotropic stiffness on the one hand (computations from subsection “Experimental Validation of the Microelastic Model For Gypsum,” plain line), and with isotropic stiffness on

the other (dashed line); the differences are negligible.

Thus, in the following, we will consider the stiffness of the single gypsum crystals as isotropic, characterized by the Young's modulus $E_s = 40$ GPa and the Poisson's ratio $\nu_s = 0.34$.

## Modeling Brittle Strength of Gypsum through an Elastic Limit Approach

In the previous section, a micromechanical model of gypsum elasticity has been proposed. The good agreement with experimental data supports the fact that the chosen morphological model is relevant with respect to the gypsum microstructure. The purpose of the present section is to reuse the morphological model, in order to investigate the strength of gypsum.

The failure mechanism of dry gypsum is brittle and seems to involve decohesion of the crystals (Meille 2001). Precise modeling of such an interfacial phenomenon is feasible for the case of spherical inclusions (Dormieux et al. 2007; Fritsch et al. 2007a), but the extension to prolate particles would involve mathematical expenditures which go far beyond the scope of this paper. Therefore, we extend our fundamental vision of the interfacial decohesion of needle-type objects, remembering that such objects are predominantly loaded through stresses related to the normal pointing in the needle direction. Remembering in addition that our approach delivers (average) stress fields which are homogeneous within the individual needles, we propose to formulate a failure criterion in terms of the stress vectors related to the needle direction (being homogeneous within the corresponding cylindrical phase), using therefore a bulk failure concept for a virtually one-dimensional (needle) structure, as to mimic interfacial failure. This concept might sound convincing from the viewpoint of engineering intuition, but daring from a strictly mathematical viewpoint. Thus, there is only one way to show its practical relevance, careful experimental validation. Before, we have to develop the theory for upscaling of brittle strength from the single needle scale to that of the porous polycrystal. Given the brittle nature of the material, its strength coincides with its elastic limit. At the individual needle scale, this limit is described through a microscopic criterion describing the set of stresses admissible within a single gypsum crystal. The failure of the polycrystal, in turn, is considered to occur, once one crystal reaches this failure criterion. Thus, it is the most adversely stressed crystal direction facing failure. In this sense, failure is indeed related to “stress peaks” in the porous gypsum polycrystal, although we employ average stresses over single crystals exhibiting the same orientation.

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Young's modulus of gypsum, from anisotropic or isotropic prolates ($r_s = 15$)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Decomposition of the stress vector acting on the cross section of a crystal</figcaption></figure>

## Criterion at the Crystals Scale

The microscopic failure criterion depends on the average stress in the crystal. Moreover, as the gypsum crystals exhibit an elongated shape, it seems to be reasonable to take into account the morphological prominence of one direction. The average stress in a crystal of axis $\mathcal{N}$ is denoted by $\boldsymbol{\sigma}^N$. We thus propose a microscopic criterion based on the stress vector acting on the cross section of the crystal, that is $\boldsymbol{\sigma}^N \cdot \mathcal{N}$. This stress vector can be decomposed into an axial stress $\sigma_N = \mathcal{N} \cdot \boldsymbol{\sigma}^N \cdot \mathcal{N}$ and a tangential stress vector $\sigma_T = \boldsymbol{\sigma}^N \cdot \mathcal{N} - \sigma_N \mathcal{N}$ (see Fig. 13) with magnitude $\sigma_T = |\sigma_T|$. We consider $\sigma_N$ and $\sigma_T$ fulfilling a Mohr-Coulomb criterion

$$ \sigma_T = \beta(\sigma_{cr} - \sigma_N) \quad (10) $$

with two parameters $\sigma_{cr}$ and $\beta$ characterizing the microscopic strength (see Fig. 14). Brittle failure is considered; the macroscopic criterion is reached once one single crystal fails.

## Homothetic (Cone-Type) Nature of Macroscopic Failure Criterion

This section is devoted to the investigation of the qualitative shape of the macroscopic elastic domain, when the failure criterion of the crystals is described by Eq. (10). The macroscopic elastic domain can be written in the form

$$ \{ \Sigma \text{ such that } \sup_{(N, T), \mathcal{N} \cdot \mathcal{T} = 0, |\mathcal{N}| = |\mathcal{T}| = 1} \beta \sigma_N + \sigma_{NT} < \beta \sigma_{cr} \} \quad (11) $$

with $\mathcal{N}$ the axis of the crystal, $\sigma_N = \mathcal{N} \cdot \boldsymbol{\sigma}^N \cdot \mathcal{N}$ and $\sigma_{NT} = \mathcal{N} \cdot \boldsymbol{\sigma}^N \cdot \mathcal{T}$.

## Localization of Average Stress in Crystals

To express brittle failure in terms of the macroscopic stress tensor $\Sigma$, a relation between the stress $\boldsymbol{\sigma}^N$ in the crystals and the macroscopic stress is due. Such a relation is standardly called localization or concentration relation, and expressed through a stress concentration tensor $\mathbb{B}$, according to

$$ \boldsymbol{\sigma}^N = \mathbb{B}(\mathcal{N}, v_s, r_s, \varphi) : \Sigma \quad (12) $$

The stress concentration tensor obviously depends on the orientation $\mathcal{N}$ of the crystal. Moreover, as dimensional analysis shows, it also depends on the Poisson's ratio $v_s$ of the crystals, on the aspect ratio $r_s$ of the latter, and on the porosity $\varphi$. In the following, we omit the dependence on the three parameters, simply denoting the concentration tensor as $\mathbb{B}(\mathcal{N})$. Considering $\mathcal{N} = \varepsilon_r$, the expression for $\mathbb{B}(\mathcal{N})$ can be derived from Eqs. (2) and (4), and (5)

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Elastic domain for gypsum crystals</figcaption></figure>

$$ \mathbb{B}(\mathcal{N}) = \frac{1}{1-\varphi} C_s: [1 + \mathbb{P}_{\text{pro}}^{\text{SCS}}(\theta, \phi):(C_s - \mathbb{C}^{\text{SCS}})]^{-1} \times \left\{ \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} [1 + \mathbb{P}_{\text{pro}}^{\text{SCS}}(\theta, \phi):(C_s - \mathbb{C}^{\text{SCS}})]^{-1} \right\}^{-1} : C_s^{-1} \quad (13) $$

Note that the stiffness tensor $C_s$ of the gypsum crystals is now regarded as isotropic (see subsection “Equivalent Isotropic Stiffness Tensor”). In Eq. (13), every tensor is isotropic, except for $\mathbb{P}_{\text{pro}}^{\text{SCS}}(\theta, \phi)$, which is transverse isotropic, with $\mathcal{N}$ being the normal to the isotropic plane. In the frame $(\varepsilon_\theta, \varepsilon_\phi, \varepsilon_r)$, the components of $\mathbb{P}_{\text{pro}}^{\text{SCS}}(\theta, \phi)$ do not depend on $(\theta, \phi)$. This independence also holds for the stress localization tensor $\mathbb{B}(\mathcal{N})$.

## Case I: Hydrostatic (Isotropic) Macroscopic Stress

To begin with, a purely spherical (isotropic and hydrostatic) macroscopic stress state $\Sigma = \mathbf{h} \mathbf{1}$ is considered. We seek the maximum value of $\mathbf{h}$ compatible with the failure criterion [Eq. (10)]. Since $\mathbb{B}(\mathcal{N})$ is a transversely isotropic fourth-order tensor with the isotropic plane perpendicular to $\mathcal{N}$, it follows that

$$ \mathbb{B}(\mathcal{N}):\mathbf{1} = b_N \mathcal{N} \otimes \mathcal{N} + b_T(\mathcal{T}_1 \otimes \mathcal{T}_1 + \mathcal{T}_2 \otimes \mathcal{T}_2) \quad (14) $$

with $\mathcal{T}_1 \cdot \mathcal{N} = 0$, $\mathcal{T}_2 \cdot \mathcal{N} = 0$, and $\mathcal{T}_1 \cdot \mathcal{T}_2 = 0$. It is emphasized that $b_N$ and $b_T$ do not depend on the orientation $\mathcal{N}$ of the crystal. The concentration rule $\boldsymbol{\sigma}^N = \mathbf{h} \mathbb{B}(\mathcal{N}):\mathbf{1}$ yields $\sigma_N = \mathbf{h} b_N$ and $\sigma_{NT} = 0$. Accordingly, for an isotropic macroscopic stress, the inequality [Eq. (11)] which characterizes the admissible macroscopic stress states, takes the form

$$ \sup_{(N, T), \mathcal{N} \cdot \mathcal{T} = 0, |\mathcal{N}| = |\mathcal{T}| = 1} \beta \sigma_N + \sigma_{NT} = \beta h b_N < \beta \sigma_{cr} \quad (15) $$

In other words, the local criterion proves to be independent of the orientation of the crystal. This is clearly due to the fact that an isotropic distribution of crystal orientation has been assumed, together with an isotropic macroscopic loading. Eq. (15) implies that the maximum value $h^{\text{max}}$ of $\mathbf{h}$, chosen such that $\mathbf{h} \mathbf{1}$ be admissible, reads $h^{\text{max}} = \sigma_{cr} / b_N$. In the limit case of needle-shaped crystals (aspect ratio $r_s \to \infty$), the expression of $h^{\text{max}}$ becomes analytical

$$ h^{\text{max}} = \frac{(1 - \varphi)k_s(3k^{\text{SCS}} + 3g^{\text{SCS}} + g_s)}{3k_s(k^{\text{SCS}} + g^{\text{SCS}}) + g_s(3k_s - 2k^{\text{SCS}})} \sigma_{cr} \quad (16) $$

Here, we recall that moduli $k^{\text{SCS}}$ and $g^{\text{SCS}}$ are numerically determined from Eqs. (7) and (8). From Eq. (16), using $k^{\text{SCS}} < k_s$, it is readily verified that $h^{\text{max}} \le \sigma_{cr}$, the equality being achieved in the asymptotic situation $\varphi = 0$ (pure gypsum). Dimensional analysis of Eq. (16) also reveals that the normalized tensile elastic limit $h^{\text{max}} / \sigma_{cr}$ only depends on $\varphi$ and $v_s$. For a constant $v_s = 0.34$, its dependence on the porosity $\varphi$ is presented in Fig. 15.
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Macroscopic elastic limit in hydrostatic tension, as a function of porosity, for needle-shaped crystals ($r_s \to \infty$) and $v_s=0.34$</figcaption></figure>

### Case II: Arbitrary Macroscopic Stress

The purpose of this section is to establish that the elastic domain $\mathcal{E}$ in the macroscopic stress space is a cone which apex lies on the straight line of isotropic stresses, namely, at point $\Sigma=h^{\max}\mathbf{1}$—and this holds for any aspect ratio of the solid inclusions.

To do so, let us assume that $\Sigma=h^{\max}\mathbf{1}+S \in \mathcal{E}$, $S$ denoting a given symmetric second-order tensor. We have to show that $\Sigma=h^{\max}\mathbf{1}+\lambda S \in \mathcal{E}$, $\forall \lambda > 0$. Owing to the fact that $\Sigma=h^{\max}\mathbf{1}+S \in \mathcal{E}$, we have

$$ \begin{gathered} \forall (\underline{N}, \underline{T}) \text{ such that } \underline{N} \cdot \underline{T} = 0 \text{ and } |\underline{N}| = |\underline{T}| = 1 \\ \beta \underline{N} \cdot [\underline{B}(\underline{N}):(h^{\max} \mathbf{1} + S)] \cdot \underline{N} + \underline{N} \cdot [\underline{B}(\underline{N}):(h^{\max} \mathbf{1} + S)] \cdot \underline{T} < \beta \sigma_{cr} \end{gathered} \quad (17) $$

Recalling that $h^{\max}$ meets the condition $\underline{N} \cdot [\underline{B}(\underline{N}):h^{\max}\mathbf{1}] \cdot \underline{N} = \sigma_{cr}$ and that $\underline{N} \cdot [\underline{B}(\underline{N}):h^{\max}\mathbf{1}] \cdot \underline{T} = 0$ [see Eq. (14)], Eq. (17) becomes

$$ \begin{gathered} \forall (\underline{N}, \underline{T}) \text{ such that } \underline{N} \cdot \underline{T} = 0 \text{ and } |\underline{N}| = |\underline{T}| = 1 \\ \beta \underline{N} \cdot [\underline{B}(\underline{N}):\underline{S}] \cdot \underline{N} + \underline{N} \cdot [\underline{B}(\underline{N}):\underline{S}] \cdot \underline{T} < 0 \end{gathered} \quad (18) $$

Multiplying Eq. (18) by $\lambda > 0$ and adding $\beta \underline{N} \cdot [\underline{B}(\underline{N}):h^{\max}\mathbf{1}] \cdot \underline{N} = \beta \sigma_{cr}$ and $\underline{N} \cdot [\underline{B}(\underline{N}):h^{\max}\mathbf{1}] \cdot \underline{T} = 0$ eventually yields

$$ \begin{gathered} \forall (\underline{N}, \underline{T}) \text{ such that } \underline{N} \cdot \underline{T} = 0 \text{ and } |\underline{N}| = |\underline{T}| = 1 \\ \beta \underline{N} \cdot [\underline{B}(\underline{N}):(h^{\max} \mathbf{1} + \lambda \underline{S})] \cdot \underline{N} + \underline{N} \cdot [\underline{B}(\underline{N}):(h^{\max} \mathbf{1} + \lambda \underline{S})] \cdot \underline{T} < \beta \sigma_{cr} \end{gathered} \quad (19) $$

This inequality states that the macroscopic stress state $\Sigma=h^{\max}\mathbf{1}+\lambda S \in \mathcal{E}$. This completes the proof that $\mathcal{E}$ is a cone of the stress space with apex $h^{\max}\mathbf{1}$. We now focus on the shape of the cross section of the cone.

### Cross-Sectional Shape of Macroscopic Failure Criterion

The purpose of this section is to characterize the shape of the cone representing the elastic domain $\mathcal{E}$ in the macroscopic stress space. Assuming an isotropic distribution of crystals orientation, $\mathcal{E}$ can be sought in the space of principal stresses ($\Sigma_1, \Sigma_2, \Sigma_3$). We therefore write the macroscopic stress state in the form

$$ \Sigma = \Sigma_1 \varepsilon_1 \otimes \varepsilon_1 + \Sigma_2 \varepsilon_2 \otimes \varepsilon_2 + \Sigma_3 \varepsilon_3 \otimes \varepsilon_3 \quad (20) $$

For the sake of simplicity, we restrict the analysis to the case of needle-shaped crystals ($r_s \to \infty$).

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Porosity-dependent coefficients $a_0$ and $a_1$ found in the localization Eq. (23), giving access to (microscopic) normal stress $\sigma_N(v_s=0.34)$</figcaption></figure>

### Case I: Failure of Single Crystals Governed by Tensile Stress Only

We first consider the case where the failure of the elementary crystal is controlled by the normal stress

$$ \sigma_N < \sigma_{cr} \quad (21) $$

The first step consists in relating the macroscopic stress $\Sigma$ to the normal (microscopic) stress in the crystal oriented in direction $\underline{N}=\underline{e}_r$, this microscopic normal stress being denoted by $\sigma_N=\underline{N} \cdot \underline{\sigma}^N \cdot \underline{N}$. Making therefore use of the stress concentration tensor defined by $\underline{B}$ [Eq. (12)], one obtains

$$ \sigma_N = \underline{N} \cdot [\underline{B}(\underline{N}):\underline{\Sigma}] \cdot \underline{N} \quad (22) $$

After some algebra, it is found that

$$ \begin{aligned} \sigma_N = a_0(\Sigma_1 + \Sigma_2 + \Sigma_3) &+ a_1[(\Sigma_1 \cos^2 \phi + \Sigma_2 \sin^2 \phi)\sin^2 \theta + \Sigma_3 \cos^2 \theta] \end{aligned} \quad (23) $$

where Eq. (13) has been used. The expressions of the coefficients $a_0$ and $a_1$ are too lengthy for being reported herein. Still, it is worth mentioning that they only depend on $v_s$ and $\phi$ with $a_1 > 0$. They are plotted in Fig. 16 as a function of porosity $\phi$.

Next, we are interested in the orientation of the single crystal which is the first to fulfill Eq. (21), i.e., of the single crystal which will fail first (critical crystal). Accordingly, we seek, for a given macroscopic stress state ($\Sigma_1, \Sigma_2, \Sigma_3$), the orientation ($\theta, \phi$) which maximizes $\sigma_N$ in the form Eq. (23). Since $a_1 > 0$, it suffices to maximize over $[0, \pi] \times [0, 2\pi]$ [the function $f(\theta, \phi)$ defined by

$$ f:(\theta, \phi) \mapsto (\Sigma_1 \cos^2 \phi + \Sigma_2 \sin^2 \phi)\sin^2 \theta + \Sigma_3 \cos^2 \theta \quad (24) $$

The maximum of $f(\theta, \phi)$ is found to depend on the way the principal stresses are ordered (see result in Table 6). As expected, the axis of the critical crystal coincides with the direction of the maximum principal stress (see Fig. 17). Moreover, the maximal value of $f$ is equal to this maximum principal stress $\Sigma_1$.

Let us order the principal stresses as $\Sigma_1 \ge \Sigma_{II} \ge \Sigma_{III}$. In the line of the previous result, the maximum of $\sigma_N$ Eq. (23) reads identically for all possible orientations of single crystals, namely

<table><caption>**Table 6.** Maximization of the Function $f$ Defined by Eq. (24)</caption><thead><tr><th rowspan="2">Case</th><th colspan="3">Maximum of f</th></tr><tr><th>θ</th><th>φ</th><th>Value of f</th></tr></thead><tbody><tr><td>Σ<sub>1</sub> &gt; Σ<sub>2</sub> and Σ<sub>1</sub> &gt; Σ<sub>3</sub></td><td>π/2</td><td>0</td><td>Σ<sub>1</sub></td></tr><tr><td>Σ<sub>2</sub> &gt; Σ<sub>1</sub> and Σ<sub>2</sub> &gt; Σ<sub>3</sub></td><td>π/2</td><td>π/2</td><td>Σ<sub>2</sub></td></tr><tr><td>Σ<sub>3</sub> &gt; Σ<sub>1</sub> and Σ<sub>3</sub> &gt; Σ<sub>2</sub></td><td>0</td><td></td><td>Σ<sub>3</sub></td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Axis of the critical crystal (represented by a bold arrow) depending on how the principal stresses are ordered</figcaption></figure>

$$ \max_{0 \le \theta \le \pi, 0 \le \phi < 2\pi} \sigma_N = a_0(\Sigma_1 + \Sigma_2 + \Sigma_3) + a_1 \Sigma_1 \quad (25) $$

Thus, when introducing the intensity of the spherical part of the macroscopic stress tensor, as $\Sigma_m=\text{tr}\Sigma/3$ (mean stress), the macroscopic failure criterion takes the form

$$ 3a_0\Sigma_m + a_1\Sigma_1 = \sigma_{cr} \quad (26) $$

where $a_0$ and $a_1$ have been plotted on Fig. 16. Specifying Eq. (26) for the stress state in the cone apex, i.e., for $\Sigma_1=\Sigma_2=\Sigma_3=h^{\text{max}}$, yields $h^{\text{max}}=\sigma_{cr}/(3a_0+a_1)$.

Note that $\varphi \to 0$ implies $a_0 \to 0$ and $a_1 \to 1$. Accordingly, the macroscopic criterion tends toward the Rankine-type criterion $\Sigma_1=\sigma_{cr}$, which can be also given in the form

$$ \Sigma_1 < \sigma_{cr} \text{ and } \Sigma_2 < \sigma_{cr} \text{ and } \Sigma_3 < \sigma_{cr} \quad (27) $$

This domain $\mathcal{E}_0$ corresponds to the eighth of space delimited by the three planes of equation $\Sigma_i=\sigma_{cr}, i \in \{1, 2, 3\}$.

In the general case $\varphi \ne 0$, it is found that $\mathcal{E} \subset \mathcal{E}_0$. This inclusion of the space of admissible (macroscopic) stress states within $\mathcal{E}_0$ is often introduced as a physically motivated requirement when the (macroscopic) failure criterion is chosen without microstructural considerations (i.e., in the framework of so-called phenomenological approaches), see Desrues (2002) for the case of the Drucker-Prager criterion. Interestingly, the present micromechanical approach meets this requirement naturally, irrespective of the chosen model parameter $\sigma_{cr}$.

The cross section of the macroscopic criterion is an equilateral triangle (Fig. 18), which equation is given by Eq. (26). Note that the macroscopic counterpart of a uniaxial (frictionless) micro-

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Macroscopic criterion depicted in the space of principal stresses, and cross section through the hydrostatic axis at $\Sigma_m = \text{tr}\Sigma/3$ constant</figcaption></figure>

scopic criterion $\sigma_N < \sigma_{cr}$ is a (conical) frictional criterion, that is, a criterion depending on the mean stress $\Sigma_m=\text{tr}\Sigma/3$.

To quantify the influence of porosity on the cone aperture, it useful to introduce two Drucker-Prager criteria [defined by $\Sigma_d = \alpha(h^{\text{max}}-\Sigma_m)$]: their cross sections are circles which are inscribed into $(\alpha_i)$ and circumscribed $(\alpha_c)$ around the equilateral triangle [Eq. (26)]. The corresponding friction coefficients $\alpha_i$ and $\alpha_c$ read

$$ \alpha_i = \frac{\sqrt{3}}{2} \left( 1 + \frac{3a_0}{a_1} \right) \quad \text{and} \quad \alpha_c = 2\alpha_i \quad (28) $$

The latter are plotted on Fig. 19 as a function of porosity: the higher the porosity, the sharper the cone.

## Case II: Failure of Single Crystals Governed by Normal and Shear Stress

We now consider the complete microscopic criterion [Eq. (10)]

$$ \sigma_N + \frac{1}{\beta} \sigma_T = \sigma_{cr} \quad (29) $$

with a finite value of $\beta$.

In addition to the normal stress $\sigma_N$ [Eq. (23)], also the shear stress $\sigma_T=|\boldsymbol{\sigma}^N \cdot \boldsymbol{N} - \sigma_N \boldsymbol{N}|$ needs to be related the macroscopic stress $\Sigma$. After some algebra, the corresponding concentration relation reads as

$$ \sigma_T = a_2 a_1 \sin \theta \sqrt{[\Sigma_3 - \Sigma_2 + (\Sigma_2 - \Sigma_1)\cos^2 \phi]^2 \cos^2 \theta + (\Sigma_2 - \Sigma_1)^2 \cos^2 \phi \sin^2 \phi} \quad (30) $$

with the coefficient $a_1>0$ already involved in $\sigma_N$ [Eq. (23)] and

$$ a_2 = \frac{4g_s g^{\text{SCS}}}{g_s + g^{\text{SCS}}} \frac{3k_s + g_s + 3g^{\text{SCS}}}{8g_s g^{\text{SCS}} + 9k_s g_s - 3k_s g^{\text{SCS}}} > 0 \quad (31) $$

Fig. 20 illustrates the latter coefficient as a function of the porosity.

The macroscopic loading is characterized by the stresses $(\Sigma_1, \Sigma_2, \Sigma_3)$. The critical crystal is the one which maximizes $\sigma_N+\sigma_T/\beta$. Referring to Eqs. (23) and (30), this expression reads

$$ \sigma_N + \frac{\sigma_T}{\beta} = a_0(\Sigma_1 + \Sigma_2 + \Sigma_3) + a_1 g(\theta, \phi) \quad (32) $$

with

$$ g(\theta, \phi) = (\Sigma_1 \cos^2 \phi + \Sigma_2 \sin^2 \phi)\sin^2 \theta + \Sigma_3 \cos^2 \theta + A \sin \theta \sqrt{[\Sigma_3 - \Sigma_2 + (\Sigma_2 - \Sigma_1)\cos^2 \phi]^2 \cos^2 \theta + (\Sigma_2 - \Sigma_1)^2 \cos^2 \phi \sin^2 \phi} \quad (33) $$
<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Friction coefficients $\alpha$ of the inscribed and circumscribed Drucker-Prager criteria ($v_s=0.34$)</figcaption></figure>

where $A=a_2/\beta > 0$. As $a_1 > 0$ and for symmetrical reasons, it suffices to maximize the function $g(\theta, \phi)$ on the domain $[0, \pi/2] \times [0, \pi/2]$, that is on one eighth of the unit sphere. It can be shown that this function does not admit a maximum on the open domain $]0, \pi/2[ \times ]0, \pi/2[$. A local maximum is thus sought on each of the “edges” of the eighth of the unit sphere, characterized by $\theta=\pi/2$, $\phi=0$, or $\phi=\pi/2$. The results are reported in Table 7. The global maximum is then determined. Eventually, the orientation of the critical crystal depends on both $\gamma_c=(\arctan A)/2$ and the way the principal stresses are ordered. Consequently, six cases are to be investigated (see Fig. 21).

As in the previous subsection, the principal stresses are ordered as $\Sigma_1 \ge \Sigma_{II} \ge \Sigma_{III}$. From Table 7, the maximum value of $g$ always reads as

$$ \max_{0 \le \theta \le \pi, 0 \le \phi < 2\pi} g(\theta, \phi) = \frac{\Sigma_1 + \Sigma_{III}}{2} + \frac{\Sigma_1 - \Sigma_{III}}{2} \sqrt{1 + A^2} \quad (34) $$

Consequently, when recalling Eq. (32) and the definition of $A$, the macroscopic elastic domain reads as

$$ 3a_0\Sigma_m + a_1 \left[ \frac{\Sigma_1 + \Sigma_{III}}{2} + \frac{\Sigma_1 - \Sigma_{III}}{2} \sqrt{1 + (a_2/\beta)^2} \right] < \sigma_{cr} \quad (35) $$

where $a_0, a_1, a_2$=still defined by Eqs. (23) and (31). Note that these three coefficients only depend on $v_s$ and $\phi$ (see Figs. 16 and 20). In the space of principal stresses, the criterion is a cone whose cross section is a nonregular hexagon (Fig. 22). The cone is the sharper, the lower the microscopic parameter $\beta$.

It is useful to derive the Drucker-Prager criteria inscribed into and circumscribed around the criterion defined by Eq. (35). The

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Coefficient $a_2$ [Eq. (31)] found in the localization Eq. (30) giving $\sigma_T$ ($v_s=0.34$)</figcaption></figure>

<table><caption>**Table 7.** Local Maximization of the Function $g$ Defined by Eq. (33)</caption><thead><tr><th colspan="3">Local maximum of g</th></tr><tr><th>case</th><th>φ</th><th>Value of g</th></tr></thead><tbody><tr><td>Σ<sub>1</sub> > Σ<sub>2</sub></td><td>(arctan A)/2</td><td>[Σ<sub>1</sub> + Σ<sub>2</sub> + (Σ<sub>1</sub> - Σ<sub>2</sub>)√(1 + A<sup>2</sup>)]/2</td></tr><tr><td>Σ<sub>2</sub> > Σ<sub>1</sub></td><td>(π - arctan A)/2</td><td>[Σ<sub>1</sub> + Σ<sub>2</sub> + (Σ<sub>2</sub> - Σ<sub>1</sub>)√(1 + A<sup>2</sup>)]/2</td></tr></tbody></table><table><thead><tr><th colspan="3">Local maximum of g</th></tr><tr><th>case</th><th>θ</th><th>Value of g</th></tr></thead><tbody><tr><td>Σ<sub>3</sub> > Σ<sub>1</sub></td><td>(arctan A)/2</td><td>[Σ<sub>1</sub> + Σ<sub>3</sub> + (Σ<sub>3</sub> - Σ<sub>1</sub>)√(1 + A<sup>2</sup>)]/2</td></tr><tr><td>Σ<sub>1</sub> > Σ<sub>3</sub></td><td>(π - arctan A)/2</td><td>[Σ<sub>1</sub> + Σ<sub>3</sub> + (Σ<sub>1</sub> - Σ<sub>3</sub>)√(1 + A<sup>2</sup>)]/2</td></tr></tbody></table><table><thead><tr><th colspan="3">Local maximum of g</th></tr><tr><th>case</th><th>θ</th><th>Value of g</th></tr></thead><tbody><tr><td>Σ<sub>3</sub> > Σ<sub>2</sub></td><td>(arctan A)/2</td><td>[Σ<sub>2</sub> + Σ<sub>3</sub> + (Σ<sub>3</sub> - Σ<sub>2</sub>)√(1 + A<sup>2</sup>)]/2</td></tr><tr><td>Σ<sub>2</sub> > Σ<sub>3</sub></td><td>(π - arctan A)/2</td><td>[Σ<sub>2</sub> + Σ<sub>3</sub> + (Σ<sub>2</sub> - Σ<sub>3</sub>)√(1 + A<sup>2</sup>)]/2</td></tr></tbody></table>

cone inscribed is tangent to the six planes delimiting the criterion. The circumscribed criterion is the one containing the half-line $\Sigma_1 = \Sigma_{II} \ge \Sigma_{III}$ of the criterion. Thus

$$ \alpha_i = \sqrt{3} \frac{1 + 3a_0/a_1}{\sqrt{1 + 3[1 + (a_2/\beta)^2]}} \quad \text{and} \quad \alpha_c = \sqrt{3} \left( 1 + \frac{3a_0}{a_1} \right) \quad (36) $$

These two parameters are represented as functions of the porosity $\phi$, in Fig. 23, for three values of $\beta$: 1.1, 1.0, and infinity [to retrieve the local criterion taking only the normal stress into account, given in Eq. (26)]. Note that $\beta=10$ nearly provides the results based on the purely tensile local criterion ($\beta \to \infty$). At high porosities ($\phi > 0.6$), the friction angles $\alpha_i$ and $\alpha_c$ hardly depend on $\beta$.

## Experimental Validation

In principal, it would be very desirable to have triaxial tests on gypsum, as to validate the micromechanical strength model developed in this subsection. Triaxial tests results exist, but on the gypsum rock, whose morphology (at the crystals scale) differ from the gypsum material studied here. Thus, we validate our model through uniaxial tensile and compressive strength tests. Again, we distinguish the case of a tensile failure criterion at the single crystal scale, given by Eq. (21) from that of the combined normal/shear stress criterion, given by Eq. (29).

<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Axis of the critical crystal (represented by a bold arrow) depending on the way the principal stresses are ordered</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 22.** Macroscopic criterion in the space of principal stresses and cross section through the hydrostatic axis at $Σ_m$=constant (the influence of the microscopic parameter $β$ is shown through three different values)</figcaption></figure>

## Case I: Failure of Single Crystals Governed by Tensile Stress Only

In a tensile (respectively compressive) test, the macroscopic stress state is defined by $Σ_1=Σ_r$, $Σ_{II}=Σ_{III}=0$, where $Σ_i>0$ (respectively $Σ_{III}=-Σ_c$ and $Σ_1=Σ_{II}=0$, where $Σ_c>0$). In these particular cases, the failure criterion [Eq. (26)] yields the following elastic limits

$$Σ_t^{\lim} = \frac{σ_{cr}}{a_0 + a_1} \quad \text{and} \quad Σ_c^{\lim} = -\frac{σ_{cr}}{a_0} \quad (37)$$

Apart from $φ$ and $v_s$, these elastic limits depend on the local normal stress threshold $σ_{cr}$. A direct measurement of this parameter, at the microscopic scale, seems to be out of reach. We consider the macroscopic tensile strength measured at a vanishing porosity, this strength amounting to around 15 MPa (see left part of Fig. 24). Since $Σ_t^{\lim}(φ→0)=σ_{cr}$ [see Eq. (37) for $a_0→0$ and $a_1→1$ according to Fig. 16], it follows that $σ_{cr}=15$ MPa. The model predictions are compared to experimental data on dried gypsum from various literature sources, illustrated by the left part of Fig. 24. The right-hand side of this figure illustrates the macroscopic strain (in the direction of loading) when the elastic limit is reached. For computing this failure strain, we consider that the macroscopic behavior remains linear elastic below the elastic limit.

The model predictions are satisfactory, except for the compressive strengths at low porosities. There, the shear stress at the single crystal scale obviously plays a role; it will be considered hereafter as Case II.

<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Friction angles $α$ of the Drucker-Prager criteria inscribed and circumscribed ($v_s=0.34$)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 24.** Uniaxial tensile and compressive elastic limits estimated ($v_s=0.34$ and $σ_{cr}=15$ MPa); tensile and compressive measured strengths [tensile experimental data from Ali and Singh (1975), Alsadi et al. (1996), Khalil (1972), and Sanad et al. (1982); and compressive experimental data from Ali and Singh (1975), Alsadi et al. (1996), Khalil (1972), Sanad et al. (1982), Gmouh et al. (2003), Tazawa (1998), Colak (2006), Schinkel (1993), and Fairhurst (1960)]; and macroscopic strain at the elastic limit</figcaption></figure>

## Case II: Failure of Single Crystals Governed by Normal and Shear Stress

From the criterion [Eq. (35)], the uniaxial macroscopic elastic limits are derived as

$$Σ_t^{\lim} = \frac{σ_{cr}}{a_0 + [\sqrt{1 + (a_2/β)^2} + 1]a_1/2}$$

and

$$Σ_c^{\lim} = \frac{σ_{cr}}{-a_0 + [\sqrt{1 + (a_2/β)^2} - 1]a_1/2} \quad (38)$$

The microscopic parameters $σ_{cr}$ and $β$ have to be estimated from experimental strength values at a vanishing porosity. Specifying Eq. (38) for $φ→0$, i.e., for $a_0→0$, $a_1→1$, and $a_2→1$, we have

$$Σ_t^{\lim}(φ → 0) = \frac{2σ_{cr}}{\sqrt{1 + 1/β^2 + 1}} \quad \text{and} \quad Σ_c^{\lim}(φ → 0) = \frac{2σ_{cr}}{\sqrt{1 + 1/β^2 - 1}} \quad (39)$$

At a quasizero porosity, the measured strength values in tension and in compression amount to 15 and 125 MPa (see experimental data plotted as dots on left part of Fig. 25). This yields $σ_{cr}=17$ MPa and $β=1.3$. The estimated elastic limits are then compared to the experimental measurements on the left part of Fig. 25. Now, the model is able to properly predict the strength in tension and in compression, across relatively dispersed experimental data coming from various sources.

The macroscopic strain reached at the elastic limit is plotted on the right part of Fig. 25 ($E_{33}^{\lim}=Σ_t^{\lim}/E^{SCS}$, for example, in tension). It is interesting to note that this strain at the elastic limit depends by far less on porosity (in tension the values stay in the interval $6.2×10^{-4}±39\%$, in compression the interval is $3.0×10^{-3}±12\%$), than does the stress at the elastic limit. This relative independence is confirmed by experimental results from (Dany 1995; Péronnet 1994) on dry gypsum, plotted as dots on the right part of Fig. 25. Knowing that this type of measurement is rather difficult to perform, it is interesting to note that the model predictions and the experimental data still fall into the same order of magnitude.
<figure><img src="image_28.png" /><figcaption>**Fig. 25.** Elastic limits in tension and compression estimated ($v_s$=0.34, $\sigma_{cr}$=17 MPa, and $\beta$=1.3); strength measured [experimental data from Ali and Singh (1975), Alsadi et al. (1996), Khalil (1972), and Sanad et al. (1982) for tension; and from Ali and Singh (1975), Alsadi et al. (1996), Khalil (1972), Sanad et al. (1982), Gmouh et al. (2003), Tazawa (1998), Colak (2006), Schinkel (1993), and Fairhurst (1960) for compression]; and strain at the elastic limit [experimental data from Dany (1995) and Péronnet (1994)]</figcaption></figure>

## Orientation of Critical (First Failing) Single Crystals

Finally, we are interested in the orientation of the first single crystal which fails, as function of different macroscopic stress states. We again distinguish two cases for the failure criterion at the scale of the single crystals.

### Case I: Failure of Single Crystals Governed by Tensile Stress Only

As readily seen from Fig. 17, the orientation of the critical crystal (that is to say the first crystal to reach its failure criterion) depends on the macroscopic stress state. In tension, the critical crystals are thus found to be parallel to the tension axis. In compression, the critical crystals are normal to the compression axis (see Fig. 26). In the latter case, there is an infinite amount of critical orientations characterized by the angles $\theta=\pi/2$ and $0 \le \phi < 2\pi$. These crystals appear to be in tension due to Poisson effects at the macroscopic scale. Our micromechanical analysis can thus provide a way to explain “axial splitting”—the cracking observed on rocks samples tested in compression which forms small columns (Wang and Shrive 1995).

### Case II: Failure of Single Crystals Governed by Normal and Shear Stress

When the macroscopic stress state is uniaxial ($\Sigma_1=\Sigma_2=0$, $\Sigma_3$), microstructure and loading are invariant with respect to rotation around $\epsilon_3$. Thus, the orientation of the critical crystal (which maximizes the local criterion) only depends on $\theta$. For a tension ($\theta_t^{cr}$) or compression ($\theta_c^{cr}$) test, the critical orientation is characterized by

$$ \theta_t^{cr} = \frac{\arctan A}{2} \quad \text{and} \quad \theta_c^{cr} = \frac{\pi - \arctan A}{2} \qquad (40) $$

referring to the results presented in Table 7. These angles depend on $\beta$ through $A=a_2/\beta$. At high porosities, the critical crystal is found to be near the axis of tension or nearly normal to the axis of compression (see Fig. 27).

It is interesting to investigate the normal $\sigma_N$ and shear $\sigma_T$ stresses in the critical crystal just when the macroscopic elastic limit is reached. When the macroscopic stress state is uniaxial of the form $\Sigma_1=\Sigma_2=0$ and $\Sigma_3$, Eqs. (23) and (30) become

<figure><img src="image_29.png" /><figcaption>**Fig. 26.** Critical crystals for tension and compression tests</figcaption></figure>

$$ \sigma_N = \Sigma_3(a_0 + a_1 \cos^2 \theta) \quad \text{and} \quad \sigma_T = |\Sigma_3 \cos \theta| a_1 a_2 \sin \theta \quad (41) $$

Substituting $\theta=\theta_t^{cr}$ or $\theta_c^{cr}$ into these expressions, the curves presented in Fig. 28 can be plotted. At high porosities ($\phi > 0.6$), the critical crystal is nearly under pure tension $\sigma_T^{cr} \approx 0$. Thus, the first case (studied in the previous subsubsection, local criterion $\sigma_N < \sigma_{cr}$), is nearly retrieved.

Note that in tension, when $\phi \to 0$ (in the asymptotic case of a material without voids, where the individual crystals cannot be distinguished), one could believe that the local criterion would be written as $\sigma_N < \sigma_{cr}$ since the shear component of the stress vector $\sigma_T$ tends toward 0 on a surface perpendicular to $N=\epsilon_3$. In fact, it is a misbelief, since the horizontal facet is not the most defavorable; from Fig. 27, it is clear that the critical facet is characterized by $\theta^{cr} \approx 20^\circ$.

## Conclusions

The classical self-consistent scheme has been successfully extended in order to deal with nonspherical crystals. Moreover, the good agreement with both numerical and experimental data (including various formulations of gypsum leading to various crystals aspect ratios), supports the fact that this extended self-consistent scheme is able to capture the influence of the shape of the particles on the effective stiffness under the assumption of an isotropic orientation distribution. Traditionally, the gypsum material is made up of elongated crystals. A simplified analysis of the effective strength has been conducted in the limit case of needle-like particles (infinite aspect ratio). The domain of elastic stress states proves to be a cone with polygonal basis, the shape of which depends on the description of the microscopic strength domain. More precisely, an equilateral triangle is obtained when

<figure><img src="image_30.png" /><figcaption>**Fig. 27.** Orientation of the critical crystal in tension and in compression ($v_s$=0.34, $\sigma_{cr}$=17 MPa, and $\beta$=1.3)</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 28.** Stresses $\sigma_N$ and $\sigma_T$ in the critical crystal at the elastic limit in tension and in compression ($v_s=0.34$, $\sigma_{cr}=17$ MPa, and $\beta=1.3$)</figcaption></figure>

the microscopic strength is controlled by the tensile axial stress in the needle. In turn, a hexagon is obtained when the local criterion involves both the axial and the tangential components of the stress acting on the cross section of the needle.

Despite the fact that the most likely local failure mechanism (failure in the interfaces) has been replaced by a bulk failure mechanism (failure in the crystals), encouraging results have been obtained. Indeed, the predictions are consistent with experimental data both in compression and in tension provided that the local criterion takes into account the tangential stress component in the needle. This strongly supports the idea that the microscopic strength criterion is both cohesive and frictional in nature.

The models proposed in this paper have not taken into account the presence of entrapped air voids (bubbles) in gypsum. Gypsum always contains a more or less large amount of bubbles, even if special care has been taken during the preparation of the samples. These defects obviously affect both elasticity and strength (Vekinis et al. 1993) of the material. A possible improvement is to incorporate these air voids, for example by resorting to a second homogenization step. It is quite straightforward as far as elasticity is concerned. The strength model would consider bubbles as defects around which stress concentration builds up.

Micromechanics approaches as developed herein are also relevant for other materials made up of interlocking nonspherical particles, such as porous hydroxapatite biomaterials (Fritsch et al. 2006, 2007b, 2009).

## Notation

The following symbols are used in this paper:

$\underline{A}$ = first-order tensor;

$\underline{A}$ = second-order tensor;

$\underline{A}$ = fourth-order tensor;

$A_d$ = magnitude of the deviatoric part of $\underline{A}$, $A_d = \sqrt{\underline{A}^{\text{dev}}:\underline{A}^{\text{dev}}/2}$;

$\underline{A}^{\text{dev}}$ = deviator of $\underline{A}$, $\underline{A}^{\text{dev}}=\underline{A}-1/3(\text{tr } \underline{A})\underline{1}$;

$A_m$ = mean part of $\underline{A}$, $A_m=(\text{tr } \underline{A})/3$;

$a_0, a_1, a_2$ = coefficients appearing in the concentration relations providing $\sigma_N$ and $\sigma_T$;

$(\underline{a}, \underline{b}, \underline{\epsilon})$ = elementary crystallographic lattice of De Jong and Bouman for gypsum;

$\langle a \rangle_j$ = average of field $a$ over the domain occupied by phase $j$;

$C_j$ = stiffness tensor of phase $j$;

$C^{\text{SCS}}$ = self-consistent estimate of the effective stiffness;

$\underline{E}$ = strain tensor applied at the boundary of the REV;

$\underline{E}^{\text{eff}}$ = effective Young's modulus;

$E_{33}^{\text{lim}}$ = macroscopic strain at failure in the direction of loading;

$E_j$ = Young's modulus of phase $j$;

$\underline{E}_0$ = reference strain;

($\underline{\epsilon}_r, \underline{\epsilon}_0, \underline{\epsilon}_\phi$) = spherical frame;

($\underline{\epsilon}_1, \underline{\epsilon}_2, \underline{\epsilon}_3$) = reference frame;

$f_j$ = volume fraction of phase $j$ in the REV;

$g_j$ = shear modulus of phase $j$;

$h$ = magnitude of a spherical macroscopic stress state ($\Sigma=h\underline{1}$);

$h^{\text{max}}$ = maximum magnitude of an admissible macroscopic stress state;

$\underline{I}$ = identity fourth-order tensor for tensors presenting the minor symmetry;

$\underline{J}$ = tensor extracting the spherical part of a symmetrical second-order tensor ($\underline{J}=1/3\underline{1}\otimes\underline{1}$);

$\underline{K}$ = tensor extracting the deviatoric part of a symmetrical second-order tensor ($\underline{K}=1-\underline{J}$);

$k_j$ = bulk modulus of phase $j$;

$\underline{N}$ = axis of a given crystal;

$\underline{P}^j_{\text{shape}}$ = Hill tensor of shape in phase $j$ if phase $j$ is isotropic $C_j=3k_j\underline{J}+2g_j\underline{K}$;

$r, \theta, \phi$ = spherical coordinates;

$r_s$ = aspect ratio;

$S^j_{\text{shape}}$ = Eshelby tensor of shape in phase $j$;

$\underline{T}$ = vector perpendicular to $\underline{N}$;

$\underline{T}^{-1}$ = inverse of a second or fourth-order tensor;

$\text{tr } \underline{A}$ = trace of $\underline{A}$;

($\underline{u}_1, \underline{u}_2, \underline{u}_3$) = crystal frame;

$\underline{v}^{\text{eff}}$ = effective Poisson's ratio;

$v_j$ = Poisson's ratio of phase $j$ with the classical relationships: $k_j=E_j/3(1-2v_j)$, $g_j=E_j/2(1+v_j)$, $E_j=9k_jg_j/3k_j+g_j$, and $v_j=3k_j-2g_j/6k_j+2g_j$;

$\underline{\alpha}$ = friction coefficient of a Drucker-Prager criterion;

$\underline{\epsilon}$ = strain tensor;

$\underline{\mathcal{E}}$ = macroscopic elastic domain;

$\underline{\mathcal{E}}_0$ = macroscopic elastic domain at a vanishing porosity;

$\underline{\theta}^{\text{cr}}_t$ = critical crystal orientation for a macroscopic tensile test;

$\underline{\theta}^{\text{cr}}_c$ = critical crystal orientation for a macroscopic compressive test;

$\underline{\theta}, \underline{\phi}, \underline{\psi}$ = Euler angles;

$\underline{\xi}$ = displacement vector;

$\underline{\sigma}$ = stress tensor;

$\underline{\sigma}_{cr}, \underline{\beta}$ = parameters characterizing the strength at the crystals scale;

$\underline{\sigma}_N^{\text{cr}}$ = axial part of the stress vector on the cross section of the critical crystal at failure;

$\underline{\sigma}_T^{\text{cr}}$ = tangential part of the stress vector on the cross section of the critical crystal at failure.

$\underline{\sigma}^N$ = average stress tensor in the crystal of axis $\underline{N}$;

$\underline{\sigma}^{N\cdot\underline{N}}$ = stress vector acting on the cross section of the crystal of axis $\underline{N}$;

$\underline{\sigma}_N$ = axial part of $\underline{\sigma}^N\cdot\underline{N}$;

$\underline{\sigma}_T$ = tangential part of $\underline{\sigma}^N\cdot\underline{N}$;

$\underline{\sigma}_T$ = magnitude of $\underline{\sigma}_T$;

$\underline{\Sigma}$ = average stress tensor in the REV;

$\underline{\Sigma}^{\text{lim}}_t$ = macroscopic tensile elastic limit ($\underline{\Sigma}^{\text{lim}}_t>0$);

$\underline{\Sigma}^{\text{lim}}_c$ = macroscopic compressive elastic limit ($\underline{\Sigma}^{\text{lim}}_c>0$);
$\Sigma_1, \Sigma_2, \Sigma_3$ = macroscopic principal stresses;

$\Sigma_I, \Sigma_{II}, \Sigma_{III}$ = reverse ordered macroscopic principal stresses;

$\varphi$ = porosity;

$\varphi^c$ = critical porosity above which the effective stiffness vanishes;

$\Omega$ = domain occupied by the REV (representative elementary volume);

$\Omega_j$ = domain occupied by phase *j* over the REV; and

$\mathbf{1}$ = identity second-order tensor.

## References

Ali, M. A., and Singh, B. (1975). "The effect of porosity on the properties of glass fibre-reinforced gypsum plaster." *J. Mater. Sci.*, 10(11), 1920–1928.

Alsadi, S., Combe, E. C., and Cheng, Y.-S. (1996). "Properties of gypsum with the addition of gum Arabic and calcium hydroxide." *J. Prosthet. Dent.*, 76(5), 530–534.

Cherkaev, A. E., Lurie, K. A., and Milton, G. W. (1992). "Invariant properties of the stress in plane elasticity and equivalence classes of composites." *Proc. R. Soc. London, Ser. A*, 438(1904), 519–529.

Colak, A. (2006). "Physical and mechanical properties of polymer-plaster composites." *Mater. Lett.*, 60(16), 1977–1982.

Dany, C. (1995). "Étude de l'influence de la microstructure sur les propriétés mécaniques du plâtre sec." MA thesis, Génie Physique des Matériaux, INSA, Lyon, France.

Day, A. R., Snyder, K. A., Garboczi, E. J., and Thorpe, M. F. (1992). "The elastic moduli of a sheet containing circular holes." *J. Mech. Phys. Solids*, 40(5), 1031–1051.

De Jong, W., and Bouman, J. (1938). *Z. Kristallogr.*, 100, 275.

Desrues, J. (2002). "Limitations du choix de l'angle de frottement pour le critère de plasticité de drucker-prager." *Revue Française de Génie Civil*, 6(5), 853–862.

Dormieux, L., Sanahuja, J., and Maalej, Y. (2007). "Résistance d'un polycristal avec interfaces intergranulaires imparfaites." *C. R. Mec.*, 335(1), 25–31.

Eshelby, J. D. (1957). "The determination of the elastic field of an ellipsoidal inclusion, and related problems." *Proc. R. Soc. London, Ser. A*, 241(1226), 376–396.

Fairhurst, C. W. (1960). "Compressive properties of dental gypsum." *J. Dent. Res.*, 39(4), 812–824.

Fritsch, A., Dormieux, L., and Hellmich, C. (2006). "Porous polycrystals built up by uniformly and axisymmetrically oriented needles: Homogenization of elastic properties." *C. R. Mec.*, 334(3), 151–157.

Fritsch, A., Dormieux, L., Hellmich, C., and Sanahuja, J. (2007a). "Micromechanics of crystal interfaces in polycrystalline solid phases of porous media: Fundamentals and application to strength of hydroxyapatite biomaterials." *J. Mater. Sci.*, 42(21), 8824–8837.

Fritsch, A., Dormieux, L., Hellmich, C., and Sanahuja, J. (2007b). "Micromechanics of hydroxyapatite-based biomaterials and tissue engineering scaffolds." *Tissue engineering using ceramics and polymers*, Woodhead Publishing Limited, 529–565.

Fritsch, A., Dormieux, L., Hellmich, C., and Sanahuja, J. (2009). "Mechanical behaviour of hydroxyapatite biomaterials: An experimentally validated micromechanical model for elasticity and strength." *J. Biomed. Mater. Res. Part A*, 88A(1), 149–161.

Garboczi, E. J., Snyder, K. A., Douglas, J. F., and Thorpe, M. F. (1995). "Geometrical percolation threshold of overlapping ellipsoids." *Phys. Rev. E*, 52(1), 819–828.

Gmouh, A., Eve, S., Samdi, A., Moussa, R., Hamel, J., and Gomina, M. (2003). "Changes in plaster microstructure by pre-stressing or by adding gypsum grains: Microstructural and mechanical investigations." *Mater. Sci. Eng., A*, 352, 325–332.

Haussühl, S. (1965). "Elastische und thermoelastiche eigenschaften von CaSO₄, 2H₂O (gips)." *Z. Kristallogr.*, 122, 311–314.

Khalil, A. A. (1972). "Gypsum plasters. II: Constitution and properties of the prepared β-hemihydrates." *J. Appl. Chem. Biotechnol.*, 22, 703–709.

Klima, W. F. (1983). "The kinetics of crystal growth of calcium sulfate dihydrate." Ph.D. thesis, State Univ. of New York at Buffalo, New York.

Kröner, E. (1977). "Bounds for effective elastic moduli of disordered materials." *J. Mech. Phys. Solids*, 25, 137–155.

Meille, S. (2001). "Étude du comportement mécanique du plâtre pris en relation avec sa microstructure." Ph.D. thesis, INSA, Lyon, France.

Meille, S., and Garboczi, E. J. (2001). "Linear elastic properties of 2D and 3D models of porous materials made from elongated objects." *Modell. Simul. Mater. Sci. Eng.*, 9(5), 371–390.

Péronnet, F. (1994). "Caractérisation du comportement mécanique des plâtres." *Rapport de Stage*, Génie Physique des Matériaux, INSA, Lyon, France.

Phani, K. K. (1986). "Young's modulus-porosity relation in gypsum systems." *Am. Ceram. Soc. Bull.*, 65(12), 1584–1586.

Sanad, M. E., Combe, E. C., and Grant, A. A. (1982). "The use of additives to improve the mechanical properties of gypsum products." *J. Dent. Res.*, 61(6), 808–810.

Sanahuja, J., Dormieux, L., and Chanvillard, G. (2008). "A reply to the discussion 'Does C-S-H particle shape matter?' by F.-J. Ulm and H. M. Jennings of the paper 'Modelling elasticity of a hydrating cement paste,' CCR 37 (2007)." *Cem. Concr. Res.*, 38(8–9), 1130–1134.

Schinkel, H. (1993). "Creep of gypsum under flexural load." *Technical Rep.*, Bouwcentrum Advies, Rotterdam, The Netherlands.

Simon, B., and Bienfait, M. (1965). "Structure et mécanisme de croissance du gypse." *Acta Crystallogr.*, 19(5), 750–756.

Tazawa, E. (1998). "Effect of self stress on flexural strength of gypsum-polymer composites." *Adv. Cem. Based Mater.*, 7(1), 1–7.

Vekinis, G., Ashby, M. F., and Beaumont, P. W. R. (1993). "Plaster of Paris as a model material for brittle porous solids." *J. Mater. Sci.*, 28(12), 3221–3227.

Wang, E. Z., and Shrive, N. G. (1995). "Brittle failure in compression: Mechanisms, models and criteria." *Eng. Fract. Mech.*, 52(6), 1107–1126.