Three-dimensional fatigue crack propagation simulation using extended finite element methods for steel grades S355 and S690 considering mean stress effects

The assessment of fatigue crack propagation of steel structures is essential and important especially to improve the application of high strength steel in construction. The load ratio R, reflecting mean stress effects, will be changed with crack extension in the steel structures with complicated geometry. In this paper, the Walker equation is employed to fit the fatigue crack propagation rate of steel grades S355 and S690 based on experimental data in the literature to incorporate the mean stress effects. The material fatigue crack propagation pa- rameters with 95%, 97.7%, and 99% guarantee of Walker equation were obtained by a stochastic analysis using the Monte Carlo method. The fatigue life was firstly predicted by the analytical method and was used as a baseline for numerical fatigue crack propagation simulation. A user-defined fatigue crack propagation subroutine based on the Walker equation was developed using phantom nodes-based extended finite element method (PN- XFEM) and Virtual Crack Closure Technique (VCCT) to consider the mean stress effects. The proposed three-dimensional fatigue crack propagation simulation subroutine is successfully validated of both steel grades, S355 and S690.


Introduction
The application of high strength steel in construction is becoming increasingly attractive during the past two decades, especially in infrastructure applications. High strength steels provide the potential for economical solutions for highly loaded and slender structures. Despite the benefits of the increased yield strength, the fatigue resistance of steel structures is one of the most important concerns [1]. The fatigue phenomenon is the process of progressive localized permanent structural change occurring in the material under cyclic loading. The fatigue process of steel structures is divided into two stages: the fatigue crack initiation period and the fatigue crack propagation period [2][3][4].
In 1961, from the fatigue experiments, Paris [5] proposed a fatigue crack growth rate equation that links the stress intensity factor (SIF) range ΔK to sub-critical fatigue crack growth rate dα/dN, known as "Paris' law". Nowadays, the Paris' law and its extensions [6][7][8][9][10] are widely used to predict the fatigue crack growth for different engineering structures. The crack growth curve is divided into three regions of generalized behaviour. At the lower end, the fatigue crack growth rate dα/dN approached zero at a threshold ΔK th , below which the crack will not grow. At the intermediate part, the crack growth curve is linear in the log-log plot, described by Paris' equation. The unstable region at the top end is characterized by rapid, unstable crack growth, which is due to two possible reasons [11]: The crack growth accelerated as maximum SIF approaches fracture toughness; ② The influence of crack tip plasticity on the true driving force for fatigue is larger at higher SIF range.
One shortcoming of the original Paris' equation is that it could not consider mean stress effects. In 1970, Elber [12] noticed that the compliance of fatigue specimens at higher loads agreed with the standard formulae for fracture mechanics specimens but at low loads, the compliance was close to that of an un-cracked specimen. Elber [12] postulated that this change in compliance was due to the contact between crack surfaces at loads that were low but greater than zero, known as crack closure. Crack closure decreased the fatigue crack growth rate by reducing the effective stress intensity range in correlating fatigue crack growth data at various R ratios. Elber [12] proposed an empirical relationship of effective stress intensity ratio U by measuring the closure stress intensity in aluminum at various load levels and R ratios. Similar empirical expressions were reported by subsequent researchers [13,14]. Walker [15] reported on the effects of the R-ratio on crack propagation for aluminum in 1970, concluding that increasing ratio R resulted in an increased growth rate. An empirical model with three material parameters was developed to address the effect of R-ratio on fatigue crack growth rate, denoted as the Walker equation. The Paris and Walker equations did not address the asymptotic behaviour in the unstable region. A theoretical model to take into account the crack closure effects based on the Walker model and including the threshold stress intensity factor range, ΔK th , was proposed by Correia et al. [16]. Forman et al. [7] proposed a model incorporating both the R-ratio and maximum SIF effects to address the above behaviour. In the Forman equation, the crack propagation rate dα/dN tends toward infinity when the maximum SIF approaches fracture toughness. Besides, the crack growth rate could also be affected by other parameters [17], such as load frequency, load sequencing, temperature, and environmental factors, etc. Hartman and Schijve reported [18] that an increase in loading frequency produces a decrease in growth rates of aluminum alloys, and the fatigue crack growth rates are quite different when it is exposed to the different aging environments. Yokobori and Kiypshi [19] investigated the influence of temperatures and stress intensity factors effects on fatigue crack propagation rate and striation spacing exposed to high temperature from 1500 K to 3500 K. Lesiuk et al. [20] have investigated the influence of heat treatment on the fatigue crack propagation for the 42CrMo4 and 41Cr4 steels based on experimental testing. Several experimental tests have been conducted to study the fatigue crack propagation behaviour in structural steels under pure-mode I and mixed-mode (I + II) loading conditions applied to old puddled irons [21], S355J0 steel [22], and S235 steel [23]. Silva et al. [23] investigated the crack propagation rates under mixed-mode (I + II) loading conditions based on an analytical-tonumerical methodology using the digital image correlation (DIC) data and finite element modelling. The fatigue crack propagation phase in the fatigue life prediction of structural details has been done based on Paris material properties obtained experimentally and supported by numerical simulation through the finite element method combined with the Virtual Crack Closure Technique (VCCT) or J-integral method [24][25][26].
Several methodologies could be used for the fatigue crack propagation simulation, such as traditional finite element methods [27], the extended element method (XFEM) [28], mesh-free method [29], discrete element method [30][31][32], Peridynamics [33], and phase-field approach [34]. Da Silva et al. [35] studied the fatigue life prediction of welded joints where a comparison between the finite element method  (FEM) and the extended element method (XFEM) applied to the fatigue crack propagation stage was presented. Some studies of 2D and 3D crack propagation modelling using the XFEM method have been presented [36][37][38][39][40][41], however, without considering the R-ratio effects. Gonzalez-Herrera et al. [42] have analyzed the crack opening and closure in 2D and 3D fatigue crack propagation modelling applied to compact tension specimens (CT) considering the meshing and element size effects. Mahmoud and Riveros [43] suggested a fatigue crack propagation analysis of a single stiffened ship hull panel based on a finite element modelling combined with a Monte Carlo simulation (MCS). Lee et al. [44] analyzed the fatigue behaviour of a welded T-joint considering the crack initiation and propagation phases supported by finite element simulation where the authors found that the residual stress effects cannot be neglected. Huang et al. [45] developed a numerical study on fatigue crack propagation behaviour of CFRP-repaired FPB specimens under mixed-mode loading conditions considering the NASGRO model. Menk et al. [46] predict the lifetime of solder joints undergoing thermal cycling for the electronic industry based on microstructural features of the joint using XFEM. In addition, the developments of accuracy and convergence of the crack driving force are very important [47]. Wyart et al. [48,49] investigated two sub-structuring methods to allow for the application of XFEM within commercial FE software without the need for modifying their kernel. Duflot and Bordas [50,51] derived an associated recoverybased error indicator based on global recovery techniques. Pierard et al. [52] evaluate the influence of the contact around crack lips on fatigue crack path and crack growth rate under multi-axial loading conditions. The results showed that both fatigue crack path and growth rate are strongly influenced by the partial contact around the crack tip. Peng et al. [53] proposed an isogeometric boundary element method based on non-uniform rational B-splines to simulate three-dimensional crack growth. Convergence studies in the crack opening displacement are performed for a penny-shaped crack and an elliptical crack. Sutula et al. [54][55][56] investigated the energy-minimal multiple crack propagation in a linear elastic solid under quasi-static conditions. The crack extension directions and lengths are calculated based on the minimization of the total energy of the mechanical system through three parts: Theory and state of the art review, discrete solution with XFEM, and XFEM computer implementation and applications.
In order to evaluate the fatigue performance of high strength steels, the fatigue crack propagation tests of steel grades S355 and S690 were reported by Jesus et al. [1], and four different groups of crack propagation parameters were fitted based on the Paris equation in terms of different loading ratio R. Besides, the commercial software ABAQUS [57] includes the extended finite element method (XFEM) [28] to predict the fatigue crack growth using original Paris' law based on Griffith energy rate (G) and the Virtual Crack Closure Techniques (VCCT). In this way, the load ratio effects on fatigue crack propagation properties of structural steels are not considered currently. However, in the engineering application, the stress ratio, R, will be changed with crack extension in terms of steel structures with complicated geometry. The crack propagation model in the commercial software ABAQUS is expected to be extended, and the material parameters of steel grades S355 and S690 are suggested to be fitted to incorporate the stress R-ratio effects to tackle the shortcoming of original Paris' equation. Noted that this paper is limited to the linear elastic fracture mechanics, the enrichment shape functions should be enhanced to represent the singularities if elastic-plastic fracture mechanics are expected. Elguedj et al. [58] proposed a new tip enrichment basis coupled with a Newton-like iteration scheme and a radial return method for plastic flow. The results presented good accuracy for numerical evaluation of standard fracture parameters such as J-integral. Elguedj et al. also [59] presented an augmented Lagrangian formulation in the XFEM framework that can deal with elastic-plastic crack growth with the treatment of contact. Pañeda et al. [60] presented a gradient enhanced framework by means of XFEM using microstructurally-motivated models.
Hence, in this paper, the Walker equation is employed to fit the fatigue crack propagation rate of steel grades, S355, and S690, based on experimental data in the literature to incorporate the R-ratio effects. The parameters fitted by normal probability distributions, with 95%, 97.7%, and 99% guarantee of the Walker equation were obtained by a stochastic analysis using Monte Carlo simulation. The fatigue life was firstly predicted by the analytical method and was used as a baseline for the numerical fatigue crack propagation model. A user-defined fatigue crack propagation subroutine was developed based on the Walker equation    using phantom nodes-based extended finite element method (PN-XFEM) and VCCT to consider the mean stress effects. The proposed threedimensional fatigue crack propagation simulation method is successfully validated of both steel grades, S355 and S690.

