Finite element simulation of fracture in a restored premolar tooth using high aspect ratio elements

Introduction: The term cracked tooth syndrome refers to an incomplete fracture of a vital posterior tooth that involves the dentin and occasionally extends into the pulp. There is a very limited number of publications trying to model dentin crack growth using numerical techniques. Therefore, it is essential to numerically model this phenomenon in order to improve the clinical procedures. Methods: A 2D finite element model is proposed to simulate crack initiation and propagation in a restored premolar tooth. The geometric model was based on computed tomography data. A special finite element technique, named mesh fragmentation technique, is used to model and analyze the behavior of the tooth. This technique was used to model cracks in quasi-brittle materials based on the use of interface solid finite elements with high aspects ratio. A tension damage constitutive relation between stresses and strains consistent with the continuous strong discontinuity approach is used to describe crack formation and propagation. Results: The main aspects of modeling technique and procedures are explained in detail as well as the whole results, including both elastic and fracture analyses of the restored tooth. Conclusion: The results of the current fracture analysis show that, under various loading conditions, there is no crack initiation in the restored tooth under typical loading magnitude. However, in the case of tooth with a pre-existing crack, which can be aroused during the restoration process, a crack propagation was observed, while they did not reach a critical fracture state.


Introduction
The term cracked tooth syndrome (CTS) refers to an incomplete fracture of a vital posterior tooth that involves the dentine and occasionally extends into the pulp (Cameron, 1964;Lynch and McConnell, 2002;Rosen, 1982).The symptoms are very variable, making it a notoriously difficult condition to diagnose.The term was first introduced by Cameron (1964), who noted a correlation between restoration size and the occurrence of CTS.Mention is made in the earlier literature of pulpal pain resulting from incomplete tooth fractures.A more recent attempt to define the nature of this condition describes it "[...] a fracture plane of unknown depth and direction passing through tooth structure that may progress to communicate with the pulp and/or periodontal ligament [...]" (Ellis, 2001, p. 428).The condition presents mainly in patients aged between 30 years and 50 years (Ellis et al., 1999;Hiatt, 1973;Snyder, 1976).Two classic patterns of crack formation exist.The first occurs when the crack is centrally located; the second when the crack is more peripherally directed and may result in cuspal fracture.Pressure applied to the crown of a cracked tooth leads to separation of the tooth components along the line of the crack causing pain.
There have been various clinical studies during last decades dedicated to fracture analysis of teeth under various conditions, from restored tooth to root canal treated tooth.Siso et al. (2007) compared the cusp fracture resistance of teeth restored with composite resins.They concluded that, for root filled maxillary premolars, adhesive resin composite restorations increased the fracture resistance of the buccal cusps.Slutzky-Goldberg et al. (2009) investigated the restoration of endodontically treated teeth and concluded that a ferrule of 1-2 mm of tooth tissue coronal to the finish line of the crown significantly Res.Biomed. Eng. 2018 March;34(1): [54][55][56][57][58][59][60][61][62][63][64] 55/64 improves the fracture resistance of the tooth and is more important than a change in core and post material type.Prosthetic rehabilitation of endodontically treated teeth usually requires the use of an endodontic post and coronal core to enhance artificial crown retention.Las Casas et al. (2014) presented a numerical predictive analysis of crack propagation, which may lead to fracture after root reconstruction.A scalar damage model based on the maximum principal stress criterion was used to predict crack propagation.They have concluded that when weakened roots of endodontically treated teeth are treated with adhesive composite reconstruction and post/core restoration, the risk of tensile damage to the root walls is higher with stronger adhesive interfaces.Apparently, localized failures of the interface corresponding to peak stress areas decrease the risk of damage to the root dentin walls.Munari et al. (2015) studied and compared the areas of stress concentration in a three-dimensional premolar tooth model with anisotropic or isotropic enamel using the finite element method.Because tooth structures are more resistant to compression, damage such as the formation of cracks and fracture of tooth tissues are likely to be caused by tensile stress from the eccentric contacts of unbalanced occlusion.For more information in this topic, see for example (Banerji et al., 2010;Lin et al., 2013;Lubisich et al., 2010;Silva et al., 2012;Yahyazadehfar et al., 2013Yahyazadehfar et al., , 2014)).
There is a limited number of publications where the researchers tried to model crack growth in teeth using the finite element method.The modeling technique can be chosen from various computational methods, which are available in the engineering applications, such as the finite element method (Cornacchia et al., 2010) or extended finite element method (Zhang et al., 2013).Both methods require techniques to track the crack path during the analysis to provide information about the position of the crack surfaces for the corresponding discontinuous kinematic enrichment.These techniques are relatively simple to represent a few cracks in 2D analyses but can be very complex and even unsuitable for problems involving multiple crack surfaces in 3D analysis (Jager et al., 2008;Manzoli et al., 2014Manzoli et al., , 2016)).
To avoid the necessity of crack tracking schemes, a technique based on the insertion of special interface elements in between regular elements of the mesh was proposed by Manzoli et al. (2016) called mesh fragmentation technique.This technique was developed for modelling cracks in quasi-brittle materials based on the use of interface solid finite elements.As explained in Manzoli et al. (2012), in the limit case when the thickness of the interface elements tends to zero (and the aspect ratio tends to infinite), these elements present the same kinematics as the Continuum Strong Discontinuity Approach -CSDA (Oliver and Huespe, 2004;Oliver et al., 1999).Consequently, the interface elements are able to describe the kinematics associated to discontinuities, so that the crack formation can develop along the boundaries of the regular elements.Therefore, the analyses can be performed integrally in the context of the continuum mechanics, which is an advantage of the method applied.
Current work presents a modeling approach for the fracture analysis in tooth structure using a finite element analysis and special crack propagation techniques, like the method proposed by Manzoli et al. (2016).

