In recent years,the health and esthetics of teeth become increasingly important and attract more and more concern from people. However,tooth fracture is always a threat to the health of human teeth. Therefore,it is of great need to understand the crack propagation and fracture behavior of human teeth comprehensively so as to help to prevent tooth disease and maintain tooth health. Human teeth are composed of enamel,dentin,and pulp. Dentin is a highly mineralized tissue and occupies most of the volumes of teeth. It locates between enamel and pulp with the function of protecting pulp. Therefore,it is imperative to establish a good understanding on the fracture properties of dentin.

Conventional strength tests cannot provide enough insight into the fracture mechanism of dentin^[1] due to the complexity of fracture in dentin. The fracture of dentin is composed of crack initiation and crack propagation which cannot be properly characterized by strength. Recently,fracture mechanics has been introduced to study the fracture behavior of dentin^{[2, 3, 4, 5, 6, 7, 8, 9]}. With the help of fracture mechanics,it is found that the fracture toughness of dentin depends on the orientation^{[2, 3]}. However,these studies are based on the assumption that the fracture toughness of dentin is constant throughout the whole dentin. Kruzic et al.^[6] found that the fracture toughness of dentin increases with crack propagation,i.e.,dentin possesses a rising resistance curve (R-curve). Multiple toughening mechanisms had been proposed to account for the rising R-curve property,e.g.,crack bridging by uncracked ligaments at the crack wake and microcracking^[6]. It is also found that hydration and aging have great influence on the fracture toughness^{[6, 7, 8]}. Although many work has been conducted on the fracture mechanism of dentin from different aspects,the critical fact is ignored that dentin is a typical heterogeneous material. The mechanical behavior is not only related to its orientation but also the distance from the dentin-enamel junction (DEJ). However,the region dependent fracture behavior has not been reported in the literature yet.

Dentin has a hierarchical structure and is composed of mineral platelet and protein matrix at nano-scale[10−11]. At micro-scale,dentin consists of numerous tubules. The deformation and fracture behaviors of dentin closely depend on the structure and the corresponding properties at different scales. Recently,an experimental study indicated that the fracture toughness of human dentin varied with the distance to the dentin-enamel junction^[12] . Therefore,it is necessary to employ the multiscale numerical simulation to study the mechanism of fracture of dentin. To our knowledge,the fracture behavior of human dentin has also not been simulated by using the multiscale numerical model so far.

The virtual internal bond (VIB) model,which is developed by Gao and Klein^[13],is a typical multiscale method and had been successfully applied to various fracture simulations^{[14, 15, 16, 17, 18]}. In the VIB model,the macro-mechanical properties are derived from the properties of bonds between the material particles. The advantage of the VIB model is that it does not require an extra fracture criterion. However,the Poisson’ ratio is fixed in the VIB model,restricting the application of the VIB model. Zhang and Ge^[19] developed the virtual multidimensional internal bond (VMIB) model to solve this limitation. The VMIB model inherits the advantages of the VIB model,while at the same time,introduces a variable Poisson’s ratio. It has been successfully employed to some crack growth simulations^{[20, 21, 22]}.

In order to take the effect of the hierarchy and heterogeneity of human dentin into account, the present work simulates crack propagation in human dentin by using the VMIB model and Monte Carlo method to study the region effect on the R-curve to provide understanding to the fracture behavior of human dentin. The primary objective is to study the fracture behaviors of dentin at different distances to the DEJ. 2 Methods

In the VMIB model,the material is regarded as the assembly of numerous material particles^{[19, 20]}. In the present study,dentin is modeled as a material consisting of numerous dentin particles. Virtual bonds exist between adjacent dentin particles. The macroscopic constitutive relationship can be derived from the stiffness of the microscopic bond. A detailed description of the VMIB model is provided here. The elastic tensor in the VMIB model can be calculated as follows^[20]:

Click Browse formula

where k is the normal stiffness,and r is the shear stiffness; ξ is the unit vector,and in the sphere coordinate system,ξ = (sin ϑ cos ϕ,sin ϑ sin ϕ,cos ϑ)^T; η is the vector perpendicular to ξ,η′ = ξ × (x₁ × ξ),η′′ = ξ × (x₂ × ξ),η′′′ = ξ × (x₃ × ξ),and x_i (i = 1,2,3) is the base vector of the coordinate; D(ϑ,ϕ) is the bond density distribution function.

