Statistical damage constitutive model based on self‐consistent Eshelby method for natural gas hydrate sediments

Natural gas hydrate (NGH) is widely distributed in marine sediments and continental permafrost, which is a promising energy resource. The sustainable and safe exploration and production of NGH requires fully understanding of its mechanical behaviors, which is still a challenge to our community. In this paper, a statistical damage constitutive model named SC‐SDCM is developed for NGH sediments. The NGH sediments are regarded as two phases: the matrix phase of solid mineral grains, pore fluid, and gas, and the inclusion phase of hydrate crystal. The effective elastic parameters of such two‐phase composite are estimated by self‐consistent method (SC) according to the equivalent inclusion principle of micromechanics. The mesoscopic element strength of the NGH sediments is described by the Weibull statistical distribution and the damage theory of composite materials. And then combined with the Drucker‐Prager failure criterion of microelement, the damage constitutive model of NGH sediments is established. The new SC‐SDCM is tested to be reliable and robust by triaxial experimental observations on artificial cores and naturally occurring samples under different confining pressures and hydrate saturations. The SC‐SDCM could properly describe the stiffness, peak strength, and strain softening properties of the NGH sediments. Moreover, the predictions of SC‐SDCM are closest to the experiment observations compared to several published constitutive models, especially for the near‐ and after‐damage stages with relatively high hydrate saturation.