Fatigue crack propagation considering r-ratio effects
Jesus et al. [1] performed crack propagation tests of steel grades S355 and S690 with four different load ratios, namely R = 0.01, R = 0.25, R = 0.5, and R = 0.75. Fig. 1 illustrates the detailed geometry of the compact tension (CT) specimen [61]. The SIF ranges were computed using the formulation proposed by ASTM E647 for the CT specimens [61], as shown in Eq.(1). The seven-points incremental polynomial technique, adopted by ASTM E647 [61], was used to determine the crack propagation rate.
where: ΔF is the applied load range, a is crack size, W is the effective width, B is the thickness. As shown in Eqs. (2) and (3), the Walker equation [15] proposed the equivalent SIF, ΔK, to consider the load ratio effects. As shown in Fig. 2, a linear regression was used in MATLAB [62] to fit the material parameters in the Walker equation after transforming the data to log-log space. The fitted coefficients of the Walker equation for steel grades S355 and S690 are listed in Table 1. The Walker equation collapses the data from R-ratios to a single line when plotting crack propagation rate dα/dN versus the equivalent SIF, ΔK, as shown in Fig. 3. It is noted that the Root Mean Square Error (RMSE) of steel grade S690 is smaller than S355.
where:C 0 , m, γ are material parameters of fatigue crack propagation in the Walker equation.  Fatigue crack growth is inherently stochastic to predict. The micro material imperfections, such as voids and secondary particles, have a statistical distribution, and it is important to investigate how these imperfections affect the structure's fatigue life in a presence of a crack to guarantee a safer and more robust design. Hence, the Monte Carlo method [63] is employed to investigate the fatigue life distribution, where the material properties are replaced by Probability Distribution Functions (PDF), random values are selected from each PDF and the analysis is run as many times as needed.
Three material constants are presented in the Walker equations and each is assumed to be normally distributed around the mean value in log-log space where the coefficients were fit, shown in Fig. 2. But the crack growth constants are not suggested to be independently selected at random from one PDF because those parameters are jointly distributed during linear regression fit [64]. A multivariate normal distribution, expressed in Eq. (4), is used as a PDF through MATLAB [62] software. A total of one million random material parameters was selected based on multivariate normal distribution using the mean value and standard derivation listed in Table 1 to meet the large numbers requirement and central limit theorem [63]. The histogram of material parameters for steel grades S355 and S690 is presented in Figs. 4 and 5. Because the random parameters from multivariate normal distribution will interact with each other, the mean value of random material parameters agreed well with the fitted results, but the standard derivation is different from the original fitted results in Table 1.
where: μ is the mean value, Σ is the non-singular symmetric covariance matrix describing the interaction between constants, n is the total number of samples. The fatigue crack growth rate is calculated based on one million random material parameters generated by a multivariate normal distribution. The fatigue crack growth rate distribution of steel S355 and S690 based on Monte Carlo Simulation is presented in Figs. 6 and 7 in terms of different loading ratios. The fatigue crack growth rate with 95%, 97.7%, and 99% guarantee rate is calculated with four different load ratios, presented in Figs. 6 and 7. In terms of different guarantee rate, the material parameters in the Walker equation are fitted again using linear regression by MATLAB [62] after transforming the data to log-log space. The material coefficients of the Walker equation based on stochastic analysis are listed in Table 2. The results showed that a higher guarantee rate leads to a larger value of material parameter C 0 and m, but the effect on material parameter γ is not obvious.

