Predicting the mechanical properties of brittle porous materials with various porosity and pore sizes

Author’s Accepted Manuscript Predicting the Mechanical Properties of Brittle Porous Materials with Various Porosity and Pore Sizes Zhiwei Cui, Yongmin Huang, Honglai Liu www.elsevier.com/locate/jmbbm

PII: DOI: Reference:

S1751-6161(17)30072-3 http://dx.doi.org/10.1016/j.jmbbm.2017.02.014 JMBBM2235

To appear in: Journal of the Mechanical Behavior of Biomedical Materials Received date: 22 November 2016 Revised date: 11 February 2017 Accepted date: 12 February 2017 Cite this article as: Zhiwei Cui, Yongmin Huang and Honglai Liu, Predicting the Mechanical Properties of Brittle Porous Materials with Various Porosity and Pore Sizes, Journal of the Mechanical Behavior of Biomedical Materials, http://dx.doi.org/10.1016/j.jmbbm.2017.02.014 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

-1-

Predicting the Mechanical Properties of Brittle Porous Materials with Various Porosity and Pore Sizes Zhiwei Cui, Yongmin Huang*, Honglai Liu Key Laboratory of Specially Functional Polymeric Materials and Related Technology, School of Chemistry and Molecular Engineering, East China University of Science and Technology, Shanghai 200237, China.

*Corresponding author at: Key Laboratory of Specially Functional Polymeric Materials and Related Technology, School of Chemistry and Molecular Engineering, East China University of Science and Technology, Shanghai 200237, China. E-mail address: [email protected]

-2-

ABSTRACT In this work, a micromechanical study using the lattice spring model (LSM) was performed to predict the mechanical properties of BPMs by simulation of the Brazilian test. Stress-strain curve and Weibull plot were analyzed for the determination of fracture strength and Weibull modulus. The presented model composed of linear elastic elements is capable of reproducing the non-linear behavior of BPMs resulting from the damage accumulation and provides consistent results which are in agreement with experimental measurements. Besides, it is also found that porosity shows significant impact on fracture strength while pore size dominates the Weibull modulus, which enables us to establish how choices made in the microstructure to meet the demand of brittle porous materials functioning in various operating conditions.

-3Keywords: Lattice spring model; Porosity; Brittle fracture; Mechanical properties; Weibull modulus

Nomenclature 𝐸𝐴 , 𝐸𝐵

Young’s modulus of the solid/void phase

𝜈

Poisson’s ratio

𝑘

the central force constant of the spring

𝑐

the non-central force constant of the spring

𝑀𝑖𝑗 , 𝑀[𝑖𝑗𝑘]

the matrix containing force constants for the spring connecting nodes 𝑖 and 𝑗, or that in the direction of [𝑖, 𝑗, 𝑘]

𝑅

the rotation matrix in three dimensions

𝐹𝑖

the force acting on the 𝑖th node

𝑢𝑖 , 𝑟𝑖

the displacement of the 𝑖th node from its original position

𝑎𝑖

the acceleration of the 𝑖th node

𝑚𝑖

the mass of the 𝑖th node

𝑣𝑖

the velocity of the 𝑖th node

Θ

the viscous damping constant

𝑡

time

𝜀𝑥

normal strain along the x-axis

ℎ

the initial distance between adjacent nodes

𝐹𝑠

the force acting on any surface 𝑠 of a cubic cell

𝐴

the area of a surface 𝑠

-4𝑠 𝑝𝑘𝑙

the unit vector either normal (𝑘𝑙 = 11,22,33) or parallel (𝑘𝑙 = 12,23,13) to a surface 𝑠

𝜎𝑘𝑙

the stress acting on a cubic cell

𝐸𝑖

the local elastic energy stored in the 𝑖th node

𝑝𝑖 (𝑡)

the rate of failure associated with node 𝑖 at time 𝑡

𝑊𝑖

the minimum value of strain at which element breakage can happen

𝐵𝑖

the arbitrary scaling parameter

𝛽

a modulus determining the relationship between the rate of failure and elastic energy field

𝑃𝑖 (𝑡)

the probability of failure associated with node 𝑖 at time 𝑡

𝑐𝑖 (𝑡)

the cumulative probability of failure associated with node 𝑖 at time 𝑡

𝑃𝑓

the probability of failure defined in Weibull statistics

𝜎𝑢

the threshold stress parameter

𝜎𝜃

the characteristic strength defined in Weibull statistics

𝑚

the Weibull modulus

𝑁

the sample size

𝑟

the radius of the cylindrical shaped specimen

𝑤

the height of the cylindrical shaped specimen

𝜔0

the contact angle between the jaws and the specimen

𝑙

the contact length between the jaws and the specimen

1. INTRODUCTION

-5Brittle porous materials such as ceramics and bioactive glasses possess low weight, high porosity, excellent chemical and thermal stability and even bioactivity.1,2 These properties have led to extensive applications including catalyst supports,3 membrane reactor,4 solid oxide fuel cell insulators and scaffolds used in bone regeneration.5-7 However, implementation of ceramic membrane as well as porous bioactive glass scaffold has been greatly limited mainly because of their low fracture toughness.8-10 Besides, physical breakage of catalyst pellets leads to the formation of fragments and fines, causes misdistribution of fluid flow and increases the pressure drop of chemical reactor,11 in spite of the relative high fracture strength of catalyst supports. Therefore, an understanding of the failure behavior of brittle porous materials is vital in developing materials for various applications. The prediction of the fracture strength has been one of the most important research areas in the field of porous materials. Generally, fracture strength decreases with the porosity,12 which is proved by both experimental and numerical researches.13-16 A number of expressions can be used to describe the porosity-strength behavior of porous ceramics, and one of the simplest form was proposed based on the minimum solid area concept.17 The porosity-strength dependence is approximated by an exponential function, 𝜎 = 𝜎0 exp(−𝑏𝑝), where 𝜎0 is zero-porosity strength, 𝜎 is the strength at pore volume fraction 𝑝, and the constant 𝑏 is related to the pore characteristics.17,18 Specifically, zero-porosity strengths are supposed to be identical for systems with same level of porosity, regardless of the pore size. Minimum solid area concept is shown to be rational for elastic properties as well as tensile, compressive and flexural strength.19,20 Nonetheless, obvious deviation from theoretical expectation was observed. Liu employed three types of PVB particles of sizes 0.093mm,

-60.188mm and 0.42mm as pore-forming material to fabricate HAP ceramics of various pore size for the investigation of porosity-strength behavior.18 The zero-porosity strengths for the ceramics were 5920Mpa, 2580Mpa and 1900Mpa for 0.093mm, 0.188mm and 0.42mm particles, respectively, which is in contradiction to theoretical prediction. Zhang and Malzbender assessed the fracture strength of nano- and micro-porous alumina by the use of ring-on-ring bending tests.21 It turns out that in spite of equal porosity, elastic moduli and fracture stress values revealed significant differences. Therefore, it seems important to characterize the effect of porosity together with the macropore size on the compressive strength for a better understanding of the porosity-strength behavior. On the other hand, the fracture strength of porous materials, especially brittle porous materials, is found scattered over a wide range of values.22,23 Therefore, the Weibull modulus, which reflects the scattering extent of fracture strength, is adopted as an important parameter for the measurement of mechanical reliability.

24,25

A decreasing Weibull modulus with

increasing porosity has been reported for niobium-and tantalum-doped zirconia,26 titania,27 alumina and hydroxyapatite.13,24,28-30 Also, Keleş et al. observed that Weibull modulus decreased more than threefold for a change in porosity from 1 to 2 vol.% by simulating the uniaxial tensile test of porous systems containing non-overlapping circular pores with identical size of 80μm on the basis of 2D finite element method.15 Meanwhile, the Weibull modulus did not show a significant change with increasing porosity for NiO-YSZ (yttria-stabilized zirconia) and MgAl2O4.31,32 However, in the works of Frandsen et al. and Fan et al., an improvement in Weibull modulus was detected during a variation of porosity from 10 ~ 15 vol.% and 55 ~ 62 vol.%, respectively, which contradicts the former results.29,33