| INTRODUCTION
Natural gas hydrate (NGH) is a kind of solid crystal composed of water and methane and widely gathered in land permafrost and marine sediments. [1][2][3] According to the US Geological Survey, the global reserve of NGH is estimated as about 100 000 to 300 000 000 trillion cubic feet, and the total amount of methane stored in NGH has far exceeded that in proved conventional natural gas reservoirs. [4][5][6] So NGH is an important and promising energy resource. However, the NGH reserves in marine sediments are usually in shallow, loose-compacted formations with low strength. Such reservoir condition leads to high risks of borehole instability and sanding during drilling and production, which bring great challenges to the safe and effective mining of NGH. For example, the Canadian Mallik 2L-38 well in 2007 and the Japanese trial production at Nankai Trough in 2013 were forced to shut off due to serious sand production after 60 hours and 6 days, respectively. [7][8][9][10] Moreover, the strength of hydrate reservoir will further decrease with the decomposition of NGH, which would result in the slip of hydrate sediments to the mining center and subsidence of the seabed, especially for the NGH production on the subsea slopes, and could even trigger geological disasters, such as submarine landslide and tsunami. 1,[11][12][13][14][15][16][17] Therefore, it is of great significance to well understand the mechanical properties of NGH sediments for the safe and sustainable exploitation of NGH.
Aiming to better understand the mechanical characteristics of NGH sediments, researchers have carried out lots of experimental tests, especially the triaxial shear tests. Due to the difficulties and high-costs for coring of the naturally occurring NGH sediments, most of the published triaxial tests are done using artificial cores with the skeleton made of Toyoura sand, [18][19][20] Ottawa sand, 21 or quartz sand, 19 and the NGH simulated by ice (eg, Masui et al 18 ) or synthetic methane hydrate (eg, Hyodo et al 22 ). These measurements on the artificial cores indicate that (i) the existence of gas hydrates increased the shear strength of hydrate sediments, and the cohesive strength and the elastic moduli of the sediments are influenced by the NGH saturation and the confining pressure. 18,19,[21][22][23] (ii) The stress-strain curve of artificial NGH cores shows strain stiffening before peak strength and strain softening after peak strength. 24,25 (iii) The type of sand that constitutes the skeleton of the artificial cores has a key effect on the stiffness of the sample, but has almost no effect on the strength of the specimen. 26 In addition, Yun et al 27 used sand soil, silt, and clay to synthesize soil specimens containing methane hydrate, and they observed that the stress-strain relation of hydrate sediments is complex function of sediment particle size, confining pressure and hydrate content. Kajiyama et al 20 used round glass beads and Toyoura sands to make artificial cores with hydrates, and then did triaxial tests, respectively. Though the glass beads and Toyoura sands have similar stiffness, these two types of artificial cores show quite different in their breakage pattern due to the shape of the grains and their roughness. Compared to numerous measurements on artificial cores, the experiments done with naturally occurring NGH samples are still far from enough. To our knowledge, Yoneda et al [28][29][30] has done series of measurements on samples from the seabed of Naikai Trough and the Krishna-Godavari Basin. The mechanical properties obtained by Yoneda et al [28][29][30] are in consistent with those of artificial specimens.
Based on experimental measurements, researchers have proposed different types of constitutive models to quantitatively describe the geomechanical behavior of NGH sediments, including nonlinear elastic model, 26 elastic-plastic model, 24,[31][32][33][34] damage statistical constitutive model 35,36 , and critical state model, 25,[37][38][39][40][41] which is an active area in recent years. There are mainly three key issues for the constitutive model of NGH sediments: the effects of saturation and cementation of NGH, the critical state of NGH, and the damage description. Klar 33 established an elastoplastic constitutive model within the state-dependent and shear dilatancy framework. The cementation strength parameter (p b ) was used to quantitatively describe the cementation contribution of NGH, and the degradation of the cementation was expressed as a function of the plastic deformation of sediments. Yan and Wei 34 introduced the concept of effective hydrate saturation to describe the influence of the occurrence mode of NGH, and the hydrate cementation effect was considered in a yield function. Another way to consider the influence of NGH on the mechanical properties of their sediments is to integrate the critical state of NGH phase behavior. Uchida et al 37 established a critical state model with five model parameters to describe the mechanical behavior of NGH sediments. Implicit (ISSF) and explicit analytical strain softening (ESSF) models were developed, respectively, by Pinkert and Grozic 42 and Pinkert et al, 43 which considered ten correlation coefficients for soil skeleton and five parameters related to NGH behavior. This kind of models was extended by introducing different parameters to consider the critical state line, cementation between hydrate and sand particles, and the like. [38][39][40][41]44 By analogy with the damage theory of permafrost, Wu et al 35 established a nonlinear damage constitutive model of NGH sediments based on the assumptions that the microelement strength followed Weibull distribution and controlled by the Drucker-Prager failure criterion. And some researchers also integrated models of the critical state soil mechanics to describe the inelastic mechanism. [45][46][47] Zhang et al 36 further developed the damage statistical constitutive model of NGH sediments to include the influence of damage threshold and residual strength. Although lots of efforts have made on the constitutive model of NGH sediments, it is still not fully solved and far from satisfactory due to its complexity.
Basically, NGH sediment is composed of solid mineral grains, hydrate crystals, pore fluid, and gas, which is a kind of multiphase composite material with remarkable heterogeneity in microstructure. Based on the idea of "homogenization" and integrated the information of the meso-structure, mesomechanics of composite materials offers methods to obtain the effective properties of macroscopic homogeneous and meso-heterogeneous media related to the Eshelby's inclusion problem. 48 Such effective inclusion methods mainly include self-consistent scheme, 49-52 differential method 53,54 and Mori-Tanaka method. 55 According to the idea of homogenization, researchers have proposed different methods to model the effective elastic properties of NGH sediments, such as Helgerud et al 56 and Nguyen et al. 57,58 In this paper, we develop a new constitutive model of NGH sediments based on the equivalent inclusion model of mesoscopic mechanics and the statistical damage theory. The NGH sediment is treated as two-phase composite material: solid mineral grains, pore fluid, and gas are the matrix phase (background material), and hydrate crystal is the inclusion phase. First, we will revisit the Eshelby's inclusion problem and predict the effective moduli of NGH sediments with different hydrate saturations by the selfconsistent scheme. Then, we will develop the framework of the constitutive model by integrating the statistical damage theory in mesomechanics and the self-consistent scheme. Finally, the performance of this new constitutive model will be tested by published experimental results and compared with other constitutive models.