The mesh fragmentation technique
The mesh fragmentation technique is based on the use of interface solid finite elements with a high aspect ratio (Manzoli et al., 2016), which are inserted in between standard finite elements (regular mesh).The crack formation and propagation is handled via these interface elements using an appropriate continuum damage model.Figure 1a presents the main steps of the algorithm implemented (in Matlab©) for the proposed mesh fragmentation technique for 2D problems.This procedure has three steps: generation of the standard FE mesh to be fragmented (Step 1), creation of new nodes in order to isolate the regular elements moving these nodes toward the center of gravity of each corresponding regular element, generating gaps between them (Step 2).Finally, insertion of pairs of interface elements with a high aspect ratio in between adjacent regular finite elements (Step 3).More details about this procedure can be found in Manzoli et al. (2016).
Standard three-node triangular finite elements are considered here to describe the features of the interface solid finite elements in 2D modeling, as illustrated in Figure 1b.The geometry of these elements can be characterized by the position of their nodes according to a local Cartesian coordinate system ( ) , n s , defining the unit vector, n, normal to the element base, and the height, h, given by the distance between the node 1 and its projection on the element base 1'.Following the standard finite element approximations, the strain field of the solid finite element can be expressed by: where B is the strain-displacement matrix and D is the nodal displacement vector of the element.Then, the phenomenon of crack initiation and propagation through the interface solid finite elements is described by a tension damage model.This model is defined by the following constitutive relation: where σ is the nominal stress; is the scalar damage variable;  is the fourth order elastic tensor; ε is the strain tensor; and the product : ε  defines the effective stress tensor σ.When the scalar damage variable, d, reaches its maximum value, i.e., 1, the interface element fails completely.The value of this variable is obtained by using the fracture toughness and tensile strength of the material.Further information on this technique and also the interface solid element can be found in (Manzoli et al., 2012(Manzoli et al., , 2016)).
The mesh fragmentation approach is completed by a continuum tension damage model formulated to describe the formation and propagation of cracks along the interfaces with high aspect ratio (ratio between the largest to the smallest dimension).Therefore, as this relation tends to zero, the aspect ratio increases, as well as the strains, approaching to the kinematics of the strong discontinuity, which is consistent with the CSDA approach (Oliver et al., 1999), as showed by Manzoli et al. (2016).
For a monotonic increase (in mode-I) of the normal displacement, the evolution of the normal stress becomes: where nn σ is the component of the stress normal to the base of the element ( ) . .= ,where t f is the tensile strength of the material and  is the softening parameter, the mode-I fracture energy, f G , i.e., the energy dissipated in a complete degradation of the interface element in mode-I, is given as: When h tends to zero, the fracture energy tends to a non-null value given as: from which one can define the softening parameter in terms of the material properties as: The values of softening parameter and tensile strength, along with material properties are used to define the interface solid finite elements inside of the mesh fragmentation procedure.As it is well known, the constitutive tangent operator may become singular in the strain softening regime, and, as a consequence, the solution of the resulting systems of nonlinear equations using a fully implicit scheme may not be achieved.According to Oliver et al. (2006) this problem may also be present in the numerical modeling of material failure even when a powerful continuation method to pass structural points is used (e.g.arc length methods to transverse limit and turning points).To address this problem, the implicit-explicit (IMPL-EX) integration scheme proposed by Oliver et al. (2006Oliver et al. ( , 2008) ) is used for the integration of the stress-strain relation of the damage model.Further explanations on the mesh fragmentation technique and the constitutive damage model enriched with IMPL-EX scheme can be found in (Manzoli et al., 2016;Rodrigues et al., 2016).