-7Furthermore, Fan et al. showed that the values of the Weibull moduli form a “U” shape in which m is higher in regions I (Vp < 10vol.%) and III (Vp > 55vol.%) than in region II (10vol.% < Vp < 55vol.%), with m values between 4 and 11, based on their extensive experiments on hydroxyapatite, together with data from 16 different research groups and eight different materials, a first-time collection of more than 1500 fractures specimens.29 Nevertheless, Weibull modulus data for highly porous materials is extremely limited in the open literature and there is almost no research systematically investigating the effect of pore size on the Weibull modulus. Consequently, methods capable of revealing the connection between microstructure and Weibull statistics are desirable. To summary, determination of mechanical properties such as fracture strength and Weibull modulus of brittle porous materials (BPMs) is very important for their application in chemical industry and tissue engineering. Nevertheless, experimental researches are short in relating the microstructure with mechanical properties, and they can barely reflect the impact of porosity and pore size separately. Meanwhile, numerical methods based on first order fracture usually focus on the incipient fracture regions, neglecting damage accumulation during load, and hence may produce erroneous results. Lattice spring model (LSM), which is a micromechanical model and algebraically equivalent to a simple finite element method (FEM), is capable of correlating microstructure with mechanical properties of materials.34 It transforms continuum into a network formed by nodes and springs that can be used to simulate the mechanical response of the simplified materials under loading. Local heterogeneity as well as interface discontinuity is represented by diversity in force constants of the springs, meanwhile fractures, microscopic or

-8macroscopic, are imitated by the removal of clusters of springs.35-37 The early lattice model contains only normal springs transmitting central forces and therefore allows a fixed Poisson’s ratio of 0.25 to be obtained.34 By introducing a non-central two-body interaction or a harmonic potential for rotation of bonds from their initial orientation,38 the Poisson’s ratio is rendered variable value with an upper bound of 0.25. However, the rotational invariance is not recovered for large deformation. Zhao et al.39 proposed a rotational invariant 3D distinct LSM (DLSM), which contains multi-body shear type springs in addition to the normal ones, and this model can enlarge the limitation of Poisson’s ratio to 0.3.40 Also, a 4D LSM adopting only central interactions while holding the rotational invariance is proposed recently and is reported to be suitable for solving complex solid mechanics problems.41 Actually, by the introduction of other methods or theories such as Cahn Hilliard method,36 mesoscopic dynamic simulation

42-44

, Zenar model,42

Volume-Compensated Particle Model and Monte Carlo simulation,45,46 LSM is capable of simulating quasi-static/dynamic deformation and fracture of binary polymer blends, block copolymers and polycrystalline materials and even the stress relaxation behavior of nanocomposites. Therefore, the lattice spring model is regarded fairly suitable for exploring mechanical properties of microstructure-sensitive materials. In this work, the effects of porosity as well as pore size on the Weibull modulus of brittle porous material was studied systematically because the two factors greatly affect the microstructure of porous materials and thus their mechanical properties. With this aim, lattice spring model is used to obtain the fracture strength and Weibull modulus by the simulation of the Brazilian test. The first goal was the validation of the model, and this was achieved by

-9comparing the elastic fields with the incipient morphology of porous materials. Initial stiffness was also measured and compared with experimental and analytical results. Moreover, stress-strain curves and Weibull plots for systems with various microstructure were reproduced for the estimation of fracture strength and Weibull modulus. Their dependence on porosity and pore size was studied as well.

2. MODELS AND METHODOLOGY 2.1 Lattice spring model A 3D lattice spring model is used to discretize continuum with infinite degree of freedom into discrete system with finite degree of freedom. The generated system is a network of springs connecting regularly distributed sites or nodes. Each node is linked with 18 other nodes and the initial length as well as alignment of the linking springs varies from each other. Springs connecting the nearest nodes are parallel to the coordinate axes and 1 in length, and those connecting the next-nearest nodes make an angle of 45 degree or 135 degree with axes and are √2 in length. Force imposed by a spring can be calculated through the corresponding matrix and matrices of springs in different alignment can be derived by similarity transformation on 𝑘 𝑀[100] = [0 0

0 𝑐 0

0 0], 𝑐

(1a)

which is the matrix of springs aligned in the positive direction of X-axis. For instance, the matrix corresponding to the [101] direction is obtained by the following equations: 𝜋

𝑅 = 𝑅 (4 ) =

1 −1 0 1 0 ], 0 0 √2