Analytical fatigue crack propagation prediction
The fatigue crack propagation rate could be described as Eq. (5) by combining Eq. (1) and Eq. (2). The fatigue crack propagation life could be calculated analytically by three methods, namely integration method, crack increment iteration method, and fatigue cycle increment iteration method.
(1) Integration method The fatigue life is calculated by solving the integral equations with the crack propagation rate listed in Eq. (6). In this paper, the integral equation is solved numerically by MATLAB [62] based on Eq. (6).

(2) Crack increment iteration method
The fatigue life is predicted iteratively, as expressed in Eq. (7), with a constant crack length increment Δa. The fatigue life increment is calculated by the forward Euler method listed in Eq. (8).
(3) Fatigue cycle iteration method The fatigue crack propagation length is predicted iteratively, as expressed in Eq. (9), with a constant fatigue life increment ΔN. The crack length increment is calculated by the forward Euler method listed in Eq. (10).
In terms of steel grade S355 with an initial crack length, a i = 11.0 mm, exposed a constant load range ΔP = 5.7 kN, the fatigue crack growth comparison with different increment size is shown in Fig. 8. The results showed that, compared with integrated value, the larger crack length increment leads to an overestimation of the fatigue cycle in terms of crack increment iteration method, but the larger fatigue cycle increment leads to an underestimation of the crack length in terms of fatigue cycle iteration methods. Tables 3 and 4 listed the error comparisons using crack increment iteration methods and fatigue cycle iteration methods. The results showed that the fatigue cycle difference at a f = 25 mm between crack increment iteration methods and integrated value is within 5% when the crack length increment is smaller than 0.3 mm and is within 10% when the crack length increment is smaller than 0.6 mm. The results also showed that the crack length difference at N f = 5.7 × 10 6 between the fatigue cycle increment method and integrated value is within 5% when the fatigue cycle increment is less than 3000 and is within 10% when the fatigue cycle increment is less than 7000.
For CT specimen with an initial crack a i = 15 mm exposed a constant load range ΔP = 5.7 kN, the fatigue life is predicted using crack increment iteration method with increment Δa = 0.001 mm based on material parameters with different guarantee rate. The fatigue crack propagation comparisons with different guarantee rates are shown in Fig. 9. The results showed that the guarantee rates affect the fatigue crack propagation a lot, and higher guarantee rates contribute to a much faster crack propagation rate. Interestingly, the fatigue cycle ratio at the same crack length between different probabilities kept as a constant when the crack propagated from 15 mm to 30 mm. The probability effects on fatigue crack propagation are quantified in Table 5. The fatigue life calculated by parameters with 95%, 97.7%, and 99% guarantee rates is only 0.03, 0.008 and 0.004 times of it predicted from the parameters fitted by experiments for steel grades S355 respectively, and the fatigue life calculated by parameters with 95%, 97.7%, and 99% guarantee rates is 0.055, 0.030 and 0.016 times of it predicted from the parameters fitted by experiments for steel grades S690 respectively.
For CT specimen made of S355 steel with an initial crack a i = 15 mm exposed a constant load range ΔP = 5.7 kN, the fatigue life is predicted using crack increment iteration method with increment Δa = 0.001 mm with different R ratio. The load ratio R effects on fatigue crack propagation are shown in Fig. 10. The results showed that a higher R ratio exposed to a constant load range leads to faster crack propagation. The R ratio effects on fatigue crack propagation exposed to a constant load range are quantified in Table 5. The fatigue cycle ratio exposed to a constant load range at the same crack length between different R ratios kept as a constant when the crack propagated from 15 mm to 30 mm.

Implementation of fatigue crack propagation model considering mean stress effects
Modeling discontinuities as an enriched feature [28] is an effective way to simulate the initiation and propagation of the discrete crack along an arbitrary, solution-dependent path without the requirement of re-meshing in the bulk materials. Based on the concept of partition of unity, XFEM [28] treated the cracks as a special enriched function in conjunction with additional degrees of freedom. The nodes are enriched with the jump function when the elements are intersected by a crack.
The jump function H(x) for a crack is given by: As illustrated in Fig. 11, phantom nodes are superposed on the original real nodes aiming to represent the discontinuity of the cracked elements. The element is "cracked" by splitting it into two parts using the level set method (LSM) [65], where two orthogonal signed distance functions are defined, one is used to describe the crack surface and the other is to build the orthogonal surface. The intersection of those two surfaces gives the crack front.
The VCCT criterion [66], under the assumption that the strain energy released when a crack is extended by a certain amount is the same as the energy required to close the crack by the same amount, is used to calculate the strain energy release rate. As shown in Fig. 12, the strain energy release rate of Mode I [57] could be calculated through the following equations under the assumption of linear elastic behavior. Similar arguments and equations could be used to calculate the strain energy release rate of Mode II and III.
where: G I is the Mode I energy release rate, b is the width, d is the length of the elements at the crack front, F v,2,5 is the vertical force between nodes 2 and 5, v 1,6 is the vertical displacement between nodes 1 and 6. A quasi-static analysis subjected to sub-critical cyclic loading is used to simulate the fatigue crack propagation. The Walker equation is implemented through user subroutine UMIXMODEFATIGUE based on commercial software ABAQUS [57] to predict the fatigue crack propagation considering mean stress effects. The detailed implementation flow chart is shown in Fig. 13. The crack length a N is extended, from the current cycle forward over a number of cycles ΔN, to a N+ΔN by fracturing at least one enriched element ahead of the crack tips. Based on the Walker equation combined with the known element length and propagation direction at the enriched elements ahead of the crack tips, the number of cycles necessary to fail each of the enriched elements ahead of the crack tips could be calculated. The minimum number ΔN is represented as the number of cycles to grow the crack equal to its element length, as the cycle jump strategies. The element is completely fractured with a zero constraint and zero stiffness at the cracked surfaces at the end of the completed cycle. As the enriched element is fractured, the load is redistributed, and a new fracture energy release rate will be calculated for the enriched elements ahead of the crack tips for the next cycle. To accelerate the fatigue crack growth analysis and to provide a smooth solution for the crack front, a nonzero tolerance,ΔD tol , for the least number of cycles to fracture an enriched element, as presented in Eq. (14).