The bond stiffness function has great influence on the macroscopic properties of material. According to the previous studies^{[13, 14, 17]},a two-parameter cohesive law,U′(l) = Al exp(−l/B), can be used to model the fracture behavior. In the equation,U′(l) is the force of bond. A is the initial modulus of material,and B is the strain which corresponds to the ultimate tensile strength. Zhang and Ge^[19] also discussed the application of the two parameter cohesive law. Considering that the two-parameter cohesive law can describe the stress-strain response of materials^[19],the uniaxial tensile response of human dentin is written as

Click Browse formula

In the present study,the plane strain condition is assumed and the following phenomenological stiffness functions are adopted:

Click Browse formula

where ε_t is the strain corresponding to the ultimate tensile strength; l is the stretch of bond given by l = |ξ_iε_ijξ_j|,and ε_ij is the strain tensor; α₁,α₂ are the rotation angles,and α₁ = |ξ_kε_klη_l′|, α₂ = |ξ_kε_klη_l′′|; µ is the Poission’s ratio; k₀ = ,and r₀ = ^[20]. C₁,C₂, and C₃ are model parameters,which can be determined by fitting the uniaxial tensile response of materials. In the present study,C₁ = 0.724 81,C₂ = 0.709 53,and C₃ = 1.465 20. The first stiffness function in (3) is used in ^[19]. Under plane strain condition,the elastic tensor can be obtained by substituting (3) into (1) as follows:

Click Browse formula

Human dentin can be regarded as an isotropic material^[23] and D(ϑ) = 1 for this case^[19]. The VMIB model can be used to model the fracture behaviors of different materials by assigning the corresponding parameters of the materials. In this case,the strain corresponding to the ultimate tensile strength of human dentin ε_t is introduced to delineate the mechanical properties of human dentin.

The bond stiffness characterizes the interactions between material particles. According to (3),the bond stiffness decreases with the increasing deformation. At the micro-scale,the decrease of bond stiffness indicates that the distance between particles is increased. To an extreme extent,when the distance is large enough,the macroscopic cracking occurs^[19]. According to the study of Staninec et al.^[24],the ultimate tensile strength of dentin can be characterized by using the following Weibull distribution^[24]:

Click Browse formula

where P is the fracture probability of dentin; σ is the ultimate tensile strength of human dentin; σ₀ is the scale parameter,and m is the Weibull modulus. In the present study,the strain corresponding to the ultimate tensile strength ε_t is considered to be a random variable which obeys the Weibull distribution.

In this study,three regions with different distances to DEJ,i.e.,the superficial dentin,the middle dentin,and the deep dentin,are considered. Considering the symmetry of the CT specimen shown in Fig. 1(a),only half of the CT specimen is simulated. The corresponding Weibull parameters are listed in Table 1^[24,25]. The simulation region is meshed with threenode triangular isoparametric elements. As shown in Fig. 1(b),the left-side region of the model is meshed with ordinary linear elastic elements while the right-side of the model is meshed with VMIB elements. Totally,11 342 elements are meshed in the model. A displacement-controlled load is applied at the hole to achieve the quasi-static loading. The effective crack length is calculated by using the compliance method^[26] . Accoding to the ASTM Standard^[27],the Rcurve of the material is determined^[28]. The R-curve of the middle dentin is firstly calculated to verify the employed method and determine the number of the random samples. Then,this numerical method is applied to predict the R-curves of the superficial dentin and the deep dentin,and the corresponding R-curves are also plotted to reveal the region dependent fracture behavior.

Fig. 1 Compact tension specimen geometry (a) and finite element model of CT specimen (b)

Figure options

Table 1. Weibull distribution parameters for three dentin locations^{[23, 24]}

Table options

Two hundred random samples are calculated in the middle dentin by using Monte Carlo simulations. The mean load-displacement curve is shown in Fig. 2(a),and the corresponding R-curve,which illustrates the fracture toughness with respect to the crack extension,is shown in Fig. 2(b).

Fig. 2 Mean load-displacement curve (a) and R-curve for crack propagation (b) in middle dentin

Figure options