√2 [1 2

(1b)

- 10 𝑘+𝑐 1 𝑀[101] = 𝑅 −1 𝑀[100] 𝑅 = 2 [ 0 𝑘−𝑐

0 2𝑐 0

𝑘−𝑐 0 ]. 𝑘+𝑐

(1c)

All matrices contain two types of force constants, central and noncentral. The central force constant energetically penalizes springs from stretching while the noncentral force constant penalizes the rotation of springs. By comparing the energy stored in a volume element before and after discretization of the continuum, the Young’s modulus, 𝐸, and Poisson’s ratio, 𝜐, are found related to the force constants by the following equations: 𝐸=

5𝑘(2𝑘+3𝑐)

𝜈=

4𝑘+𝑐 𝑐−𝑘 4𝑘+𝑐

,

(2)

,

(3)

where 𝑘 and 𝑐 are the central and noncentral force constants respectively.34 The elastic force due to the change in position of a neighboring node is a linear function of the local displacement. The resultant force acting on the 𝑖th node, 𝐹𝑖 , is given by ∑𝑗 𝑀𝑖𝑗 ⋅ (𝑢𝑖 − 𝑢𝑗 ) = 𝑚𝑖

𝜕2 𝑢𝑖 𝜕𝑡 2

= 𝐹𝑖 ,

(4)

where the vector 𝑢𝑖 is the displacement of the 𝑖th node from its original position, 𝑀𝑖𝑗 is a symmetric 3 × 3 symmetric tensor which describes the interaction between various nodes through central and noncentral force constants, ∑ 𝑗 represent a sum over the 18 nearest- and next-nearest nodes and 𝑚𝑖 is the mass of the 𝑖th node. As all nodes fit the equation, a global stiffness matrix is integrated to describe the relationship between the local displacements of all nodes and all the forces applied to the entire system. The right side of the equations indicates boundary conditions under which the system is loaded; for instance, uniaxial compression. By adjustments to components in the global stiffness matrix as well as the right vector, various boundary conditions is able to be imposed on the studied system to simulate a variety of mechanical-strength measurements.

- 11 When the equations are set up, the conjugate gradient (CG) method is adopted to solve the equations, giving local displacements of all nodes forming the system. The strain fields, 37 stress fields and energy fields can be calculated by the following equations34,38,47 𝜀𝑥,𝑢(𝑖,𝑗,𝑘) =

−𝑢(𝑖+2,𝑗,𝑘) +8𝑢(𝑖+1,𝑗,𝑘) −8𝑢(𝑖−1,𝑗,𝑘) +𝑢(𝑖−2,𝑗,𝑘) 12ℎ

,

1

𝑠 𝜎𝑘𝑙 = 𝐴 × ∑𝑠 𝐹𝑠 ⋅ 𝑝𝑘𝑙 , 𝑇

1

𝐸𝑖 = 2 ∑𝑗(𝑢𝑖 − 𝑢𝑗 ) ⋅ 𝑀𝑖𝑗 ⋅ (𝑢𝑖 − 𝑢𝑗 ),

(5a) (5b) (5c)

where 𝑢(𝑖,𝑗,𝑘) is the displacement field at coordinates 𝑖, 𝑗 and 𝑘, ℎ = 1 is the initial distance between adjacent nodes, 𝜎𝑘𝑙 is the stress acting on the cubic cell, ∑ 𝑠 represents a 𝑠 sum over the cubic surface, 𝐹𝑠 is the force on any surface 𝑠 of the cell, 𝑝𝑘𝑙 is a unit vector

either normal (𝑘𝑙 = 11,22,33) or parallel (𝑘𝑙 = 12,23,13) to the surface 𝑠, 𝐴 is the surface area and 𝐸𝑖 is the elastic energy stored in the 𝑖th cubic lattice. In order to simulate the formation and propagation of fractures, three hypotheses were proposed in this paper: springs were expected to deform in a linear elastic way before breakage in a tensile mode; material degradation could be introduced by reassigning a low elastic modulus to the failed elements rather than moving them from the system48, 49; the total failure of the system resulted from the accumulation of element damages. The rate of failure 𝑝𝑖 (𝑡) is introduced to determine whether spring clusters across the surface are apt to break as a result of tensile stress and the rate of a surface at time t is calculated by the given equation 𝑝𝑖 (𝑡) = [

𝜀𝑖 (𝑡)−𝑊𝑖 𝛽 𝐵𝑖

] ,

(6)

where 𝜀𝑖 (𝑡) is the local strain of a certain surface, 𝑊𝑖 is the minimum value of strain at which element breakage can happen, 𝐵𝑖 is the arbitrary scaling parameter which can be

- 12 chose as 𝐵𝑖 = 𝑊𝑖 and the modulus 𝛽 allows for a nonlinear relationship between failure rate and local strain. The global rate of damage is introduced to estimate if a single element breakage is to happen somewhere in the system and it can be obtained by the summation of rates of failure over all nodes forming the material. If the counted global rate of damage is over 1, a single element breakage is guaranteed to take place inside the material and the location of the impending fracture is determined by both the probability of failure, 𝑃𝑖 (𝑡), and the cumulative probability, 𝑐𝑖 (𝑡), defined as 𝑝 (𝑡)

𝑃𝑖 (𝑡) = ∑𝑁 𝑖 𝑝 (𝑡),

(7a)

𝑐𝑖 (𝑡) = ∑𝑘≤𝑖 𝑃𝑘 (𝑡).

(7b)

𝑖=1 𝑖

Assuming that an element breakage is to occur somewhere in the material, a random number ranging from 0 to 1 will be produced and compared with the cumulative probabilities of all nodes. The surface for which 𝑐𝑖 (𝑡) < 𝑅𝑁𝐷[0,1] < 𝑐𝑖+1 (𝑡) will be chosen to fail and springs across the surface will be reassigned a relatively low elastic modulus by multiply the stiffness of the intact element to a ratio called the residual strength coefficient.48,49 The formation of a local fracture surface is shown in Fig. 2. The elastic energy gradually stored during loading will be transferred into surface energy for newly created cracks or released in the form of a sound emission on account of fracture for real elastic-brittle materials in most cases. However, in this paper, the elastic energy released due to element fracture were mainly dissipated through viscous damping of the neighboring nodes. In order to capture the dynamics of these nodes as well as the entire system, the velocity Verlet algorithm was adopted and it took the positions, velocities, and

- 13 accelerations at time 𝑡 to obtain the same quantities at time 𝑡 + ∆𝑡 in the following way: 𝑟𝑖 (𝑡 + ∆𝑡) = 𝑟𝑖 (𝑡) + 𝑣𝑖 (𝑡)∆𝑡 + ∆𝑡

𝑣𝑖 (𝑡 + 2 ) = 𝑣𝑖 (𝑡) + 𝑎𝑖 (𝑡 + ∆𝑡) = −

𝐹𝑖 (𝑡) 𝑚𝑖

𝐹𝑖 (𝑡) 2𝑚𝑖

𝐹𝑖 (𝑡) 2𝑚𝑖

∆𝑡 2 ,

∆𝑡,

(8b) ∆𝑡

− Θ𝑣𝑖 (𝑡 + 2 ), ∆𝑡

𝑣𝑖 (𝑡 + Δ𝑡) = 𝑣𝑖 (𝑡 + 2 ) +

(8a)

𝐹𝑖 (𝑡+Δ𝑡) 2𝑚𝑖

∆𝑡,

(8c) (8d)

where 𝑣𝑖 , 𝑎𝑖 , and 𝑚𝑖 are the velocity, acceleration and mass of node 𝑖 respectively. The viscous damping constant, namely 𝛩 , was assigned a proper value to result in an underdamped system and the damping would continue throughout the system until all nodes found their balancing positions. The energy consumed through dissipation in the simulation is theoretically equivalent to that in the form of surface energy and acoustic energy. Stress-strain curves were plotted for better understanding of the relationship between local structure and global mechanical properties. The strain of the system was calculated from the displacement of the free boundary and the stress of a certain section was obtained by the summation of local stresses along the section. In the early stage of loading, it is expected that increment in strain is in proportion to that in stress and the slope of the tress-strain curve exhibits the initial stiffness of the entire system. As the material deforms due to the increasing loading, a large number of springs stretch to an extent beyond the fracture threshold and stop functioning, which leads to degradation in the global stiffness and finally the total failure. Previous studies determined the end of loading when a crack propagated the length of the material, rendering it two separate pieces, and it is reasonable as well as practical for homogeneous systems or those with simple structure. However, it does not suit materials with complicated structure especially for porous materials originally containing percolating

- 14 cavities. Therefore, in this paper, the end of loading was determined by stress and the system was considered failed when the section stress for timestep 𝑡 was lower than 85% of the peak stress37.

2.2 Weibull distribution In practice, the Weibull modulus of brittle materials such as catalyst supports is recognized only after their fracture strength and the mechanical strength distributions are usually estimated by means of Weibull statistics derived from the “weakest-link hypothesis”, which states that mechanical failure is recognized to be microstructure-sensitive, and the weakest spot under stress determines the strength of a brittle material. Moreover, the Weibull distribution is usually considered the best choice for brittle materials because it is bounded (the lowest possible fracture strength is zero), it can provide reasonably accurate failure forecasts with small numbers of test specimens and it provides a simple and useful graphic plot. In the specific case of Weibull fracture strength analysis, the cumulative probability function is written such that the probability of failure, 𝑃𝑓 , increases with the fracture stress variable, 𝜎,: 𝜎−𝜎𝑢 𝑚

𝑃𝑓 = 1 − exp [− (

𝜎𝜃

) ].

(9)

This is known as the Weibull three parameter strength distribution, which is directly derived from weakest link theory. The threshold stress parameter, 𝜎𝑢 , represents a minimum stress below which a test specimen will not break. 𝜎𝜃 is called the scale parameter or characteristic strength, which is dependent on the stress configuration and test specimen size,

- 15 and 𝑚 the Weibull modulus, or the distribution shape parameter. However, it is very risky to assume a finite threshold strength exists without careful screening or nondestructive evaluation. Hence, the two-parameter form is most commonly used for simplicity. Setting 𝜎𝑢 to zero in Eq. (2) and taking the double logarithm of the resulting two-parameter Weibull distribution yields: 𝜎

𝑚

𝑃𝑓 = 1 − exp [− (𝜎 ) ], 𝜃

ln [ln (

1 1−𝑃𝑓

)] = 𝑚(ln𝜎 − ln𝜎𝜃 ).

(10a) (10b)

The reason the double logarithm of the Weibull equation is used in strength analysis is the ease of accessing information. In this regard, the Weibull modulus determines the slope of the distribution function and characterize the scattering extent of failure strength. One can obtain the Weibull modulus by calculating the slope of the Weibull plots. A high Weibull modulus (steep slope, less spread) is desirable for all materials since it indicates an increased homogeneity in the flaw population and a more predictable failure. Conversely, a low Weibull modulus is indicative of a large spread within the group and a less predictable failure behavior. Besides, The characteristic strength, 𝜎𝜃 , also is a location parameter that indicates distribution location along the abscissa X axis: a large 𝜎𝜃 shifts the data to the right, while a small 𝜎𝜃 shifts the data to the left. Furthermore, the characteristic strength is the strength value, 𝜎, at 𝑃𝑓 = 63.2%, when the left side of Eq. (10a) = 0, thus reported Weibull characteristic strength values (𝑃𝑓 = 63.2%) are slightly greater than the mean strength values (𝑃𝑓 = 50%). In order to obtain adequate Weibull strength distribution parameters, it is common to require no fewer than ten test specimens and preferably 30 to obtain good estimates of the

- 16 Weibull modulus, with more test specimens contributing little towards better uncertainty estimates.50 After simulation, the collected stress values in terms of the same test configuration were ranked in ascending order, 𝑖 = 1, 2, 3, … , 𝑁, where 𝑁 is the total number of test specimens and 𝑖 is the 𝑖th datum. Thus, the lowest stress represents the first value (𝑖 = 1), the next lowest stress value is the second datum (𝑖 = 2), etc., and the highest stress is represented by the Nth datum. As linear regression (LR) analysis was used in the fitting process with the true value of 𝑃𝑓 (𝜎𝑖 ) for each 𝜎𝑖 unknown, a probability estimator was necessary to assign a ranked probability of failure to each datum. A commonly used estimator that has low bias when with linear regression analyses is50 𝑃𝑓 (𝜎𝑖 ) =

𝑖−0.5 𝑁

,

(11)

where 𝑖 is the 𝑖th datum and 𝑁 is the total number of data points. This is also the estimator used in the paper. Many studies have shown that for 𝑁 > 20, this estimator produces the least biased estimates of the Weibull parameters.

3. RESULTS 3.1 Morphology The Brazilian test was adopted in this work to obtain the mechanical property of brittle porous materials instead of the uniaxial compressive strength (UCS) test. This choice was based on the conclusion drawn by Hallbauer et al.51 and Tang52 that the formation of macroscopic fracture plane is more likely the result of a tensile process rather than a shear one, where the shear stress is found small during uniaxial compression53, and both splitting

- 17 and shear failures usually observed in the axial compressive experiments of brittle materials may be the final manifestation of earlier tensile crack growth. Therefore, the Brazilian disc test was simulated to focus only on the fracture process promoted by tensile stress. Fig. 3 shows a system composed of two linear elastic bodies i.e. the cylinder-shaped specimen (radius 𝑟 = 32, height 𝑤 = 64) and plane jaws (or curved ones of predefined curvature) proposed for the Brazilian test. The contact angle, 𝜔0 , was set to be 8.8° and as a result the contact length between the sample and the metallic jaw, 𝑙 = 2𝜔0 𝑟, was close to 10. To simplify the problem, the central part of the specimen, with its size being 10×64×64, was isolated as an object for investigation (shown as the inset in the upper right corner of Fig.3). This simplification can be justified by the fact that the displacement field of the central part is insensitive to the exact distribution of the external load along the contact arc54 and pellets as well as discs breaking at small loads typically split off along the longitudinal section during the Brazilian test.23 To investigate the influence of porosity on the mechanical response of brittle porous materials, the plates were embedded with pores, which were uniform in size as 2, and the porosity could be adjusted by the number of pores introduced. Besides, the impact of pore size was studied by embedding pores of different sizes, namely 1, 2, and 3, while making the porosity fixed at 20vol.%. It should also be noted that in addition to pore size, grain size also has eminent impact on elastic properties such as Young’s modulus. As reported by Cleveland and Bradt,55 a sharp decrease in stiffness is detected when the average grain size surpass a critical value. This value is about 1𝜇𝑚 for Al2TiO5, 3𝜇𝑚 for Fe2TiO5 and 5𝜇𝑚 for MgTi2O5. Moreover, Lu et al.56 confirmed by a 2-D FEM analysis that the stress intensity factor shows

- 18 obvious increase when the ratio of grain size to pore size is under 0.4. Therefore, in this work the grain size is comparable to the pore size and is set to be about 15𝜇𝑚, assumed that the diameter of the specimen is about 1𝑐𝑚. Systems considered in this paper are listed in Table. I.

Table 1. Materials Evaluated and Their Specifications material porosity pore size (vol.%) (a.u.) P0-S0 0 0 P10-S1 10 1 P10-S2 10 2 P10-S3 10 3 P20-S1 20 1 P20-S2 20 2 P20-S3 20 3 P30-S1 30 1 P30-S2 30 2 P30-S3 30 3 P40-S1 40 1 P40-S2 40 2 P40-S3 40 3 P50-S1 50 1 P50-S2 50 2 P50-S3 50 3

Fig. 4 shows the morphology and structure along the longitudinal section of plates of different porosity. The solid phase is depicted in white and pores are depicted in blue. The average distances between two pores decreases as porosity increases and a transition from dispersed to bicontinuous structure is observable. The produced morphologies can then be adopted as inputs for simulating the mechanical response of porous materials through LSM. The global stiffness, local strain/stress field and even the fracture and degradation process of complex brittle porous materials can be examined by assigning different force constants to

- 19 the two phases.

3.2 Normal Strain/Stress Field At the beginning of loading, each sample was deformed by the application of a global strain along the Z-axis at the boundary to examine the elastic fields throughout the systems. Within the produced lattice model, there are three kinds of springs depending on the two grids they link. The linked pair of grids can be identical as solids (A-A) or pores (B-B), or different with a solid element and a void on each side (A-B). All the springs obey linear elastic relation during the mechanic simulation regardless of the grids they link. To obtain the properties of the entire system rationally, proper parameter assignment is required before simulation. In the work of Tang et al., the residual strength coefficient was set to be 0.1.48 Analogously, the Young’s modulus of the solid phase, 𝐸A , was set to be 12.5, which was 5 times that of the void phase, 𝐸B , and the stiffness of the interfacial springs, 𝐸I , was half the value of the geometric means of 𝐸A and 𝐸B . Although the Young’s modulus of voids should be zero theoretically, application of the approximation in stiffness is sufficient to reflect disparity of the two phases while avoiding singularity during solution and improving convergence as well as computing speed. As a result, difference in elastic moduli undoubtedly causes complexity in elastic fields within the system. Fig. 5 depicts the strain fields of the longitudinal section for plates with various porosity during initial loading. The blue regions and the green regions corresponded to the stiff solid phase of low strain and the much-more-ductile macro-pore phase of relatively large deformation, respectively. Indeed, the pattern of the strain field reflects the morphology of the

- 20 section investigated and is in accord with the structure information shown in Fig. 3. Besides, a reduction in tensile strain was observable for solid elements adjacent to pores compared with those within a solid bulk and this relaxation mostly occurred along the direction parallel to the Z-axis, which is considered as a result of the blocking effect by pores. However, it should also be noted that for real materials, strain at the boundary between pores and solid must be zero rather than just reduce. This inconformity lies in the adoption of residual strength and the difference can be eliminated by assigning an infinitesimal value to the residual strength coefficient. The probability distribution function of normal strain field is shown in Fig. 6. All functions were plotted from data averaged over 5 independent runs. As for the P10-S2 group, a large and narrow peak is observed at a relatively low value of strain, corresponding to the deformation of the major stiff solid phase, and the tailing is attributed to that of the ductile minor pore phase. As porosity increases, the peak shifts toward low strain value, dwindling in height and growing in peak width. This conversion can be explained by the space distribution of pores. For a plate with relatively low porosity, pores are usually separated from each other, as is the case shown in Fig.4a, and their deformations are readily inhibited by the ambient solid phase, thus resulting in a relatively low strain value compared with that of a homogeneous system built up merely by the pore elements, which is also referred to in the work of Buxton and Balazs.57,58 However, as porosity grows, the distribution of the pores becomes much more complicated. Pores can stay isolated, attach each other to form local clusters of closely spaced pores, or even penetrate the whole system when the porosity is high enough. In fact, as Caiulo and Kachanov report, the maximal value of stress intensity factors

- 21 (SIFs) was noticeably elevated by clustering and this elevation had substantial fluctuation from one crack array to another, in the range of 10%–30%.59 It may also be the case for clusters of pores and therefore lead to the disparity in the local deformations and stress of the elements and cause the conversion of the peaks mentioned above. Fig .7 depicts the normal stress fields of the longitudinal section for plates with different porosity. Similarly, the pattern of the stress field matches the origin structure of the corresponding system. In contrast, the blue regions are attributed to pores of low stress, which should be zero for an ideal cavity, though, and the green regions correspond to the highly stressed solid phase. Cannillo et al.60,61 carried out a simulation of uniaxial tensile test on model porous glasses with closed, isolated circular pores by utilizing a 2-D simulation code called OOF as well as the finite element method. The effect of pore size, pore clustering, volume fraction and distance between two pores on the stress field were investigated thoroughly. It was confirmed that if two pores are close enough, the two stress field interact and stress concentration can be detected at pore tips aligned perpendicular to the direction of loading applied. Moreover, the stress concentration is intensified as the distance between pores decreases. This phenomenon can also be detected in the region marked by a red circle in Fig.7. Specifically, the solid elements situated in the narrow gaps surrounded by pores usually possess stress much higher than the other solid elements, indicating the presence of stress concentration. Fig. 8 shows the probability and cumulative distribution functions of local stress field. Regarding the P10-S2 group, there was a large peak at a relatively high value of stress

- 22 attributed to the stress inside the major solid phase, and a shoulder peak at a low value of stress corresponding to that inside the cavity. A scattering was detected at the low stress region due to the difference in the environment where pores were situated. As more pores were introduced into the system, stress transfer of the solid phase was blocked so that the stress as well as the strain inside the majority of the solid decreased. This is indicated by the shift in the leading peak to a lower stress value. Besides, the shape of the curve changed into a trapezium for the P20-S2 and P30-S2 sample and finally into a tailing peak for the P50-S2 sample, indicating an essential change in stress distribution within the range investigated. It has to be pointed out that despite the decline in the section stress, the disparity between the maximum and the most probable value in local stress is enlarged along with increased porosity, leading to a transition from instantaneous fracture mode to progressive fracture mode during late stage loading.

3.3 Constitutive relation To study the fracture and degradation process in the ensuing simulation, the porous plates were further loaded at a strain rate of 6.24×10-4. A strain-based fracture criterion was utilized, and the lower bound fracture threshold was 0.03, 100 and 0.04 for the A-A, B-B and A-B spring, respectively. This guarantees that local fracture predominately take place in the solid phase and never occur in the pores. Each spring fractures just for once at most, and the broken bond possesses parameters the same as those of the springs linking two voids. At each timestep, the entire strain field was checked and springs possessing strain exceeding the lower bound threshold were screened and selected to break according to Eqs. (6), (7a) and

- 23 7(b). Section stress was measured at each step, and the system was considered failed once the section stress was less than 85% the value of peak stress.

Table 2. Mean Young’s Modulus, Fracture Strength and Weibull Parameters of the Porous Materials Determined by the Brazilian Test material P10-S1a P10-S2a P10-S3a P20-S1a P20-S2a P20-S3a P30-S1a P30-S2a P30-S3a P40-S1a P40-S2a P40-S3a P50-S1a P50-S2a P50-S3a P10-S1b P10-S2b P10-S3b P20-S1b P20-S2b P20-S3b P30-S1b P30-S2b P40-S1b P40-S2b P50-S1b P50-S2b

sample capacity 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 20 20 20 20 20 20 20 20 20 20 20 20

Young’s modulusa 10.10(0.05) 9.9(0.1) 9.8(0.1) 8.93(0.04) 8.57(0.08) 8.5(0.1) 7.84(0.02) 7.37(0.09) 7.3(0.1) 6.8(0.04) 6.35(0.06) 6.3(0.1) 5.88(0.05) 5.39(0.07) 5.3(0.1) – – – – – – – – – – – –

fracture strengtha – – – – – – – – – – – – – – – 0.173(0.002) 0.159(0.002) 0.143(0.003) 0.143(0.002) 0.132(0.003) 0.130(0.004) 0.113(0.001) 0.113(0.002) 0.096(0.001) 0.098(0.002) 0.082(0.001) 0.086(0.002)

a

Standard deviation from the mean in parentheses.

b

Standard error of the mean in parentheses.

Weibull modulusb – – – – – – – – – – – – – – – 128(6) 82(5) 62(3) 96(6) 65(3) 36(2) 103(7) 64(5) 106(6) 71(5) 118(7) 64(4)

characteristic strength – – – – – – – – – – – – – – – 0.173 0.160 0.145 0.143 0.133 0.133 0.114 0.114 0.096 0.099 0.082 0.087

- 24 -

Fig. 9 indicates the stress-strain curves of plates with various pore distribution. The slope of the linear part is a measurement of the incipient stiffness of the entire porous system, and the specific values are given in Table II. As is obvious in Fig. 10, the global Young’s modulus varies almost linearly with growing in the volume fraction of more deformable phase, which is consistent with previous findings observed by the utilization of LSM.37,57,58 However, it should be noted that in terms of porous ceramics or other BPMs the porosity dependence of Young’s modulus seems to follow a power law type rather than a linear one. Actually, the empirical equation first coined by Gibson and Ashby62, 𝐸𝑝 = 𝐸𝑑 ∗ (1 − 𝑝)𝑐 , is used to describe the relationship between the Young’s modulus 𝐸𝑝 of a porous body and its dense material properties 𝐸𝑑 .53, 63 In this relationship the pore shape factor 𝑐 is adopted in addition to the porosity, which is a scalar parameter, for an adequate description of the influence of the pore structure on elastic properties. According to Bruno et al., the pore shape factor of systems containing spherical pores (open or closed) is about 2. 53, 63 The value of Young’s modulus in Fig. 10a is generally larger than those predicted by the former empirical equation and it may result from the relatively large residual strength coefficient. Therefore, Young’s moduli obtained under the condition that 𝐸𝐴 = 10𝐸𝐵 were fitted by the power law type relationship and the result is depicted in Fig.10b. The good consistency certifies the capability of LSM to producing reliable prediction of BPMs. Contrary to the porosity dependency, pore size doesn’t show great impact on the Young’s modulus. The stiffness does reduce when the pore size grows from 1 to 2, and the discrepancy enlarges as porosity increase. However further growth in pore size doesn’t make

- 25 recognizable difference. As the plate was further squeezed along the Z-direction, the local strain (𝑢𝑥𝑥 , where 𝑢𝑖𝑗 is the strain tensor) inside elements reached the threshold bond value and incipient bond breakage occurred. The larger the porosity, the earlier the initial fractures took place. Fig.11 and Fig .12 show the number of broken bonds as well as its cumulative distribution at each timestep in the form of strain. Fig.11a depicts the fracture process of a pure solid (P0-S0) system and Fig.11b represents that of a P20-S2 sample. Stress drop is observable in both plots as a result of the brittle fracture nature of elements. When a cluster of springs across the fracture surface was judged to fail, the previous stored elastic energy would be released, accompanying drastic degradation in the spring stiffness. Besides, grids originally linked to each other would move to new positions to maintain force balance, which ultimately led to the stress drop in section stress. Therefore, it is reasonable to anticipate that the stress drop will be in proportional to the number of broken bonds across the section and this assumption is proved to be rational by the plot in Fig.11b. Especially, for an ideal homogeneous elastic-brittle material, the stress drop ensuing incipient fracture is much more pronounced compared with that of porous systems. Moreover, the fact that the cumulative distribution of the homogeneous system is a step function (see Fig.12a) indicates the occurrence of instantaneous fatal failure, which is in accord with experimental findings.49 In the recent work of Pandey et al., which investigated the uniaxial tensile response of porous ceramic materials with porosity over 50vol.% i.e. cordierite, it was found that these materials displayed non-linear tensile behavior such that the secant modulus decreased with increasing strain and the degree of nonlinearity is related to the initial microcrack density and

- 26 evolution of damage in the material.64 By the usage of simple linear elastic elements, LSM can reproduce linear behavior of homogeneous materials as well as non-linear behavior of composites or porous system, which is a significant advantage for simulating the entire fracture process. For all systems considered in this study, after initial fracture occurred, the section stress continued growing with increase in load within a certain range. Meanwhile, the defects inside the system kept accumulating, resulting in a deviation from linear behavior, and the deviation is most evident around the peak of the curve (see Fig. 9), which accords with the experimental results in the work of Pandey et al.64 Actually, the method of obtaining non-linear behavior by means of linear elements is not a theoretical assumption only,65,66 it is experimental proved by the usage of Electron Backscattered Diffraction as well. 67 As Bruno et al. report, the lattice deformation of aluminum titanate (AT) is approximately linear elastic during uniaxial compression tests, compared with the non-linear macroscopic behavior under high stresses.67 The value of fracture stress is impacted by porosity together with void size according to Fig. 13 and Table II. It is clear that for a range of operation conditions, increment in porosity leads to weakening in the peak stress and thus, the fracture stress. Besides, the obtained data are compared to the empirical equation, 𝜎𝑝 = 𝜎𝑑 ∗ (1 − 𝑝)𝑐 , where the pore shape factor 𝑐 equals to 2.53,

63

Deviation from the curve is detected from Fig. 13b, which means that

porosity dependence of fracture stress differs from that of stiffness and the former is stronger. This phenomenon was also found by Shyam et al.68 In fact, the quantitative correlation between the porosity dependence of elastic properties and fracture properties is not clear, and the underlying reason for the absence of correlation lies in the entirely different sensitivity to

- 27 local features of the crack field geometry.59, 63 In terms of samples with identical porosity, systems with large-sized pores embedded exhibit degraded strength, while those containing small-sized pores is more sensitive to changes in pore size. This trend is qualitatively in accord with the findings of Liu et al.18 In their work porous hydroxyapatite (HAp) ceramics were prepared utilizing pore-forming agent of different particle size (corresponding to different macropore sizes). The as-prepared samples, whose porosity ranged from 33vol.% to 78vol.%, were then investigated for their compressive strength and it was found that the compressive strength decreases linearly with increase macropore size. On the other hand, fracture strength stays almost the same when the pore size increase from 2 to 3, which is similar to the trend regarding the stiffness. After the peak, section stress decreases regardless of further increment in load and the relation between strain and stress maintains non-linear. In addition to the pre-peak behavior of porous materials, the post-peak behavior is also influenced by porosity, as can be seen in Fig. 9. In terms of the systems whose porosity is less than 30vol.%, slight increase in load leads to dramatic decline in section stress, and defects can propagate and penetrate the entire system to cause fatal failure in a short period of time, which is characteristic of the brittle fracture mode. However, highly porous systems are different. Smolin et al. carried out a 3D computer simulation of mechanical behavior of a brittle porous material under uniaxial by the use of movable cellular automata (MCA) method in their work. They report that porous 3D samples can demonstrate quasi-viscous regime of fracture and such regime of fracture is determined entirely by the pore space structure and observed for 3D samples having porosity greater than 40%.69 Similarly, in this simulation

- 28 section stress is not sensitive to load during the post-peak region and it remain almost the same despite an obvious increment in tensile strain along the X-direction. It means that ceramics with large porosity stochastically distributed in space can miss the stage of rapid propagation of the main crack and a sharp drop in the strength, which is typical for brittle materials.69 This phenomenon derives from both the introduction of voids and the adoption of elements with residual strength. Introduced pores entangle the elastic fields, causing more element breakages throughout the whole system while diminishing the average breakage number of each section at each timstep, which inhibits the formation of instant fatal failure along a certain section as is the case in homogeneous system.

3.4 Weibull statistics for fracture data In this section, seven groups of systems with various pore distribution were investigated, and for each group, 20 similar samples were tested to obtain the Weibull parameters. In previous works of simulations, the analysis is limited to the first failure on basis of the weakest link assumption.15,70 However, in this paper, statistics of high order failures can be analyzed by the use of LSM, which allows damage accumulation and localized compliance reduction during loading.70 Fig. 14 depicts the Weibull probability plots for each of the 7 groups of specimens. The Young’s modulus of the solid phase is set to be 12.5, and stresses for solids of different stiffness can be determined by scaling with linear elasticity. All the plots fabricated in this study are relatively linear due to the fact that the synthetic microstructure possesses only one defect population. In Fig. 14a, the characteristic strength, together with average strength,

- 29 decreases along with increase in porosity, which is represented by the shift to lower value along the lateral axis. In Fig. 14b, a slight degradation in characteristic strength is observed when the pore size grows from 1 to 2, whereas no difference is detected when the pore size further grow to 3.These findings are in consistence with those revealed in Fig.13. In addition to the characteristic strength, the Weibull modulus, which is the slope of the Weibull plot, is also an important measurement of the reliability of a group of identical specimens. Fig. 15 shows the Weibull modulus as a function of porosity and pore size. As the pore size is fixed at 2, the Weibull modulus declines from 82 (P10-S2) to 65 (P20-S2), and remains almost unchanged as porosity increases. Similar trend can be detected for Weibull moduli of systems whose pore size is 1. Although the value of the Weibull moduli calculated in this study is much larger, the extreme difference is quite close to that measured by experiments.24,29 Additionally, the tendency to degrade in Weibull modulus also consists with the 𝑚~𝑃 behavior in Region I and II of the “U” shaped plot revealed by Fan et al.29 On the other hand, the Weibull modulus decreases from 96 (P20-S1b) to 36 (P20-S3b), and the disparity (about 59) in Weibull modulus caused by changes in pore size is twice more than that caused by variation in porosity. By examining the Weibull modulus, it helps to distinguish samples with higher uniformity from those possessing similar characteristic strength. For instance, the characteristic strengths are approximately 0.133 for both the P20-S2b and P20-S3b group. However, the Weibull modulus differs a lot from each other, which means the range and dispersion of the data for the P20-S2b group is much narrower, and hence the better mechanical performance.

- 30 -

4. DISCUSSION Our findings suggest that LSM is capable of reproducing the damage accumulation of brittle porous materials and predicting their fracture strength as well as Weibull modulus. Porosity shows significant impact on fracture strength while pore size dominates the Weibull modulus. We show that lattice spring model is an efficient method for interrelating the microstructure, the mechanical characteristics of the components and the macroscopic behavior of heterogeneous material. It therefore enables us to establish how choices made in the microstructure to meet the demand of materials functioning in various operating conditions. For example, solid catalysts with a high Weibull modulus is beneficial to their applications in chemical industry,22 and the optimization can be attained by usage of small-sized pore-forming agent; whereas fracture strength is the main obstacle for commercial implementation of inorganic membrane reactor utilized in chemical synthesis and separation,9 and reducing the porosity will be effective in improving the strength. However, it should also be noted that reduction in porosity will reduce the catalytic or the filtration activity. Therefore, investigation regarding the dependence of fracture stress on other factors, such as pore shape,53 is of great significance. In previous reported works, the Weibull modulus ranges from 4 to 25,16,24,29,71 which is different from that obtained in this study (over 60 for most of the samples). This phenomenon is caused by two reasons. Firstly, the Weibull modulus calculated in this paper is obtained from statistics of high order failures. The fracture stress, as well as the peak stress, is the result of a sequence of local fractures rather than exclusively determined by the first failure,

- 31 which is predominantly affected by the specific morphology and thus possesses relatively large numerical dispersion. Secondly, the system simulated is under quasi-static conditions and in equilibrium during the total process regardless of the occurrence of local fracture, while unstable crack propagation often occurs ensuing the formation of micro-crack as is usually the case in experimental tests, which also leads to relatively large scattering in the fracture stress measured compared to that calculated in this study. Samples are assumed to fail in the form of longitudinal splitting in this paper. However, this assumption may not be applicable for highly porous materials. Therefore, attention should be paid when using LSM to predict the mechanical properties of materials whose porosity is over 50%. Further research should include simulation of other tests suitable for high porosity materials (e. g. biaxial flexural testing). Besides, thermal-mechanical coupling LSM is expected to provide further insight into the operational life span of brittle porous materials functioning under elevated temperature or thermal gradient, as is the case for catalyst supports and membrane reactors.9

5. SUMMARY AND CONCLUSIONS In this study, a micromechanical model, namely LSM, is adopted to simulate the Brazilian test of brittle porous materials and investigate the impact of porosity and pore size on the global properties such as stiffness, fracture strength and Weibull modulus. The system, which is taken from the central part of the cylinder employed in the Brazilian test, is described by a three dimensional network comprised of grids and springs, and the force constants of the spring are adjustable to reflect local heterogeneity. The main findings are

- 32 summarized as follows: 1) Both fields, where stiffer solid phase is found to possess relatively small strain and large stress, accord with the morphology of the corresponding system. 2) The increment in porosity increases the complexity of elastic fields, expanding the range of the local strain value, and even changes the shape of the probability distribution function of local normal strains/stresses. 3) LSM is able to reproduce the stress drop, which is characteristic of brittle materials and caused by local fracture, during the loading process. In addition, it is also proved to be capable of modeling the non-linear behavior resulting from the accumulation of local damage through the use of simple linear elements. 4) Increment in either porosity or pore size can lead to degradation in global stiffness and fracture stress, and changes in porosity have a more pronounced impact. 5) Weibull modulus declines a lot when the porosity is as low as 10vol.%, and remains almost unchanged for the 20~50vol.% range. Besides, the Weibull modulus is much more sensitive to increment in pore size for systems of relatively low porosity, which is quite different from the facts when it comes to the characteristic strength. 6) Porosity have pronounced effect on the average fracture strength while pore size dominates the Weibull modulus of the system. Therefore, it is suggested that porosity and pore size are investigated separately according to the specific operating environment when designing new porous materials or composites.