OF NGH SEDIMENTS BY SELF-CONSISTENT SCHEME
The NGH sediments are shallow and loosely compacted under the seabed, which is a typical composite material with four components: solid mineral grains, hydrate crystals, pore fluid, and gas. From the perspective of mesomechanics of composite materials, the NGH sediment is a kind of medium with heterogeneous meso-structure and homogeneous at the macrolevel. Hence, the effective macroscopic properties of NGH sediments can be estimated with the properties of their components and the corresponding volume fractions and the geometric information by self-consistent method based on the theory of equivalent inclusions (as shown in Figure 1). The elastic response caused by an eigenstrain within its matrix is the so-called "inclusion problem," which is the basis of mesomechanics. The concept of eigenstrain 59 is referred to a small inelastic strain caused by some physical or chemical reasons (such as thermal expansion, phase transition, prestrain, or plastic strain) at a subdomain Ω within an infinite isotropic medium R (as shown in Figure 2). Such eigenstrain results in a self-balancing stress, such as thermal stress, phase transformation stress, prestress, and residual stress, which is the eigenstress.
Eshelby 48 decomposed the eigenstrain problem and transformed it to the problem with initial stress (-σ * ij ) in subdomain Ω and distribution force (p i ) on its boundary, which actually became an elastic mechanics problem of infinite elastic medium with a concentrated force in it. The distributed force p i on the surface can be regarded as a concentrated force in the infinite medium, and its displacement field can be expressed by the Kelvin solution with Green's function. For such classic inclusion problem, Kinoshita and Mura 60 gave the expression of uniform strain tensor in inhomogeneous matrix phase as follows: where G ikjl is the Green's function of the corresponding displacement tensor caused by unit force applied in a given direction, C klmn is the elastic stiffness tensor of the matrix phase, ε * mn is the eigenstrain tensor. The symmetric tensor of Green's function by Mura 59 is as follows: where K ip (x) = C ijpl x j x l is the Christoffel stiffness tensor in the x direction; x 1 = sinθcosθ/a 1 , x 2 = sinθcosθ/a 2 , x 3 = cosθ/a 3 , A transformation diagram of the mesoheterogeneous NGH sediments (gray as mineral particles, yellow as hydrate crystals, blue as pore fluids and gases) to macroscopic homogeneous media based on the concept of equivalent inclusion a 1 , a 2 , and a 3 are the semiaxes of the ellipsoidal inclusion. θ and φ are the spherical coordinates of the main axis vector x of the ellipsoidal inclusion. According to the Eshelby equivalent inclusion theory, each component of the microscopic heterogeneity structure in the composite material is regarded as an inclusion phase. And based on the volume average of Hooke's law, the effective elastic coefficient (C) of the composite material can be expressed as follows: According to Willis (1977), 61 the strain concentration factor A i (the ratio of inclusion strain to matrix strain) is defined as follows: where G is Green's function tensor, I is the fourth-order symmetric unit tensor, I ijkl = (δ ik δ jl + δ il δ jk )/2, and δ ik is the Kronecker symbol. Combining Equation (3) with Equation (4), the effective elastic coefficient (C) can be expressed as follows: where V i is the volume fraction of the ith inclusion; C i is the elastic coefficient of the ith inclusion. It should be noted that C appears on both sides of the equation and needs to be solved iteratively.
For the self-consistent (SC) approximation of mixtures with N phases in Equation (5), Berryman 62-64 gave a general form as follows: where V i , K i , and μ i are the volume fraction, the bulk modulus, and the shear modulus of the ith component, respectively. The geometric factors P and Q are functions of the effective aspect ratio (α, the ratio of the minor axis to the major axis for elliptical shape) of the components, which could be calculated by the T-matrix method given in the Appendix. The superscript *i on P and Q refers to the factors for the ith component in the background medium with self-consistent equivalent moduli K sc and μ sc . To highlight the contributions of NGH, we choose a two-phase model to describe the NGH sediments: the solid mineral grains, liquid, and gas together is one phase as the matrix, and the NGH is the other phase as inclusion. The volume fraction of NGH is determined according to the hydrate saturation (S h ) and porosity.
The two formulas in Equation (6) are coupled, and therefore, K sc and μ sc must be solved by simultaneous iteration. There are only two independent elastic constants for isotropic materials according to the elastic mechanics theory. Therefore, the effective Young's modulus (E sc ) and Poisson's ratio (v sc ) can be calculated once the effective bulk modulus and shear modulus of the NGH sediments are determined by the SC scheme in Equation (6)