Modelling procedure
A geometric 3D finite element model of a maxillary premolar tooth is built.The model also includes a representation of the mandible bone, which is based on tomography images.Six different materials are considered: enamel, dentin, periodontal ligament, trabecular bone, cortical bone and composite resin (as the restoration material).Enamel and dentin are the basic tissues that constitute a sound human tooth.Resin and ceramic are the restorative materials.
The periodontal ligament, which makes the link of the tooth with the mandible bone, is also considered in the modeling process.Figure 2a illustrates the distribution of these materials in the domain.
The shape of the restoration is in accordance with the practical samples which is similar to the one used in the work of Hamouda and Shehata (2011), Figure 2b and 2c.The overall loading considered is decomposed of multiple concentrated loads.The lateral external surfaces of the mandible bone are clamped as the tooth boundary conditions.There are various investigation reporting the loading magnitude for a tooth under different loading conditions: a maximum load of 522 N was proposed in (Lyons et al., 1996;Proeschel and Morneburg, 2002) using the experimental techniques; a FE study on a mandibular premolar tooth with a loading of 400 N was conducted in (Palamara et al., 2002); a maximum biting force of 453 N was reported in (Litonjua et al., 2004) based on experimental evaluation technique; a load of 380 N was applied on an implanted premolar tooth (Hisam et al., 2015); and, a maximum load of 600 N was considered in (Misch, 2015) applied on premolar tooth.Therefore, a maximum loading magnitude of 600 N is considered for this study, to be applied on the tooth structures.The mechanical properties of each material and tooth tissue are given in Table 1.The source chosen as reference directly affects the obtained results.
A 2D model is created from a section plane along the middle part of the 3D model.Figure 3 shows the geometry, boundary conditions, loading, and mesh schematic for the 2D analysis.There are 13,903 triangular elements (linear 3-node triangular elements) in the discretized model, with 5,419 elements for dentin; 774 elements for restoration; 1,017 elements for enamel; 1,261 elements for periodontal ligament; 3,778 elements for cortical bone; and 1,654 elements for trabecular elements.As it was explained in the previous section and considering the literature, the loading magnitude is set up to a maximum of 600 N. Figure 3c shows the location of some initial cracks to be used for fracture analysis.

Elastic analysis
The stress distributions for this model, from an elastic analysis, are shown in Figure 4.The loading orientation has an important impact on both displacement and stress distributions.As it can be seen from Figures 4c and 4d as well as Figures 4g and 4h, elastic analysis from mesh fragmentation technique has good agreement with the elastic analysis obtained from ABAQUS software.
In the first loading condition (Figure 4a which is lingual load), the applied load is predominant in the vertical direction, as can be seen in Figure 3a (for the load L1).Thus, well-distributed compressive stress is expected (S1 is very small).For the second load case (Figure 4b) for the load L2 from Figure 3a, which is buccal load), there is a significant component of load in the horizontal direction, in which it produces a bending effect on the tooth, and hence localized tension stress in the left side and compressive stress in the right side.