ACKNOWLEDGEMENTS

- 33 This work was supported by the National Natural Science Foundation of China (No. 21476071) and Shanghai Leading Academic Discipline Project (B502).

REFERENCES (1)

Rámila, A.; Padilla, S.; Muñoz, B.; Vallet-Regí, M. A New Hydroxyapatite/Glass Biphasic Material: In Vitro Bioactivity. Chem. Mater. 2002, 14, 2439-2443.

(2)

Li, X.; Wang,X.; Chen, H.; Jiang, P.; Dong, X.; Shi, J. Hierarchically Porous Bioactive Glass Scaffolds Synthesized With a PUF and P123 Cotemplated Approach. Chem. Mater. 2007, 19, 4322-4326.

(3)

Studart, A. R.; Gonzenbach, U. T.; Tervoort, E.; Gauckler, L. J. Processing Routes to Macroporous Ceramics: A Review. J. Am. Ceram. Soc. 2006, 89, 1771-1789.

(4)

Armor, J. N. Applications of Catalytic Inorganic Membrane Reactors to Refinery Products. J. Membr. Sci. 1998, 147, 217-233.

(5)

Lo, Y. W.; Wei, W. C. J.; Hsueh, C. H. Low Thermal Conductivity of Porous Al 2O3 Foams for SOFC Insulation. Mater. Chem. Phys. 2011, 129, 326-330.

(6)

Lacroix, J.; Jallot, E.; Nedelec, J. M.; Lao, J. Influence of Glass Scaffolds Macroporosity on the Bioactive Process. J. Phys. Chem. B 2013, 117, 510-517.

- 34 -

(7)

Chung, J. J.; Li, S.; Stevens, M. M.; Georgiou, T. K.; Jones, J. R. Tailoring Mechanical Properties of Sol–Gel Hybrids for Bone Regeneration through Polymer Structure. Chem. Mater. 2016, 28, 6127-6135.

(8)

Wang, C.; Shen, H.; Tian, Y.; Xie, Y.; Li, A.; Ji, L,; Niu, Z.; Wu, D.; Qiu, D. Bioactive Nanoparticle–Gelatin Composite Scaffold with Mechanical Performance Comparable to Cancellous Bones. ACS Appl. Mater. Inter. 2014, 6, 13061-13068.

(9)

Li, H.; Song, J.; Tan, X.; J, Y.; Liu, S. Preparation of Spiral Porous Stainless Steel Hollow Fiber Membranes by a Modified Phase Inversion–Sintering Technique. J. Membr. Sci. 2015, 489, 292-298.

(10) Negahi Shirazi, A.; Fathi, A.; Suarez, F. G.; Wang, Y.; Maitz, P. K.; Dehghani, F. A Novel Strategy for Softening Gelatin–Bioactive-Glass Hybrids. ACS Appl. Mater. Inter. 2016, 8, 1676-1686.