| Meso statistical damage constitutive model of NGH sediments
It is assumed that the NGH sediments without considering damage are isotropic and linear elastic, so the stress tensor and strain tensor conform to Hooke's law as follows: where σ ij and ε ij are the stress tensor and strain tensor, respectively; λ is the Lame constant equal to (K sc -2μ sc /3); δ ij is Diagram for the Eshelby inclusion problem the Kronecker symbol, δ ij = 1 when i = j, and δ ij = 0 when i ≠ j. For conventional "drained" triaxial test, the pore pressure can be ignored, so the maximum principal stress (σ 1 ) is the axial compression pressure applied at the two ends of the core, and the middle principal stress is equal to the minimum principal stress (σ 2 = σ 3 ), which are both equal to the confining pressure (P c ). According to the stressstrain relationship in the principal stress space (Equation (9)), the triaxial stress-strain relationship of NGH sediments can be written as follows: where ε 1 is the axial strain; σ 3 is the confining pressure; Δσ = σ 1 -σ 3 is the axial deviation stress; E sc and v sc are the effective Young's modulus and Poisson's ratio, respectively, which are estimated by the SC scheme in section 2.
According to the equivalent strain hypothesis of Lemaitre (1984), 65 the strain ε caused by the nominal stress σ applied on the damaged material is equal to that caused by the effective stress σ * applied on the nondamaged material, that is: where D is the damage variable; I is unit tensor; C * and C are the elastic moduli of the damaged and nondamaged material, respectively. The elastic moduli C can use its selfconsistent estimations. So when describing the mechanical behavior after damage, the equivalent stress after damage can be used to replace the nominal stress in the constitutive relation of nondamaged materials. Based on the theory of damage mechanics, the constitutive relation of NGH sediments considering damage can be obtained based on Equations (10) and (11) as follows: It should be noted that E sc and v sc here are same to those in Equation (10) before damage. Similar to other naturally occurring materials (eg, rocks), the mechanical properties of NGH sediments also show obvious randomness, which is difficult to be well described with a single variable or parameter. Therefore, we discrete the NGH sediments to infinite microelements, the strength of which is a random variable (F). Considering that the Weibull distribution 66 can cover a variety of random distributions by changing its shape parameter (m) and has been widely used in different engineering areas, here we choose the Weibull distribution to describe the statistical distribution of F for NGH sediments. The probability density function of F with Weibull distribution is as follows: where F is the random variable describing the Weibull distribution of the element strength; m and F 0 are parameters of Weibull distribution.
Under the action of external temperature and pressure, the damage of NGH sediments grows gradually with the generation of microdefects within it. Assuming that the damage of NGH sediments is a continuous process and the damage variable D is a macrodescription of the cumulative effect of the microdefects, then the damage evolution equation can be obtained by integrating the probability density function within the interval of [0, F] as follows: Physically, the damage variable D is the fraction of the damaged microelements among the total number of elements (N t ), which varies between 0 and 1. D = 0 represents no damage or the initial state, and D = 1 represents the complete damage state, and 0 < D < 1 when partial damage occurs.
By substituting Equation (14) into Equation (12), we obtain: This is a general form of the constitutive model of NGH sediments considering damage. When there is no damage at relatively low load condition (D = 0), it reduces to the classic constitutive model of linear elasticity as Equation (10). When damage occurs (0 < D ≤ 1), additional formula is required to determine F, F 0 , and m, which will be deduced in the following section based on the Drucker-Prager criterion.

| Parameters of constitutive model for NGH sediments
Drucker and Prager 67 proposed the Drucker-Prager criterion that can reflect the influence of average principal stress (or the first stress tensor invariant I 1 ) on the shear strength of geotechnical materials, which is an extended form of Mohr-Colomb failure criterion and has been widely used. [68][69][70][71] Here, for the microelements of NGH sediments with the Weibull-distributed strength of F, its Drucker-Prager failure criterion be expressed as follows: where I 1 is the first stress tensor invariant, J 2 is the second deviation stress invariant, α 0 is the material parameter associating Drucker-Prager and Mohr-Coulomb criteria in the main stress space, ϕ is the internal friction angle. σ 1 * , σ 2 * , and σ 3 * are the effective stresses. σ 1 , σ 2 , and σ 3 are nominal stress.
According to the equivalent strain hypothesis (Equation (11)), we have: According to the generalized Hooke's law, the axial strain can be expressed as follows: So Equation (17) can be rewritten as follows: In the conventional triaxial compression test, confining pressure σ 2 =σ 3 , substituting Equation (19) into Equation (18), F can be expressed as follows: As a general observation in the conventional triaxial compression test of NGH sediments that the measured stress-strain variation trend changes at the peak point which is mathematically a stationary point for the stressstrain relation. Therefore, in Equation (15), the peak of deviation stress (Δσ max ) and the corresponding peak axial strain (ε max ) under specific confining pressure and hydrate saturation should meet: Based on Equations (15), (20), and (21), the parameters m and F 0 can be expressed as follows: and where F max is the microelement strength of hydrate sediment corresponding to the peak of the deviatoric stress.
The Equations (6), (15), (20), (22), and (23) together give the constitutive model of NGH sediments considering the damage effect. We name this constitutive model as SC-SDCM in the following section for convenience. The model contains five parameters, which are Young's modulus E sc , Poisson's ratio v sc , the microelements strength of NGH sediments F, and Weibull distribution parameters F 0 and m. The elastic parameters (E sc and v sc ) of NGH sediments linked to the hydrate content are calculated by the self-consistent method, the Weibull distribution parameters F 0 and m are determined by the peak points of triaxial test curve, and F is calculated by Weibull distribution function.