Three-dimensional fatigue crack propagation using XFEM method
This section will investigate the effects of boundary condition, mesh size, damage tolerance, and R-ratio on fatigue crack propagation of CT specimens based on the XFEM model using commercial FEM software ABAQUS. One cycle in the numerical simulation is assumed as 1 s, and the time increment is defined as 0.01 s. The material is assumed as elastic, with a modulus 210.0 GPa and the Poisson's ratio 0.3. The  Fig. 11. The schematic of the phantom node method [2,3]. element type of CT model is C3D8. To accelerate the fatigue crack growth analysis and to provide a smooth solution for the crack front, a nonzero tolerance, ΔD tol (expressed in Eq. (14)), for the least number of cycles to fracture an enriched element. The nonzero tolerance,ΔD tol , is assumed to be 0.1.

Boundary condition effects
As shown in Fig. 14, the fatigue crack propagation is predicted by three-dimensional XFEM methods with five different boundary conditions to investigate the effects of the area tied to reference point, ② constrains of horizontal displacement, ③ constrains of rotation. The boundary condition effects on fatigue crack propagation of steel grade S355 is shown in Fig. 15. The results of XFEM simulation with BC3, BC4, and BC5 agreed well with the integration methods, but the simulation results with BC1 and BC2 lead to an overestimation of fatigue crack propagation life because the loading point is no longer pin support due to over constrains. In all the boundary conditions, BC3 agreed best with the analytical results.