(11) Wu, D.; Song, L.; Zhang, B.; Li, Y. Effect of the Mechanical Failure of Catalyst Pellets On the Pressure Drop of a Reactor. Chem. Eng. Sci. 2003, 58, 3995-4004.

(12) Ryshkewitch, E. Compression Strength of Porous Sintered Alumina And Zirconia. J. Am. Ceram. Soc. 1953, 36, 65-68.

- 35 -

(13) Fan, X.; Case, E. D.; Ren, F.; Shu, Y.; Baumann, M. J. Part II: Fracture Strength And Elastic Modulus as a Function of Porosity for Hydroxyapatite and Other Brittle Materials. J. Mech. Behav. Biomed. 2012, 8, 99-110.

(14) Le Huec, J. C.; Schaeverbeke, T.; Clement, D.; Faber, J.; Le Rebeller, A. Influence of Porosity on the Mechanical Resistance of Hydroxyapatite Ceramics Under Compressive Stress. Biomaterials 1995, 16, 113-118.

(15) Keleş, Ö.; García, R. E.; Bowman, K. J. Stochastic Failure of Isotropic, Brittle Materials with Uniform Porosity. Acta Mater. 2013, 61, 2853-2862.

(16) Keleş, Ö.; García, R. E.; Bowman, K. J. Deviations from Weibull Statistics in Brittle Porous Materials. Acta Mater. 2013, 61, 7207-7215.

(17) Rice, R. W. Comparison of Stress Concentration Versus Minimum Solid Area Based Mechanical Property-Porosity Relations. J. Mater. Sci. 1993, 28, 2187-2190.

(18) Liu, D. Influence of Porosity and Pore Size on the Compressive Strength of Porous Hydroxyapatite Ceramic. Ceram. Int. 1997, 23, 135-139.

(19) Li, J.; Lin, H.; Li, J. Factors that Influence the Flexural Strength of Sic-Based Porous

- 36 Ceramics Used for Hot Gas Filter Support. J. Eur. Ceram. Soc. 2011, 31, 825-831.

(20) Rice, R. W. Comparison of Physical Property-Porosity Behavior with Minimum Solid Area Models. J. Mater. Sci. 1996, 31, 1509-1528.

(21) Zhang, J.; Malzbender, J. Mechanical Characterization of Micro- and Nano-Porous Alumina. Ceram. Int. 2015, 41, 10725-10729.

(22) Wu, D.; Zhou, J.; Li, Y. Distribution of the Mechanical Strength of Solid Catalysts. Chem. Eng. Res. Des. 2006, 84, 1152-1157.

(23) Wu, D.; Zhou, J.; Li, Y. Mechanical Strength of Solid Catalysts: Recent Developments and Future Prospects. AIChE J. 2007, 53, 2618-2629.

(24) Fan, X.; Case, E. D.; Gheorghita, I.; Baumann, M. J. Weibull Modulus and Fracture Strength of Highly Porous Hydroxyapatite. J. Mech. Behav. Biomed. 2013, 20, 283-295.

(25) Fu, Q.; Saiz, E.; Rahaman, M. N.; Tomsia, A. P. Toward Strong and Tough Glass and Ceramic Scaffolds for Bone Repair. Adv. Funct. Mater. 2013, 23, 5461-5476.

(26) Pissenberger, A.; Gritzne, G. Preparation and Properties of Niobia-and Tantala-Doped Orthorhombic Zirconial. J. Mater. Sci. Lett. 1995, 14, 1580-1582.

- 37 -

(27) Kishimoto, A.; Koumoto, K.; Yanagida, H.; Nameki, M. Microstructure Dependence of Mechanical and Dielectric Strengths - I. Porosity. Eng. Fract. Mech. 1991, 40, 927-930.

(28) Nanjangud, S. C.; Brezny, R.; Green, D.J. Strength and Young's Modulus Behavior of a Partially Sintered Porous Alumina. J. Am. Ceram. Soc. 1995, 78, 266-268.

(29) Fan, X.; Case, E. D.; Ren, F.; Shu, Y.; Baumann, M. J. Part I: Porosity Dependence of the Weibull Modulus for Hydroxyapatite and Other Brittle Materials. J. Mech. Behav. Biomed. 2012, 8, 21-36.

(30) Ruys, A. J.; Wei, M.; Sorrell, C. C.; Dickson, M. R.; Brandwood, A.; Milthorpe, B. K. Sintering Effects on the Strength of Hydroxyapatite. Biomaterials 1995, 16, 409-415.

(31) Radovic, M.; Lara-Curzio, E. Mechanical Properties of Tape Cast Nickel-Based Anode Materials for Solid Oxide Fuel Cells Before and After Reduction in Hydrogen. Acta Mater. 2004, 52, 5747-5756.

(32) Atkinson, A.; Bastid, P.; Liu, Q. Mechanical Properties of Magnesia–Spinel Composites. J. Am. Ceram. Soc. 2007, 90, 2489-2496.

(33) Frandsen, H. L.; Ramos, T.; Faes, A.; Pihlatie, M.; Brodersen, K. Optimization of the

- 38 Strength of SOFC Anode Supports. J. Eur. Ceram. Soc. 2012, 32, 1041-1052.

(34) Buxton, G. A.; Care, C. M.; Cleaver, D. J. A Lattice Spring Model of Heterogeneous Materials with Plasticity. Modell. Simul. Mater. Sci. Eng. 2001, 9, 485-497.

(35) Buxton, G. A.; Balazs, A.C. Modeling the Dynamic Fracture of Polymer Blends Processed under Shear. Phys. Rev. B 2004, 69, 054101.

(36) Buxton, G. A.; Balazs, A.C. Micromechanical Simulation of the Deformation and Fracture of Polymer Blends. Macromolecules 2004, 38, 488-500.

(37) Deng, S.; Zhao, X.; Huang, Y.; Han, X.; Liu, H.; Hu, Y. Deformation and Fracture of Polystyrene/Polypropylene Blends: a Simulation Study. Polymer 2011, 52, 5681-5694.

(38) Hassold, G. N.; Srolovitz, D. J. Brittle Fracture in Materials with Random Defects. Phys. Rev. B 1989, 39, 9273-9281.

(39) Zhao. G.; Fang, J.; Zhao. J. A 3D Distinct Lattice Spring Model for Elasticity and Dynamic Failure. Int. J. Numer. Anal. Meth. Geomech. 2011, 35, 859-885.

(40) Zhao, G.; Khalili, N.; Fang, J.; Zhao, J. A Coupled Distinct Lattice Spring Model for Rock Failure under Dynamic Loads. Comput. Geotech. 2012, 42, 1-20.

(41) Zhao, G. Developing a Four-Dimensional Lattice Spring Model for Mechanical

- 39 Responses of Solids. Comput. Methods Appl. Mech. Engrg. 2017, 315, 881-895.

(42) Deng, S.; Huang, Y.; Xu, S.; Lin, S.; Liu, H.; Hu, Y. Mechanical Properties of High-Performance Elastomeric Nanocomposites: A Sequential Mesoscale Simulation Approach. RSC Adv. 2014, 4, 63586-63595.

(43) Deng, S.; Payyappilly, S. S.; Huang, Y.; Liu, H. Multiscale Simulation of Shear-Induced Mechanical Anisotropy of Binary Polymer Blends. RSC Adv. 2016, 6, 41734-41742.

(44) Qi, Y.; Liang, Y. Mesoscale Modeling of the Influence of Morphology on the Mechanical Properties of Proton Exchange Membranes. Polymer 2011, 52, 201-210.

(45) Chen, H.; Jiao, Y.; Liu, Y. Investigating the Microstructural Effect on Elastic and Fracture Behavior of Polycrystals Using a Nonlocal Lattice Particle Model. Mat. Sci. Eng. A-Struct. 2015, 631, 173-180.

(46) Chen, H.; Xu, Y.; Jiao, Y.; Liu, Y. A Novel Discrete Computational Tool for Microstructure-Sensitive Mechanical Analysis of Composite Materials. Mat. Sci. Eng. A-Struct. 2016, 659, 234-241.

(47) Monette, L.; Anderson, M. P. Elastic and Fracture Properties of the Two-Dimensional Triangular and Square Lattices. Modell. Simul. Mater. Sci. Eng. 1994, 2, 53-66.

- 40 -

(48) Wong, R. H. C.; Lin, P.; Tang, C. A. Experimental and Numerical Study on Splitting Failure of Brittle Solids Containing Single Pore Under Uniaxial Compression. Mech. Mater. 2006, 38, 142-159.

(49) Tang, C.; Zhang, Y.; Liang, Z.; Xu, T.; Tham, L.; Lindqvist, P. A.; Kou, S.; Liu, H. Y. Fracture Spacing in Layered Materials and Pattern Transition from Parallel to Polygonal Fractures. Phys. Rev. E 2006, 73, 056120.

(50) Wu, D.; Li, Y.; Zhang, J.; Chang, L.; Wu, D.; Fang, Z.; Shi, Y. Effects of the Number of Testing Specimens and the Estimation Methods on the Weibull Parameters of Solid Catalysts. Chem. Eng. Sci. 2001, 56, 7035-7044.

(51) Hallbauer, D. K.; Wagner, H.; Cook, N. G. W. Some Observations Concerning the Microscopic and Mechanical Behaviour of Quartzite Specimens in Stiff, Triaxial Compression Tests. Int. J. Rock Mech. Min. Sci. 1973, 10, 713-715.

(52) Tang, C.; Kou, S. Crack Propagation and Coalescence in Brittle Materials under Compression. Eng. Fract. Mech. 1998, 61, 311-324.

(53) Bruno, G.; Efremov, A. M.; Levandovskyi, A. N.; Clausen, B. Connecting the Macro-and Microstrain Responses in Technical Porous Ceramics: Modeling and Experimental Validations. J. Mater. Sci. 2011, 46, 161-173.

- 41 (54) Kourkoulis, S. K.; Markides, Ch. F.; Chatzistergos, P. E. The Brazilian Disc under Parabolically Varying Load: Theoretical and Experimental Study of the Displacement Field. Int. J. Solids Struct. 2012, 49, 959-972.

(55) Cleveland, J. J.; Bradt, R. C. Grain Size/Microcracking Relations for Pseudobrookite Oxides. J. Am. Ceram. Soc. 1978, 61, 478-481.

(56) Lu, C.; Danzer, R.; Fischer, F. D. Scaling of Fracture Strength in Zno: Effects of Pore/Grain-Size Interaction and Porosity. J. Eur. Ceram. Soc. 2004, 24, 3643-3651.

(57) Buxton, G. A.; Balazs, A.C. Simulating the Morphology and Mechanical Properties of Filled Diblock Copolymers. Phys. Rev. E 2003, 67, 031802.

(58) Buxton, G. A.; Balazs, A.C. Predicting the Mechanical Properties of Binary Blends of Immiscible Polymers. Interface Sci. 2003, 11, 175-186.

(59) Caiulo A.; Kachanov, M. On Absence of Quantitative Correlations Between Strength and Stiffness in Microcracking Materials. Int. J. Fract. 2010, 164, 155-158.

(60) Cannillo, V.; Leonelli, C.; Manfredini, T.; Montorsi, M.; Boccaccini, A. R. Computational Simulations for the Assessment of the Mechanical Properties of Glass with Controlled Porosity. J. Porous Mat. 2003, 10, 189-200.

(61) Cannillo, V.; Manfredini, T.; Montorsi, M.; Boccaccini, A. R. Use of Numerical Approaches to Predict Mechanical Properties of Brittle Bodies Containing Controlled

- 42 Porosity. J. Mater. Sci. 2004, 39, 4335-4337.

(62) Gibson, L. J.; Ashby, M. F. The Mechanics of Three-Dimensional Cellular Materials. Proc. R. Soc. Lond. 1982, 382, 43-59.

(63) Bruno, G.; Kachanov, M. Microstructure–Property Connections for Porous Ceramics: The Possibilities Offered by Micromechanics. J. Am. Ceram. Soc. 2016, 99, 3829–3852.

(64) Pandey, A.; Shyam, A.; Watkins, T. R.; Lara-Curzio, E.; Stafford, R. J.; Hemker, K. J. The Uniaxial Tensile Response of Porous and Microcracked Ceramic Materials. J. Am. Ceram. Soc. 2014, 97, 899-906.

(65) Tang, C. Numerical Simulation of Progressive Rock Failure and Associated Seismicity. Int. J. Rock. Mech. Min. Sci. 1997, 34, 249-261.

(66) Tang, C.; Xu, X.; Kou, S.; Lindqvist, P. A.; Liu, H. Numerical Investigation of Particle Breakage As Applied to Mechanical Crushing—Part I: Single-Particle Breakage. Int. J. Rock. Mech. Min. Sci. 2001, 38, 1147-1162.

(67) Bruno, G.; Kachanov, M. Porous Microcracked Ceramics under Compression: Micromechanical Model of Non-Linear Behavior. J. Eur. Ceram. Soc. 2013, 33, 2073-2085.

(68) Shyam, A.; Bruno, G.; Watkins, T. R.; Pandey, A.; Lara-Curzio, E.; Parish, C. M.;

- 43 Stafford, R. J. The Effect of Porosity and Microcracking on the Thermomechanical Properties of Cordierite. J. Eur. Ceram. Soc. 2015, 35, 4557-4566.

(69) Smolin, A. Y.; Roman, N. V.; Konovalenko, I. S.; Eremina, G. M.; Buyakova, S. P.; Psakhie, S. G. 3D Simulation of Dependence of Mechanical Properties of Porous Ceramics on Porosity. Eng. Fract. Mech. 2014, 130, 53-64.

(70) Cannillo, V.; Carter, W. C. Computation and Simulation of Reliability Parameters and Their Variations in Heterogeneous Materials. Acta Mater. 2000, 48, 3593-3605.

(71) Keleş, Ö.; García, R. E.; Bowman, K. J. Failure Variability in Porous Glasses: Stress Interactions, Crack Orientation, and Crack Size Distributions. J. Am. Ceram. Soc. 2014, 97, 3967-3972.

- 44 -

FIG.1. The scheme of 3D lattice spring model as a network consists of springs and nodes. Nearest {100} and next-nearest {110} neighboring interactions are considered.

FIG.2. Creation of a local fracture surface. The fracture surface is depicted in dark gray. Springs that would be attenuated on the formation of a fracture surface are colored black while neighboring springs that remain unbroken are in light gray. When a spring breaks, the two nodes linked by the spring will undergo damped vibration until the elastic energy stored in the spring is consumed by the surrounding viscous

- 45 fluid.

FIG.3. Scheme of the investigated system. The central part of the cylinder with 𝑟 = 32 arbitrary units is isolated for simulation of the Brazilian test. The load imposed on the edge is considered to be uniformly distributed.

- 46 -

FIG. 4. Morphology of the longitudinal section of cylindrical pellets containing pores with size fixed at 2. The porosity is (a) 10vol.%, (b) 20vol.%, (c) 30vol.%, (d) 40vol.% and (e) 50vol.% and a variation from a dispersed to a bicontinuous structure is observable.

- 47 -

FIG. 5. Normal strain fields of the longitudinal section of cylindrical pellets when compressed to a global strain of 5% in Z axis. The porosity is (a) 10vol.%, (b) 20vol.%, (c) 30vol.%, (d) 40vol.% and (e) 50vol.% and pore size is fixed at 2. The Young’s modulus of the solid phase is 5 times as great as that of pores.

FIG. 6. The probability distribution function of the local normal strains as a function of porosity.

- 48 -

FIG. 7. Normal stress fields of the longitudinal section of cylindrical pellets when compressed to a global strain of 5% in Z axis. The porosity is (a) 10vol.%, (b) 20vol.%, (c) 30vol.%, (d) 40vol.% and (e) 50vol.% and pore size is fixed at 2. The Young’s modulus of the solid phase is 5 times as great as that of the pores. The region marked by red circle indicates stress concentration caused by pore clustering.

FIG. 8. The probability distribution function of the local normal stresses as a function of porosity.

- 49 -

FIG. 9. Stress-strain plots for systems of different porosity and pore size. (a) The porosity of the plates ranges from 10~50vol.% with pore size fixed at 2 arbitrary units. (b) The pore size varies from 1 to 3 for plates whose porosity are 10vol.%, 30vol.% and 50vol.% respectively.

FIG. 10. Young’s modulus as a function of porosity and pore size. (a) EA=5EB. (b) EA =10EB. Error bars indicate one standard deviation from the mean.

- 50 -

FIG. 11. Number of element breakages as a function of strain. The porosity of the system is (a) 0vol.%, (b) 10vol.% and the pore size is 2 arbitrary units.

FIG. 12. Cumulative distribution function of the number of broken bonds as a function of strain. The porosity of the system is (a) 0vol.%, (b) 10vol.% and the pore size is 2 arbitrary units.

- 51 FIG. 13. Fracture stress as a function of porosity as well as pore size. Error bars indicate one standard deviation from the mean.

FIG. 14. Cumulative probability distributions of the crushing strength (in tensile) data on Weibull plot. (a) The porosity ranges from 10vol.% to 50vol.% and the pore size is 1 - 2 arbitrary units. (b) The pore size grows from 1 to 3 with porosity being 10 -20vol.%.

FIG. 15. Weibull modulus, m, as a function of porosity as well as pore size for cylindrical pellets under the Brazilian disc test. Error bars indicate the standard error

- 52 of the mean.