| MODEL EVALUATION
To evaluate the model performance, we compare the predictions of our SC-SDCM with published experimental data and other models, respectively. The main procedures are as the followings. Such procedures and details are also shown in Figure 3.
(i) Identify the matrix properties (E 1 , v 1 ) of NGH sediments according to the measurements of cores at different confining pressures when the hydrate saturation is zero. (ii) Estimate the effective elastic moduli (E sc , v sc ) of NGH sediments with different saturations (S h ) by SC scheme according to Equation (6). (iii) Based on the results of (ii) and the peak point (ε max , Δσ max ) of the measured stress-stain curves, calculate the parameters F 0 , m of the SC-SDCM according to Equations (22) and (23). (iv) Based on the results of (ii) and (iii), calculate the corresponding deviation stresses at different axial strains under different hydrate saturations (S h ) and confining pressures (P c ) according to the SC-SDCM and then compare the results of SC-SDCM with experimental data and other models.

| Comparison with experimental data
For the comparison with experimental data, we choose representative measurements on artificial cores by Masui et al 18 (Figures 4 and 5). To quantify the reliability of the SC-SDCM, we calculated the relative error of the predictions to the experiment measurements. The statistical distributions of relative errors are shown in Table B1 (see in Appendix B). For all the data in Figures 4 and 5, the mean values of the relative errors are 7.3% and 6.8% (Table B1), respectively. The predictions of SC-SDCM with relative error smaller than 5% make up to about 67.8% of all the measurements. To obtain cores with different hydrate saturations, the methane gas was injected into the skeletons at different confining pressures with temperature below −5.15℃ for 24 hours. After that, triaxial tests were carried out at pore pressure of 8 MPa and confining pressure of 8.5 MPa, 9 MPa, 10 MPa, and 11 MPa, respectively. And the axial load was applied at the stain rate of 0.1% per minute. According to effective stress principle of porous media, the effective confining pressure for these measurements are correspondingly 0.5 MPa, 1 MPa, 2 MPa, and 3 MPa, respectively. The effective elastic moduli and Weibull distribution parameters for the artificial NGH cores with different hydrate saturations by Miyazaki et al 26 are shown in Table 3.
The comparison results between the predictions of SC-SDCM and experiments observations of Miyazaki et al 26 using Toyoura sand, No.7 silica sand, and No.8 silica sand are shown in Figures 6-8, respectively. In these figures, the measurements by Miyazaki et al 26 are shown with mark points, while the predictions of SC-SDCM are shown by solid curves. The results indicate that the SC-SDCM works well regardless of the sand type of the core skeleton. All the stress-strain curves (Δσ -ε 1 ) at different confining pressures and hydrate saturations are properly predicted by SC-SDCM with robust reliability. The relative errors were also calculated (See Table B1 in Appendix B). The performances of the SC-SDCM for the Toyoura sand cores ( Figure 6) and the No.7 sand cores ( Figure 7) are better than that for the No.8 sand cores (Figure 8). For NGH cores with skeleton of No.8 sand, the larger deviations between predictions of SC-SDCM and experimental measurements are mainly at relative low hydrate saturations (eg, S h = 0% and 27%-28%) when the effective confining T A B L E 2 Effective elastic moduli and Weibull distribution parameters for the artificial specimens in Masui et al (2005) 18  pressure is 3MPa ( Figure 8A, B). The average of the relative errors for all the data is about 10% (see details in Table B1 in the Appendix B).