Mesh size effects
The crack length is extended from the current cycle forward over a number of cycles by fracturing at least one enriched element ahead of the crack tips. In order to alleviate numerical error induced by the forward Euler method, similar to the crack increment size effects reported in Fig. 8, the mesh size is suggested to smaller than a certain value. As shown in Fig. 16, fatigue crack propagation simulation is conducted with different mesh sizes and boundary conditions of BC3 and load ratio R = 0.01. The mesh size effect on fatigue crack propagation is shown in Fig. 18. The crack shape at 27 mm of CT specimen with different mesh size is shown in Fig. 19. The results showed that XFEM simulation with mesh size 0.5 mm and 1.0 mm agreed well with the integration methods, but the simulation results with mesh size 1.5 mm and 2 mm lead to an overestimation of fatigue crack propagation life. Hence, it is recommended that the mesh size of fatigue crack propagation is smaller than 1.0 mm for the CT crack propagation simulation using XFEM method.

Load ratio effects
As shown in Fig. 17, fatigue crack propagation simulation of CT specimens of steel grade S355 with an initial crack length 15 mm (Fig. 1a) and S690 with an initial crack length 12 mm (Fig. 1b) is conducted with different R-ratio. The maximum load F is 5.7kN and 3.3kN of CT specimens made of S355 and S690 respectively for the four different loading ratios. A good agreement is observed between XFEM simulation and analytical integration results for both S355 and S690 steel. The fatigue cycle ratio at crack length a = 30 mm for N R=0. 25 Fig. 14) and load ratio R = 0.25 made of S355 and S690 respectively. The crack gradually propagated from 15 mm (start from loading point) to around 35 mm when the fatigue cycle increased to 4.2e5 for CT specimen made of S355 steel, and the crack gradually propagated from 10 mm(start from loading point) to around 20 mm when the fatigue cycle increased to 1.3e5 for CT specimen made  of S690 steel. The crack tip is not strictly propagated uniformly in terms of thickness direction caused by the non-uniform distribution of SIFs. It is noted that the threshold effects are not considered during the fatigue crack propagation simulation due to a lack of experimental data.

Conclusions and future work
In this paper, the Walker equation is employed to fit the fatigue crack propagation rate of steel grades S355 and S690 based on experimental data in the literature to incorporate the load ratio effects. The fatigue life was firstly predicted by the analytical method and was used as a baseline for numerical fatigue crack propagation simulation. A user-defined fatigue crack propagation subroutine was developed based on the Walker equation using phantom nodes-based extended finite element method (PN-XFEM) and Virtual Crack Closure Technique (VCCT). The following conclusions are drawn below: (1) In terms of steel grade S355 with an initial crack a i = 11.0 mm exposed a constant load range ΔP = 5.7 kN, the fatigue cycle difference between crack increment iteration methods and integrated value is within 5% when the crack length increment is smaller than 0.3 mm and is within 10% when the crack length increment is smaller than 0.6 mm. The crack length difference between the fatigue cycle increment method and integrated value is within 5% when the fatigue cycle increment is less than 3000 and is within 10% when the fatigue cycle increment is less than 7000.
(2) The parameters with 95%, 97.7%, and 99% guarantee of Walker equation were obtained by a stochastic analysis using the Monte Carlo analysis. For CT specimen with an initial crack a i = 15 mm exposed a constant load range ΔP = 5.7 kN, the fatigue life calculated by parameters with 95%, 97.7%, and 99% guarantee rate is only 0.03, 0.008 and 0.004 of it predicted from the parameters fitted by experiments respectively of steel grades S355, and is 0.055, 0.030 and 0.016 times of it predicted from the parameters fitted by experiments respectively of steel grades S690. The fatigue crack propagation rate with 95%, 97.7%, and 99% guarantee aims to provide a safe design. (3) The fatigue crack propagation of CT specimens is predicted by three-dimensional XFEM methods with five different boundary conditions. The results showed that XFEM simulation with boundaries applied by two half-cylinders (BC3 in Fig. 14c), applied by reference points tied to half cycles with translation constrain free (BC4 in Fig. 14d), and applied by reference points tied to full cycles with translation constrain free (BC5 in Fig. 14e), agreed well with the integration methods, but the simulation results with boundaries applied by reference points tied to half cycles with all degree freedom fixed (BC1 in Fig. 14a) and applied by reference points tied to full cycles with all degree freedom fixed (BC2 in Fig. 14b) lead to an overestimation of fatigue crack propagation life because the loading point is no longer pin support due to over constrains. In all the boundary conditions, BC3 agreed best with the analytical results. (4) XFEM simulation with mesh size 0.5 mm and 1.0 mm agreed well with the integration methods, but the simulation results with   rate will be investigated in the future study due to a lack of experimental data.

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