Nazari et al.^[8] and Zhang et al.^[9] have studied the aging effect on the fracture toughness of human dentin by using an experimental approach,showing that the initial toughness,i.e., the stress intensity factor at the onset of the crack extension,is in the range of 1 MPa·m- 1.4 ,and the plateau toughness,i.e.,the maximum fracture toughness,is in the range of 1.08 MPa·m-1.75 . The initial and plateau toughnesses of the middle age group (356age655) are (1.22±0.06) and (1.43±0.1) ,respectively. The present numerical simulations indicate that the initial toughness of the middle dentin is 1.155 , and the plateau toughnesses of the middle dentin is 1.46 ,consistent with the experimental results for middle-aged patients.

The results of all Monte Carlo simulations are shown in Fig. 3. The crack extension is expressed in the dimensionless form ∆a/W,where W is the specimen width. Figure 3(a) shows that the load-carrying capacity of dentin varies in the different regions. The mean values of the peak loads for crack growth in the superficial,the meddle,and the deep dentin are 12.83 N, 10.11 N,and 7.12 N,respectively.

Fig. 3 Load-displacement curves for all Monte Carlo simulations where black curves are mean loaddisplacement curves for three dentin locations (a) and R-curves for crack propagation in three dentin locations (b)

Figure options

Among the three dentin regions,the superficial dentin has the greatest load-carrying capacity while the deep dentin possesses the smallest one. Figure 3(b) shows that all the superficial,the middle,and the deep dentin exhibit rising R-curves. The initial toughness of the superficial,the meddle,and the deep dentin are 1.4 ,1.155 ,and 0.79 ,respectively. The plateau toughness are 1.86 MPa·m 12 ,1.46 ,and 1.03 ,respectively. 4 Discussion

The numerical simulations show that the fracture toughness of human dentin increases with the crack growth,i.e.,human dentin possesses a rising R-curve. The rising R-curve behavior of dentin is attributed to various toughening mechanism. Kruzic et al.^[6] found that the crack bridging by uncracked ligaments in the crack wake was responsible for the rising Rcurve. Microcracking was also reported as a toughening mechanism^[2]. In the VMIB model,the interaction between material particles is modeled by the cohesive law. Since the cohesive law represents the effects of the crack bridging and fracture processing zone near the crack tip^[29], the rising R-curve behavior predicted by the present model is also capable of describing the contributions of crack bridging and microcracking. For the consistency between the numerical simulation and the experimental results,the fracture behaviors of the deep and superficial dentin are also analyzed by changing the parameters listed in Table 1.

The obtained results clearly indicate that the fracture toughness of human dentin closely depends on its location. Superficial dentin has the strongest resistance to fracture,while deep dentin is the region most prone to crack. The previous study conducted by Bajaj and Arola^[30] identified that the fracture toughness of human enamel was higher than that of the middle dentin. Imbeni et al.^[31] observed a crack propagated from enamel to dentin,and found that the crack traversed the DEJ and stopped in the superficial dentin. It is reasonable to conclude that there should be a region which possesses very high fracture toughness and locates very closely to the DEJ. In such condition,the propagating crack can be arrested near the DEJ. From the numerical analysis,the maximum fracture toughness of superficial dentin reaches 1.86 ,providing reliable resistance for the crack propagation from the outer enamel.

The simulations provide more comprehensive understanding of fracture toughness of human dentin. Due to the difficulties in preparing specimens,previous experimental studies were focused on the fracture toughness of middle dentin. The numerical study provides a way to characterize the dependence of the fracture behavior of dentin on its location,i.e.,the superficial, the middle,and the deep dentins,developing an understanding of the fracture behavior of dentin as a non-homogeneous material. The findings in this study also provide guidelines for the tooth restoration. It is found that tooth fracture is the major threat to the restored teeth. Based on the present study,the materials with graded fracture toughness are required to reduce the propensity of fracture of the restored teeth. 5 Conclusions

Numerical simulations of fracture of human dentin are carried out to study the fracture toughness of dentin. The simulations show that human dentin exhibits rising R-curve behaviors. The R-curve property of human dentin depends on its location. Superficial dentin has the largest resistance to fracture,while deep dentin possesses the weakest fracture toughness. The present study also demonstrates that the multiscale numerical simulations based on VMIB theory and Monte Carlo method could be used to characterize the fracture behaviors of human dentin as well as other biological materials.