| Comparison with Yoneda et al (2019)
Yoneda et al 30 Figure 9 where the measurements by Yoneda et al are shown with mark points, and the predictions of SC-SDCM are shown by solid curves. Again, the SC-SDCM gives satisfactory predictions for these naturally occurring samples as it works for the artificial cores. The mean value of the relative errors is only 8.1% (see details in Table B1 in the Appendix B), which is a quite good result considering the complexity of naturally occurring materials.

| Comparison with other published constitutive models
Besides the verifications of the SC-SDCM by experimental data above, we also compared its performance with other published models. We made two sets of comparisons based on the experiment observations of Masui et al 18 and Miyazaki et al, 26 respectively. That is: (i) Using the experiment measurements of Masui et al 18 at different hydrate saturations when the confining pressure is 1MPa, the predictions of our SC-SDCM and other four constitutive models, including Uchida et al, 37 Yan and Wei, 34 26 Pinkert and Grozic 42, and Uchida et al, 40 are compared. The results are shown in Figure 11.
These comparisons indicate that the SC-SDCM uses relatively fewer model parameters with clear physical meanings and give similar or even better predictions than other models. As can be seen in Figures 10 and 11, all these constitutive models could basically capture the variation trends of the stress-strain curves (Δσ -ε 1 ) at different hydrate saturations. Compared to other models, the predictions of SC-SDCM are closest to the experiment

| CONCLUSION
In this paper, a new constitutive model, SC-SDCM, is developed for natural gas hydrate (NGH) sediments. The gas hydrate sediment is considered as a composite composed of matrix phase and inclusion phase. The matrix phase is conceptually equivalent elastic "solid" made of solid mineral grains, pore fluid, and gas with their volume fractions similar to the real state of NGH sediments, while the inclusion phase is the hydrate crystal. Under such work frame, the effective elastic parameters of the hydrate sediment are estimated by self-consistent method according to the equivalent inclusion theory of micromechanics. The mesoscopic element strength of the gas hydrate sediment is described by the Weibull statistical distribution and the damage theory of composite materials. And then combined with the Drucker-Prager failure criterion of microelement, the damage constitutive model of NGH sediments (SC-SDCM) was established. The SC-SDCM considers the influence of mineral composition, hydrate content, and effective pressure on the mechanical properties of NGH sediments, but it has not included the effect of other factors, such as pore structure, hydrate type, and hydrate distribution. The new SC-SDCM is tested to be reliable and robust within a wide range of hydrate saturations and stress conditions. The predictions of SC-SDCM are consistent with the triaxial experimental observations on artificial cores by Masui et al 18  The SC-SDCM established by the concept and method of composite materials has clear physical meaning and requires fewer model parameters which are easy to obtain based on experimental data or geophysical logging data. So, this new model is convenient for practical engineering application. Meanwhile, based on Eshelby equivalent inclusion theory, the relationship between the mesocomposition of natural gas hydrate sediments and their macroscopic mechanical properties is established by the SC-SDCM model. With the further research on the mesocomposition and structure of hydrate sediments, more factors effecting the mechanical properties of hydrate sediments will be understood better, such as pore structure, hydrate type, and hydrate distribution. Their effect on the mechanical properties can be involved in the estimation of the effective moduli of hydrate sediments. Then, the SC-SDCM model can be optimized and expanded to describe the mechanical properties of hydrate sediments better. F I G U R E 1 1 Comparison between SC-SDCM and three published models of Miyazaki et al, 26 Pinkert and Grozic (2014) 42 , and Uchida et al (2016) 40 based on the measurements by Miyazaki et al 26 using artificial NGH cores. The experimental data were measured at the hydrate saturations of 0%, 27%-34%, and 41%-45% when the confining pressure was 1 MPa, which are indicated by circle points here

APPENDIX A
The values of P and Q for ellipsoidal inclusions with arbitrary aspect ratio (α) are expressed as follows: Tensor T ijkl links up the strain in uniform far-field and the strain in the ellipsoidal inclusions (Wu, 1966). 63 Berrymann (1980) gave the calculation of relative scalars required for P and Q: where, with A, B, and R given by: where K i and i are the volume modulus and shear modulus of the ith phase inclusion, and v m is the Poisson's ratio of the matrix. The functions β and f are given as follows: