Stochastic finite element analysis of the diametral compression test of powder compacts and porous materials

This paper presents a stochastic finite element approach for modeling the mechanical behavior of powder compacts and porous materials under diametral compression test conditions. The main goal is assessing the validity of the diametral compression test as an indirect technique to estimate tensile strengths of brittle or quasi-brittle materials exhibiting porosity heterogeneity. Thus, the study seeks to predict the influence of porosity randomness on stress distributions and the spatial location of the highest tensile stress on thin disc-shaped specimens. The proposed formulation uses a stochastic framework that couples a random spatial field to the finite element analysis to include non-deterministic features. Two case studies consider comparable targets for the mean porosity but different coefficients of variations. For each case study, a total of 1000 realizations are conducted under identical loading and boundary conditions. The predicted stress distributions are compared to the ones from homogenous closed-form solutions from the literature. Then, the expected magnitude and location of the maximum tensile stress are evaluated by statistical means. Findings from the stochastic model show that porosity randomness induces stress concentration around less dense volumes and location deviation of the maximum tensile stress from the center of the discs. Likewise, porosity heterogeneity could affect the accuracy of experimental diametral compression tests even for small variance cases; and so, the reliability of the mechanical properties derived from models based exclusively on the classic assumption of material homogeneity.


INTRODUCTION
Porous materials have found ordinary and extensive uses in diverse fields. From biomedical to aerospace applications, porous materials prove successful when biomimetic properties are of particular importance or when high mechanical strengths are required with simultaneous weight reduction [1][2][3]. Porous surfaces ensure a proper bond between the prosthesis and the host bone in implants made of biocompatible materials such as calcium phosphate ceramics [4][5][6]. Porous materials have also found widespread uses in the aerospace industry, going from simple filtration systems to applications such as gas sensors and separators, thermal isolators, or heat exchangers [7][8][9]. For this reason, a precise mechanical characterization of porous materials is of significance to ensure the highest possible performance of the intended engineering application.
Porous materials are naturally found or industrially produced by partially densifying powder mixes using compaction techniques [10]. Properties of porous materials are not only influenced by the compaction pressures and temperatures but also by the physicochemical properties of the raw powders [11]. Among other features, particle size, particle shape, particle size distribution, roughness, and surface impurities govern the kinetics of consolidation of powder materials and so final mechanical properties. For instance, powders with low flowability and strongly agglomerated can lead to materials with broad pore size distributions or even multimodal distributions [12].
Compressive and tensile strengths, stiffness, and fracture toughness are among the most important parameters to characterize porous materials [13]. Thus, different experimental approaches are available to assess these properties. In recent decades, the diametral compression test has become an accepted experimental technique to gauge indirectly tensile strengths of brittle or quasi-brittle materials. The diametral compression test, a.k.a. Brazilian test, was initially designed to estimate tensile strengths of concretes [14].
Like the experimental setups schemed in Fig. 1, the Brazilian test consists of a disc-shaped specimen subjected to diametral forces applied by two flat plates or curved jaws. [15]. In a thin disc, loading conditions induce a plane stress state with nearly constant tensile stress along the loaded diameter and compressive stress along with the perpendicular one. 4 Numerous studies are available on the Brazilian test to evaluate properties of rocks [16] along with the ones of other homogeneous and heterogeneous media [17][18][19][20][21][22]. Among those works, different modeling approaches have been successfully implemented in commercial or inhouse finite element codes; however, to the best authors' knowledge, none of them has considered modeling diametral compression of powder compacts from a stochastic viewpoint.
Always there has been an accumulating interest in the Brazilian test to assess the tensile strengths of diverse heterogeneous materials. For instance, Yue et al. [23] used digital image techniques to carry out finite element simulations of non-homogeneous geomaterials under diametral compression test conditions, finding that stress distributions and initial crack location are significantly affected by the internal structure of those materials. Jonsén et al. [24] employed a Brazilian test setup to study the mechanical strength, fracture instability onset, and fracture energy of green iron powder compacts of different densities. They suggested a procedure to estimate tensile strengths by analyzing the linear to non-linear transition from load-displacement curves. From the experimental outcomes, Jonsén et al. [24] discussed the large spread observed in the correlation between green strengths and densities.
Jonsén and Haggblad [25] studied the tensile fracture of metal powder compacts by considering a constitutive model formulated in terms of non-linear fracture mechanics and coupled to a perfect plastic failure envelope with a hardening cap. Then, Korachkin et al. [26] Fig. 1 Experimental setups using (a) standard flat loading platens and (b) curved loading jaws employed an experimental diametral compressive setup to assess elastic moduli and tensile strengths of alumina and Distaloy AE compacts. They compared the attained resistances and failure modes to the ones obtained from three-point bending tests, suggesting that result differences are because heterogeneity introduces local weaknesses that prompt failure.
Seeber et al. [27] investigated the suitability of the Brazilian test alongside the 3-point bending and compression tests as characterization techniques of the mechanical response of highly porous alumina with diverse microstructures. Cui et al. [28] presented a micromechanical model to predict properties of brittle materials with different porosities and pore sizes in diametral compression. In that work, they evaluated the microscale mechanical properties and correlated the impact of porosity and pore size on the fracture strength, Weibull modulus, and degradation processes of the specimens.
Uncertainties in powder material design are unavoidable. Parameters such as powder mix properties or compaction processes involve significant levels of uncertainty.
Consequently, some numerical simulations are non-representative studies when based exclusively on averaged properties and limit isolated phenomena. Thus, this paper presents a computational approach in the behavior of disc-shaped specimens of porous materials in diametral compression to gain insight into inhomogeneity's significance on test results. The parameter of interest is porosity, which defines the pore volume to total volume fraction. In a probabilistic framework, a spatial stochastic field describes material heterogeneity, and a finite element model estimates displacement and stress fields. The spatial random field employs a discretization method that guarantees a truthful representation of the material together with a computational cost reduction.
The paper organizes as follows. The "Materials and methods" section highlights the most relevant equations of the proposed model. First, the constitutive formulation is presented along with an introduction to the approach used to approximate porosity random fields. Next, the finite element stiffness matrix is expressed in terms of the uncertain parameter. In the last part of the paper, the "Numerical example and discussion" section presents two case studies and how the model is used to simulate the mechanical behavior of disc-shaped specimens with non-homogeneous porosity distributions. Results are compared with the ones obtained from analytical and numerical models under the assumption of material homogeneity, and finally, conclusions are drawn.

MATERIALS AND METHODS
This section deals with the procedures that translate the material behavior into the form of constitutive relationships and mathematical terms. In particular, attention is on the modeling strategy of porous disc-shaped specimens in diametral compression and on the uncertainty propagation through the outcomes. Considering porosity as a time-independent parameter, it turns out that pore collapse does not affect deformations and that the deformation history does not govern stress states. Consequently, the material behaves as isotropic linear elastic with no local viscous-plastic effects. Thus, only two independent elastic properties: the elastic modulus and Poisson's ratio, describe the linear behavior. In the formulation, the material is thought of being a porous skeleton containing open and closed pores, which distribute randomly in shape and size. Likewise, the total material volume is the volume sum of substance and pores. Thus, using continuum mechanics, the so-called effective local elastic modulus is porosity-dependent [29,30]. In contrast, Poisson's ratio is an intrinsic property that does not vary significantly within the mesoscale.

Constitutive formulation
The properties of non-homogeneous materials depend on the links between the different length scales and the complexity of each substructure [31][32][33]. For instance, heterogeneities induce randomness and affect both the local and the overall mechanical behavior of materials.
However, to cast the problem on a single scale, the current formulation does not involve all small-scale features on the macroscopic response. Hence, field equations describing the elastic behavior are based upon continuum mechanics and given by the motion equation, the kinematic relationships, and the constitutive law. For the case studies, specimens are subjected to in-plane loads, and the thickness is considered much smaller than the diameter, as illustrated in Fig. 2a. In such circumstances, it is quite convenient to treat the problem as a plane stress case defined on the domain D with boundaries ∂D , as shown in Fig. 2b.
Accordingly, (i) the xy − plane coincides with the middle plane of the disc, (ii) the contact between the disc and the loading jaw induces load distributions, and (iii) all z − stress components vanish. The residual boundaries u ∂ ∂ D D are stress-free.
Thus, from the above hypotheses, the functional relationship between stress σ and strain ε components for isotropic materials is given as where E and ν are respectively the elastic modulus and Poisson's ratio. The out-of-plane strain component is written as ε Material heterogeneities introduce local weaknesses that trigger stress concentrations and strength drops. As pointed out, the material is non-homogeneous due to the inherent porosity variability. Thus, the constitutive model captures the underlying variation via a spatial random field function. This function contains essential local information from the mesoscale and provides a realistic approximation of the material through the effective property. Therefore, upscale simulations are performed, introducing the material stiffness as

Description of the random porosity field
In this work, porosity variability is described by using spatial statistics, where the local quantity is defined by a probability distribution and correlated with values at adjacent locations. Thus, let x be a material point within a circular domain R is the specimen's radius. Then, let us define ( ) Θ x as a real-valued random function representing porosity at each point. From the theory of probabilities [34], ( , , ) Ω F P represents the probability space components, i.e., the outcome sample space , Ω the event collection F over the subsets of , Ω and the probability measure P of each member of . F.
Thus, let us think of the random spatial field as a mapping defined from the product space where ω is a realization sampled from the probability space Ω according to . P Realizations are independent of each other and equally likely. Figure 3 depicts the proposed probability space with a selection of porosity realizations. Accordingly, ( , ) Θ • ω refers to the random variable resulting from realization ω, while ( , ) Θ • x represents the random variable for a fixed . x Material variability can be quantified by a probability distribution when sufficient statistical information is available. In what follows, let us suppose that the probability measure is known and that statistical parameters are estimated. In the case of random fields, the classic definition of mean and variance can be extended directly [35]; therefore, the and the deviation from this expected value, i.e., the variance, as is the mathematical expectation operator.
Likewise, for each ( , ), ′ x x the covariance function is expressed as For the case that concerns us, porosity sample fields are regarded as uncorrelated random variables with zero mean and finite variance, arranged along with a spatial space in For the sake of simplicity, porosity samples are supposed to be independent and identically distributed (i.i.d.). Accordingly, for each pair ( , ), fields are assumed to be statistically Fig. 3 A selection of realizations in the porosity space independent following a Gaussian probability density function. Therefore, the covariance matrix K is described by a white noise random field proportional to the Dirac-delta function δ Overcoming the issue that the white noise variance process is infinite, a grid-based approach is used to approximate the corresponding random field. A concise review on white noise modeling can found in [36].

Finite element formulation
The mechanical performance of the porous compacts is approximated numerically as an elastic continuum. Therefore, the problem is described by the linear constitutive model is the inverse of the stochastic stiffness matrix. For the sake of simplicity, the force vector F is supposed to be deterministic. From Eq. (9), it seems that uncertainty propagates through the response variability of the system ( , ) Θ U x [37,38].
A Monte Carlo approach uses a sampling method to generate a set of conditional realizations of the field Θ over the probability space. Next, realizations are used as the inputs of the finite element analyses to estimate nodal displacement's variability, i.e., ( , ).
ω U x Finally, the stochastic response defines both the strain and stress fields through the kinematic and constitutive equations, respectively.

Random field discretization
Within the stochastic finite element framework, it requires converting porosity into a discrete spatial distribution. Consequently, the random field ( , ) characterizes porosity values at discrete points x on .

( , )
i j x y + + depicted in Fig. 5, the simplest surface between the four vertices is given as where 1 c and 2 c represent bilinear weights, given as

Analytical model
The indirect tensile test of homogeneous and isotropic elastic materials has been analytically treated for some authors. Frocht [40] and Timoshenko and Goodier [41] proposed a solution for the stress distributions in discs subjected to a pair of concentrated compressive loads applied across the diameter. According to their model, the , x− , y− and xy − stress components at any point within the discs are described by the following equations where P is the concentrated load, R and t are the radius and thickness, and and are the coordinates of the point of interest. Some other closed-form solutions have been proposed in the literature under different boundary conditions [42][43][44]. It is seen from Eq. (12) that the analytical solution is independent of the elastic properties of the material and is exclusively a function of the force magnitude and disc's dimensions. Likewise, it is obvious from Eq.
(12) that the maximum tensile stress takes place at the center of the disc, where the stress state is given as The compression test concerns itself with the maximum tensile stress reaching a critical value known as indirect tensile strength. It is upon the assumption that the initial crack appears roughly near the center of the disc and propagates along the loading direction. From Eq. (13), the stress state at the center point corresponds to a biaxial one where the compressive component is three times larger than the tensile one and where the shear stress vanishes. Therefore, it is expected that due to opposite sign stresses, indirect tensile strengths determined from the Brazilian test could be slightly higher than the ones obtained from direct tensile tests [45,46].

Deterministic numerical analysis
The numerical setup validity is first assessed through the predicted stress distributions and maximum stress values, neglecting uncertainties. Thus, let us analyze the test under the assumption that porous material properties correspond to a case restricted to average values and that boundary conditions match the test operating conditions. This deterministic case qualitatively compares the numerical stress distributions with the analytical ones proposed by Frocht [40] and normalizes the results obtained in the stochastic analysis described in the last part of this section. An unstructured mesh with Lagrange triangular elements discretizes the problem domain, as shown in Fig. 4a. Respectively, a total of 1016 and 140 elements are used for the specimen and jaw meshes. Figure

Stochastic analysis
Analytical stress distributions do not suit heterogeneous materials. The closed-form solution in the "Analytical model" section is independent of the elastic properties; though, it is expected stress variability to intensify with heterogeneity. Hence, different realizations are generated to study the relations between porosity and behavior at the macroscopic scale. In this section, the formulation introduces porosity according to its spatial fluctuation. The analysis considers two case studies to investigate solution dependency on porosity. Each analysis with 1000 realizations represents correspondingly a somewhat homogeneous material and a more heterogeneous one. Following Wachtman [47] and Ashkin [48], the relationship between porosity and the elastic modulus is treated as linear and quantitatively deterministic; therefore, where 0 E is the elastic modulus of the fully dense material and θ is the local porosity. Other empirical relationships can be found in [49][50][51][52].
The numerical simulations use the two-dimensional fields described in the "Description of the random porosity field" section. The mean θ µ and coefficient of variation ω θ x is a vector of uncorrelated random porosity values such as [53] ( , ) where ( , ) ω ς x is a zero-mean vector of uncorrelated Gaussian random variables; thus, represents porosity fluctuation about the mean value θ µ . Figure 8   highest tensile stress is equal to 13.9 MPa, and its location is on the loading plane; even so, it deviates from the horizontal diameter. For the second case, the maximum tensile stress is equal to 14.7 MPa, and its location departs from the loading plane and more significantly from the horizontal diameter.
Scatter plots in Fig. 12

CONCLUSIONS
This work presented a stochastic assessment of the impact of material heterogeneity on the mechanical performance of the diametral compression test of powder compacts and porous materials. In the proposed stochastic formulation, a random spatial field described porosity