Fracture analysis
After the elastic analyses, interface elements are inserted in between the elements of dentin to model the crack propagation in the dentin.These interface elements were inserted also between dentin and restoration boundaries in order to model the bounding failure between these two parts during the crack propagation process.In addition, two different cases are modeled: dentin with and without a pre-existing 0.3 mm crack at different positions: lingual, buccal, and central as shown in Figure 3c.Similarly to the elastic analyses, three different loading cases are also considered here to see the behavior of crack propagation in the dentin.
Besides elastic modulus and Poisson's ratio, there are other parameters that must be defined for those parts that contain interface elements, such as tensile strength.The considered tensile strength (to be used for damage evaluation in interface elements) is listed in Table 1.Since the interface elements were only inserted in the dentin region, let's consider an average value for tensile strength and fracture energy from Table 1, with a value of 70 MPa for tensile strength and 648 J/m 2 for the fracture energy.Then, first start to model the crack initiation and propagation in the dentin with this parameter, with different loading conditions.
According to Figure 5, for the tensile strength value of 70 MPa, which is the average values obtained from Table 1, there is no crack propagation for both lingual and buccal loading cases until loading value of 600N, while for buccal case there is small crack propagation under loading value of 750N.For a tensile strength value of 50 MPa, there is a small crack propagation for lingual loading of 600N, while for buccal loading of 600N there is no crack propagation.Finally, for the tensile strength value of 32 MPa, we can see a quite large crack propagation for lingual loading equal to 600N, but a small amount of crack propagation ocuured for buccal loading with the same loading value.Therefore, a value of 32 MPa will be considered as the tensile strength value of the dentin in this research, for the failure point of the interface elements in the dentin part.It is the smallest value from Table 1, and used to study the fracture path under various loading conditions, while the results for the larger values are also shown in the previous figures.
The fracture path from different loading cases along with various pre-existing crack positions, as well as model without any crack are shown in Figure 6 and Figure 7 (with tensile strength of 32-80 MPa, with an average fracture energy of 648 J/m 2 ).The models without cracks are analyzed in order to verify whether any crack would be initiated under different loading orientations and various loading magnitudes.Then, for each loading cases, a corresponding crack position is inserted in the model.As an example, for lingual loading the lingual crack is placed in the model.For lingual+buccal loading case, all the three pre-existing cracks are placed in the model and their propagation behavior is studied separately.

Discussion
This work presented an elastic as well as fracture analyses of a premolar tooth under various loading conditions.The fracture modeling was done using the mesh fragmentation technique.The main aspects of this technique were described in detail and the whole steps to conduct this study were brought and discussed.
As shown in Figure 4, the stress concentrations for the buccal loading are a little far from the bottom side of the restoration, closer to the pulp boundaries, while the stress concentrations are close to the restoration bottom edges for lingual and lingual+buccal loading cases.In the case of lingual+buccal loading, the stresses are mainly concentrated near to the loading points, which is expected.Other than these points, the stress distributions are almost the same for the lower bottom corners of the restoration area.
In addition, as it can be seen from Figure 6 and Figure 7, the lingual and lingual+buccal loading cases have resulted in bigger crack propagation paths, mainly due to their loading configuration as well orientation regarding the applied load surface, i.e., the restoration surface.The realistic loading value for a tooth is chosen up to a maximum of 600 N (Misch, 2015).Therefore, all the crack propagation until the loading of 600 N is acceptable and can be considered closer to the realistic situation.However, the larger load magnitudes were shown as well to demonstrate the crack propagation behavior under other less common loading conditions.
Based on authors' investigations, there is no numerical or experimental study similar to the current work, therefore no comparison with any previous results was brought in this work.Finally, here are the main conclusions from this research: -The results indicate that the tooth restored without any discontinuity or cracks did not reach critical stresses under usual mastication loads; -The pre-existence of small initial cracks, due to the restoration process, can nevertheless lead to critical crack propagation; -The loading type also plays a very important role in the crack propagation driving force.The lingual load leads to some crack propagation while the buccal load leads to very small crack propagation.
In addition, the combination of these two loads leads to a meaningful crack propagation path; -Another important issue is the loading magnitude, a parameter for which a large variation is possible; -The material properties also play a very critical role here, as the interface elements behavior mainly depends on the material properties of the region of interest.Therefore, they must be selected/defined in such a way that provides useful and reliable results.
Assuming an exponential softening law of the form:

Figure 1 .
Figure 1.Steps of the mesh fragmentation technique and the definition of the interface element: (a) 2D mesh fragmentation process, from the generation of the standard FE mesh (step 1) to the insertion of interface elements (step 3) (Rodrigues et al., 2016); (b) Three-node triangular element, a 2D interface solid finite element.

Figure 2 .
Figure 2. Representation of the 3D geometric model and the schematic of the restoration area: (a) 3D geometric model of the tooth; (b) Real geometry (Hamouda and Shehata, 2011); (c) Model geometry used in this study.

Figure 3 .
Figure 3. Schematic of the boundary conditions, loading, meshing model, and locations of the initial cracks.(a) Geometry and loading; (b) Mesh discretization; (c) Locations of the initial cracks.

Figure 5 . 64 Figure 6 .
Figure 5. Fracture results for lingual and buccal cases with different tensile strength values.(a) Lingual loading case; (b) Buccal loading case.

Figure 7 .
Figure 7. Crack propagation path for lingual+buccal loading case.(a) Without any pre-existing crack; (b) With a pre-existing crack in both sides; (c) With a pre-existing crack located at the bottom center of the cavity.

Table 1 .
Mechanical properties of each tooth parts.