7717


FACTA UNIVERSITATIS  

Series: Mechanical Engineering  

https://doi.org/10.22190/FUME210619072C 

© 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND 

Original scientific paper 

FRACTURE PARAMETERS EVALUATION FOR THE CRACKED 

NONHOMOGENEOUS ENAMEL BASED ON THE FINITE 

ELEMENT METHOD AND VIRTUAL CRACK CLOSURE 

TECHNIQUE 

Yan Cai1,2, Liming Zhou1, Ming Li1 

1School of Mechanical and Aerospace Engineering, Jilin University, China 
2Hospital of Stomatology, Jilin University, China 

Abstract. To accurately solve the fracture parameters of enamel, we have established 

computational nonhomogeneous enamel models and constructed the fracture element of 

enamel dumb nodes, based on the enamel mineral concentration, nonhomogeneous 

mechanical properties, and virtual crack closure technique. Through the commercial 

finite element software ABAQUS and the fracture element of the enamel dumb nodes, we 

have established the user subroutines UMAT and UEL, which enabled solving of the 

energy release rates of the nonhomogeneous enamel structure with cracks. The stress 

intensity factors of central cracks, three-point bend and compact stretched enamels, and 

double-edge notched stretched enamels are determined. By comparing them with 

analytical solutions, we have proved that the fracture element of the enamel dumb nodes 

is highly accurate, simple, and convenient. In addition, the cracks can be other elements 

rather than singular or special elements; they show versatility and other advantages. 

The stress intensity factor of the dental enamel can be solved more realistically. Thus, a 

new numerical method for prevention and treatment of dental diseases is provided. 

Key Words: Nonhomogeneous Enamel Models, Crack, Virtual Crack Closure Technique, 

Finite Element Method, Fracture Element of Enamel Dumb Nodes 

1. INTRODUCTION 

Enamel is the most important component of the outermost tissue of human teeth [1-3], 

which serves as a barrier for protecting the vital pulp and the dentin. On a microstructural 

level, a sequence of crystal rods (or “prisms”) with diameters of 4–6 μm compose the 

enamel tissue. The crystal rods stretch from the dentin-enamel-junction (DEJ) to the 

                                                           
Received June 19, 2021 / Accepted November 28, 2021 

Corresponding author: Liming Zhou, Ming Li 

School of Mechanical and Aerospace Engineering, Jilin University, Changchun, Jilin 130025, China 

E-mail: lmzhou@jlu.edu.cn (Liming Zhou), liming940411@gmail.com (Ming Li). 

mailto:lmzhou@jlu.edu.cn
mailto:liming940411@gmail.com


2 Y. CAI, L. ZHOU, M. LI 

occlusal surface and are arranged in a parallel manner [4]. The interface between the rods is 

a protein-rich district considered as the rod sheath. Enamel is a composite material that has 

a unique and complex hierarchical structure. The unique properties of enamel distinguish it 

from other biological tissues like bone or dentin [5]. In the body, enamel is the most highly 

mineralized tissue, being composed of approximately 96% minerals (primarily apatite 

crystals containing Ca2+ and P3-, small amounts of other phosphate crystals, and small 

amounts of other minerals), and only 4% organic matrix and water (3% water and 1% 

proteins) by weight [6]. Teeth are well-designed to endure repeated chewing loads during a 

lifetime. To support oral function, the enamel should be equipped with adequate hardness, 

stiffness, and fracture resistance. Although the toughness of the protective enamel coat is 

comparable to that of glass, this material is highly brittle [7]. 

In dental clinic, incomplete tooth or tooth fracture is a general pathological feature [8, 9]. 

Considering the physiological characteristic of enamel and the circle force contact incurred 

during mastication, fatigue can lead to failure. Fractures that run as ribbon cracks around the 

enamel shell walls may occur as a result of macro-contacts in the cusp regions with hard 

millimeter-scale food objects (e.g., nuts, seeds, and hard particulates) [10-12]. These cracks 

are distributed along a normal direction to the tooth surface orientation and remain entirely 

contained within the enamel. Some cracks develop rapidly in clinic, and in severe cases, 

they can cause deterioration or even ultimate loss of tooth function. However, some 

fractured teeth can be used during life. Therefore, considering how to effectively protect the 

tooth activity and reduce or even block the propagation of cracks, further analysis of the 

mechanical characteristics of the dental enamel is required. 

Over the last decade, enamel has been modeled as a homogenous material [13, 14], but 

recently, it has been explored with variations in its mechanical characteristics. Enamel 

hardness and elastic modulus decrease from the occlusal surface of enamel to the DEJ [15, 

16]. The variations are probably due to chemistry, microstructure, and prism alignment. 

However, Cuy et al. [17] suggested that enamel is significantly correlated with changes in 

average chemistry: the concentrations of hydroxyapatites P2O5 and CaO are larger over the 

chewing surface and gradually decrease to the softer DEJ, but the concentrations of Na2O 

and MgO change in the opposite direction. Braly [18] and Park [19] also reported similar 

variations of mechanical properties on the axial section of the enamel [20]. More recently, 

nanoindentation tests showed that the enamel elastic modulus and hardness decreased from 

≥6 GPa and ≥120 GPa in the occlusal surface to less than 3 GPa and 70 GPa in the DEJ, 

respectively [21-23]. As reported, the crack growth resistance of enamel increases (rising 

R-curve) with crack extension from the outer to inner part; the law of variation of toughness 

is associated with the distance to DEJ [24]. 

Due to the small size of cracks, it is difficult to build an in-vitro experimental model for 

early cracked teeth. Thus, the finite element method (FEM) becomes relatively practical for 

studies on cracked teeth [25-27]. Recently, FEM has been widely used to analyze the causes 

of cracked teeth, but the existing studies have focused on the stress distribution, 

biomechanical factors, crack expanding behavior, and fatigue life of cracked teeth [28-30]. 

Moreover, due to methodological limitations, previous studies mostly regarded enamel as 

an isotropic material for mechanical analyses. 

Using the analytical method, the fracture parameters of a simple geometric structure or a 

special load form of homogenous materials can be derived [31-33]. However, this method 

cannot be applied to complex systems or boundary conditions, as the calculation process 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 3 

becomes very expensive and time-consuming. Nevertheless, these limitations can be 

overcome by using other numerical methods. 

Nowadays, FEM is extensively employed in fracture analysis of nonhomogeneous 

materials. The elastic modulus coefficient of nonhomogeneous materials is a smoothing 

function of specific coordinates, also associated with the establishment of element stiffness 

matrix. For two-dimensional (2D) problems, the stress intensity factors for modes (modes I 

and II) can be comparatively evaluated through three nonhomogeneous material-tailored 

approaches: path-independent J-integral, virtual crack closure technique (VCCT) [34-37], 

and displacement correlation. 

VCCT was used for 2D cracked specimens and thereafter expanded to 

three-dimensional crack problems [38-40]. Then, a method combining FEM and VCCT for 

the study of kinking cracks in a homogeneous material was introduced [41-43]. The energy 

release rates for cracks were computed by VCCT [44], and the interface cracks were studied 

with the virtual crack extension method and VCCT [45-47]. Krueger [48] presented an 

overview of the VCCT historical developments and discussed its different applications. 

Leski [49] provided the general conditions for implementation of VCCT in 

MSC.PATRAN. 

Azimi et al. [50] probed into the special cracked lap shear specimen deboning using 

VCCT. Zeng et al. [51] extended VCCT to investigate the fracture mechanics parameters 

and simulate crack propagation within the framework of cell-based smoothed FEM. 

Banks-Sills and Farkash [52] explored VCCT to calculate stress intensity factors of two 

modes (modes I and II) for an interface crack between two linear elastic, isotropic, and 

homogeneous materials. VCCT was also extended to the investigation of the element-free 

Galerkin method [53, 54]. Based on graded FEM, a 2D-VCCT interface element was 

proposed for the dynamic fracture analysis issues [55]. 

Despite the secondary development of the commercial finite element software (ANSYS 

or ABAQUS), the crack problem analysis is mainly focused on homogeneous materials, and 

rarely on nonhomogeneous materials [56-58]. 

In this study, we built nonhomogeneous enamel models based on the commercial finite 

element analysis (FEA) software ABAQUS and VCCT. A dummy node enamel fracture 

element was established, and the energy release rates of nonhomogeneous cracked enamel 

were solved. The method provided better results compared to the analytical method in 

addressing nonhomogeneous mechanical properties. 

2. MATERIALS AND METHODS 

2.1. Nonhomogeneous enamel model 

Previous studies mostly considered the enamel as an isotropic material for mechanical 

analysis. However, the tooth mineral concentrations gradually change from tooth surfaces 

to the DEJ [17] (Fig. 1). Based on Ref. [17], enamel is a nonhomogeneous material (from 

tooth surfaces to DEJ, the percentages of P2O5, CaO, Na2O, MgO and other components 

change in certain functions). For the model proposed, two materials were considered: 

Material C, with a composition of 94% P2O5 + CaO, 6% Na2O + MgO and other 

components, and material M, with a composition of 89% P2O5 + CaO, 11% Na2O + MgO 



4 Y. CAI, L. ZHOU, M. LI 

and other components. Fig. 2 shows the gradient model of nonhomogeneous enamel 

materials. The enamels analyzed were composed of materials M and C, whose volume 

fractions continuously changed along certain dimensions of the model (Fig. 2). The 

Mori-Tanaka models and the rule of mixture are common methods to evaluate the effective 

elastic properties of grade composites. The effective moduli of the two materials were 

homogenized through these two methods. The effective property associated with volume 

fraction exponent follows a power law: 

   m m c cP z P V PV  , 1m cV V   (1) 

  0.5
n

c
V z h  ,  2 2,  0h z h n       (2) 

where subscripts m and c in Eq. (1) denote materials M and C, respectively; P denotes the 

Poisson’s ratio or Young’s modulus; z denotes the thickness coordinate; n is the power 

exponent factor. Eq. (1) presents the volume fraction variation along non-dimensional 

thickness (z/h). Volume fraction Vc variation with different power exponent factors n is shown 

in Fig. 3. 

 

Fig. 1 Map of Young’s modulus throughout the enamel by nanoindentation 

Material C

Material M

75% Material C + 25% Material M

50% Material C + 50% Material M

25% Material C + 75% Material M

 

Fig. 2 Gradient model of nonhomogeneous enamel materials 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 5 

-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5

0.0

0.2

0.4

0.6

0.8

1.0

n=15.0

n=5.0

n=2.0

n=1.0

n=0.5

n=0.2

 

 

V
c

z/h

n=0.0

 

Fig. 3 Variation of volume fraction Vc function along the non-dimensional thickness with 

varying n 

2.2 Finite element method 

The elements were established by directly adopting the true parameter at the Gauss 

integration points of each element [59]. The equilibrium equations are given as: 

 e e eK U F  (3) 

where Ue represents the nodal displacement vector and Fe represents the load vector. The 

stiffness matrix of FEM can be expresses as: 

 
e

e e T e e e
( ) ( ) d


 K B D x B  (4) 

where e represents the problem domain of finite element e, Be and De(x) represent the 

strain-displacement matrix and constitutive matrix variable, respectively. The 

strain-displacement is related to the gradients of the interpolating functions. In this study, the 

elasticity matrix De(x) = De(x, y) was assumed to be related with spatial coordinates and 

specified at each Gaussian integration point. 

Elastic modulus E and Poisson’s ratio ν in elastic matrix De(x) can be established by 

interpolation for the nonhomogeneous materials [60]: 

 

1

ε

i i

i

E N E


   (5) 

 

1

ε

i i

i

v N v


   (6) 

where Ni denotes the FEM shape function corresponding to node i. Ei and vi denote the 

material properties at node i of the element. ε is the node number of the element. By using Eqs. 

(5) and (6), the actual variations of the material properties in an element can be approximated 

by polynomial forms.  



6 Y. CAI, L. ZHOU, M. LI 

2.3 Fracture element of enamel dumb nodes 

The enamel cracks commonly seen in the clinic are closely related to the acute bite force, 

temperature change, and congenital developmental defects of the enamel. The properties of 

the mechanical field near the crack tip in nonhomogeneous enamel are the same as in 

homogeneous materials [61]. However, the stress intensity factors and fracture toughness 

parameters follow gradient functions. For VCCT, the strain energy was assumed as the 

same when releasing from a certain degree of crack extension or closing the crack of the 

same severity [54]. This reveals the similarity between crack extension from node B to C 

along one element and crack closure at node C (Fig. 4), which is the case of only Mode I 

loading. 

C

C' 

BA

A'  

ΔaΔc

x

y

B' 

C

C' 

BA

A'  

ΔaΔc

x

y

B' 

a b

 

Fig. 4 Mesh configuration used for VCCT element, which describes the growth process of 

the crack. Crack closure at points (a) C and (b) B 

Fig. 5 shows a crack inclination and crack tip with respect to the axes in the local 

coordinate system. 

θC
C' 

B

A

A'  

Crack 

element

Δa

Δc

X

Y

x

y

 

Fig. 5 Inclined crack in a plane with specific coordinates system. The crack inclination and 

crack tip in the local coordinate system are displayed 

The details of VCCT elements can be found in Refs. [45, 62]. In particular, each element 

has five nodes. Nodes C and C' are located at the crack tip, node B is located ahead of the 

crack tip, and nodes A and A' are located behind the crack tip. The present element has two 

groups of nodes: a top group (nodes C, A, and B) and a bottom group (nodes C' and A'). A 

very stiff spring is placed between nodes C and C', and the crack-tip nodal forces related 

with node displacement are given as follows: 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 7 

  Cx Cx C CF K u u   ,  Cy Cy C CF K v v    (7) 

where (uC, vC) and (uC’, vC’) represent the displacement of nodes C and C' at the specific 

coordinate system (X, Y), respectively; KCx and KCy
 
represent the X- and Y-direction 

stiffness, respectively. Formerly, they are established with large values [45, 62], but when 

the crack growth is predicted, they reduce to zero. 

Dummy nodes A, A', B are only employed to record the information on the displacement of 

the opening node behind the crack tip and the crack jump length before the tip, and do not 

contribute to the stiffness matrix. For nodes A and A', the displacement openings are 

calculated as: 

 
A A

u u u    , A Av v v     (8) 

where (uA, vA) and (uA’, vA’) represent the displacement of nodes A and A' under the specific 

coordinate system (X, Y), respectively. Hence, the crack jump length is the distance between 

nodes C and B, whose relation is given as follows: 

    
2 2

B C B C
a x x y y      (9) 

where (xC, yC) and (xB, yB) are the coordinate values of nodes C and B in the global 

coordinate system, respectively. 

The strain energy release rates (GI and GII) under the specific coordinate system (x, y) 

attached to the crack tip was computed. Thus, the fracture modes I and II can be separated 

(Fig. 5). The angle between X and x is given as: 

 arccos arcsin
B C B C

x x y y

a a


    
    

    
 (10) 

By projecting the corresponding nodal forces and crack opening displacement into 

specific coordinates, we can obtain: 

 
cos sin

sin cos

CxCx

CyCy

FF

FF

 

 

    
         

 (11) 

and 

 
' '

' '

cos sin

sin cos

AA AA

AA AA

u u

v v

     
    

     

 

 
 (12) 

Considering 2D-VCCT, we can express the approximate SERRs as follows: 

 
I

2

Cy
F v

G
B a





, 

II
2

Cx
F u

G
B a





, (13) 

where B represents the structure thickness. 

For the analysis of the dynamic loading fixed crack problem, we can describe the 

relationship between dynamic stress intensity factors KI (KII) and strain energy release rates 

GI (GII) as follows: 



8 Y. CAI, L. ZHOU, M. LI 

  I ICyK sign F EG  (14) 

  II IICxK sign F EG  (15) 

where E  represents material matrix and has different expressions in plane stress and plane 

strain conditions (plane stress 
tip

E E , plane strain condition: 

 
2

1

tip

tip

E
E 

 

); Etip and 

μtip denote the crack tip elasticity modulus and crack tip Poisson’s ratio, respectively [63]. 

3. RESULTS AND DISCUSSION 

3.1 Example 1 

Given the problem of central enamel cracks, we have selected finite-width central crack 

tension trials. The geometric model and the loading are shown in Fig. 6, where L, W, and B 

represent length, width, and thickness of the model, respectively (L=2 mm, W=2 mm, B=0.1 

mm). The crack length was 2a and the tension load is σ=30 MPa. The enamel is composed 

of materials M and C. The two materials were subjected to Eqs. (1) and (2), where the 

power exponent factors were n=0, 0.5, 1.0, 5.0 and the elastic moduli were Mm=90 GPa and 

Mc=120 GPa. These two materials had the same Poisson's ratio ν =0.3 (the properties of the 

enamel are the same in other examples). 

L

σ 

σ 

W

2a

 
Fig. 6 Geometric model of an enamel specimen containing a central crack 

To illustrate the convergence of this method, we have simulated a discrete model of 4 

element meshes to the normalized DSIFs KI (40×40 elements, 80×80 elements, 120×120 

elements). The model was geometrically established with 2D CPE4 (2D standard plane 

strain elements) in ABAQUS®. We have employed no special means, such as the collapsed 

element technique at the crack tip or special singular elements. Table 1 shows energy 

release rate GI from the computation model with power exponent factors n=0.0, 0.5, 1.0, 

and 5.0, four grids, and crack lengths a=0.2, 0.3, 0.4, 0.5, and 0.6 mm. By comparing with 

the analytical solutions, we validated the accuracy and effectiveness of the new method. 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 9 

Fig. 7 shows energy release rate GI computed from the model with power exponent 

factors n=0.0, 0.5, 1.0, and 5.0, crack lengths a=0.3, 0.4, 0.5, and 0.6 mm, and cell division 

= 100×100. Compared with the analytical solutions, the maximum error was only 3.2%, 

which indicates that the nonhomogeneous enamel calculation model and the calculation 

method proposed in this study is accurate and effective. 

Table 1 Energy release rate GI (N/mm) computed from the new model at power exponent 

factors of n=0.0, 0.5, 1.0, and 5.0 

n Method 
a (mm) 

0.2 0.3 0.4 0.5 0.6 

0.0 

VCCT 120120 0.00485 0.00792 0.01212 0.01856 0.02664 

VCCT 8080 0.00484 0.00791 0.01202 0.01863 0.02687 

VCCT 4040 0.00480 0.00790 0.01208 0.01838 0.02681 

Anal. solu. 0.00496 0.00807 0.01235 0.01896 0.02731 

0.5 

VCCT 120120 0.00535 0.00854 0.01301 0.02006 0.02887 

VCCT 8080 0.00526 0.00852 0.01303 0.01999 0.02901 

VCCT 4040 0.00524 0.00851 0.01297 0.02016 0.02898 

Anal. solu. 0.00535 0.00871 0.01333 0.02046 0.02947 

1.0 

VCCT 120120 0.00552 0.00901 0.01376 0.02116 0.03056 

VCCT 8080 0.00555 0.00903 0.01372 0.02113 0.03053 

VCCT 4040 0.00553 0.00904 0.01381 0.02130 0.03062 

Anal. solu. 0.00566 0.00922 0.01411 0.02167 0.03121 

5.0 

VCCT 120120 0.00637 0.01062 0.01598 0.02435 0.03511 

VCCT 8080 0.00633 0.01043 0.01596 0.02433 0.03503 

VCCT 4040 0.00637 0.01039 0.01599 0.02448 0.03532 

Anal. solu. 0.00651 0.01065 0.01632 0.02504 0.03608 

0.3 0.4 0.5 0.6

0.010

0.015

0.020

0.025

0.030

0.035

G
I 
( 

N
/m

m
)

a (mm)

 Anal. solu. n=0.0

 Anal. solu. n=0.5

 Anal. solu. n=1.0

 Anal. solu. n=5.0

 VCCT n=0.0

 VCCT n=0.5

 VCCT n=1.0

 VCCT n=5.0

 

Fig. 7 Energy release rate GI at crack lengths a= 0.3, 0.4, 0.5, and 0.6 mm with different 

power exponent factors n=0.0, 0.5, 1.0, and 5.0 



10 Y. CAI, L. ZHOU, M. LI 

3.2 Example 2 

Given the three-point bending problems of enamel, we have conducted three-point 

bending trials. The geometric model with length W=4.4 mm, width H=1 mm, span S=4 mm, 

thickness B=0.25 mm is shown in Fig. 8, where a denotes the length of the crack. The model 

is subject to a concentrated load P of 1.0 N. 

S

P

a

W

H

 

Fig. 8 Enamel specimens for three-point bending trials 

Fig. 9 shows energy release rate GI with power exponent factors n=0.0, 0.5, 1.0, and 5.0 

and crack lengths a=0.45, 0.50, 0.55, and 0.6 mm. Compared with the analytical solutions, 

the maximum error was only 3.9%, indicating that the nonhomogeneous enamel 

computation model and the method proposed are accurate and effective. 

0.45 0.50 0.55 0.60

0.010

0.015

0.020

0.025

0.030

0.035

G
I 
( 

N
/m

m
)

a (mm)

 Anal. solu. n=0.0

 Anal. solu. n=0.5

 Anal. solu. n=1.0

 Anal. solu. n=5.0

 VCCT n=0.0

 VCCT n=0.5

 VCCT n=1.0

 VCCT n=5.0

 

Fig. 9 Energy release rate GI of three-point bending enamel specimens at crack lengths 

a=0.45, 0.50, 0.55, and 0.6 mm with different power exponent factors n=0.0, 0.5, 1.0, and 

5.0 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 11 

3.3 Example 3 

In the enamel compact tensile trials, the geometric model and loading are shown in Fig. 

10, where L, W, and B represent length, width, and thickness of the model, respectively (L=1 

mm, W=1.25 mm, B=0.1 mm). The crack length is a and P=1.0 N is a concentrated load. 

L

P

a

P
W

 

Fig. 10 Enamel specimen for compact tensile trials 

Fig. 11 demonstrates energy release rate GI computed from power exponent factors 

n=0.0, 0.5, 1.0, and 5.0, and crack lengths a=0.3, 0.4, 0.5, and 0.6 mm. Compared with the 

analytical solutions, the maximum error was only 3.3%, indicating that the 

nonhomogeneous enamel computation model and the method proposed are accurate and 

effective. 

0.3 0.4 0.5 0.6

0.010

0.015

0.020

0.025

0.030

0.035

0.040

0.045

G
I 
( 

N
/m

m
)

a (mm)

 Anal. solu. n=0.0

 Anal. solu. n=0.5

 Anal. solu. n=1.0

 Anal. solu. n=5.0

 VCCT n=0.0

 VCCT n=0.5

 VCCT n=1.0

 VCCT n=5.0

 

Fig. 11 Energy release rate GI from compact tension enamel specimens at crack lengths 

a=0.3, 0.4, 0.5, and 0.6 mm with different power exponent factors n=0.0, 0.5, 1.0, and 5.0 



12 Y. CAI, L. ZHOU, M. LI 

3.4 Example 4 

The geometric model and loading for double-edge-notched tension tests are shown in 

Fig. 12, where L, W, and B represent length, width, and thickness of the model, respectively 

(L=2 mm, W=2 mm, B=0.1 mm). The crack length is a, and the load σ=40 MPa is uniformly 

distributed. 

L

σ 

a a

σ 

W

 

Fig. 12 Enamel specimens for double-edge-notched tension tests 

Fig. 13 shows energy release rate GI computed at power exponent factors n=0.0, 0.5, 

1.0, and 5.0 and crack lengths a= 0.2, 0.3, 0.4, and 0.5 mm. Compared with the analytical 

solutions, the maximum error was only 4.1%, indicating that the nonhomogeneous enamel 

computation model and the method proposed are accurate and effective. 

0.2 0.3 0.4 0.5

0.010

0.015

0.020

0.025

0.030

0.035

0.040
 Anal. solu. n=0.0

 Anal. solu. n=0.5

 Anal. solu. n=1.0

 Anal. solu. n=5.0

 VCCT n=0.0

 VCCT n=0.5

 VCCT n=1.0

 VCCT n=5.0

G
I 
( 

N
/m

m
)

a (mm)  

Fig. 13 Energy release rate GI from compact tension enamel specimens at crack lengths a= 

0.2, 0.3, 0.4, and 0.5 mm 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 13 

4. CONCLUSIONS 

Based on FEM, the 2D-VCCT enamel dumb element was constructed, which can be 

employed in the fracture analysis of enamel. We implemented the element into ABAQUS 

by the user element subroutine (UMAT and UEL). In view of the numerical examples, the 

results were in agreement with the analytical solutions, which verified that this proposed 

interface element can directly apply fracture mechanics to solve fracture problems on 

ABAQUS. The proposed element has some advantages, such as the ability to extract 

fracture parameters without extra post-processing while the definition of body mesh is 

unconstrained. The proposed element can also be integrated with FEM. Thus, as the 

reliability of the proposed element was proved in various cases, it can be potentially 

employed in nonhomogeneous material fractures. In summary, the proposed method is 

simple, efficient, universal, and steady, and can be applied to the structure-level analysis of 

nonhomogeneous fracturing problems by resorting to the commercial FEA codes. 

Acknowledgments：The authors are grateful for the financial support to the National Natural 
Science Foundation of China (grant no. 51975243), Scientific Research Project of Jilin Provincial 

Department of Education (grant no. JJKH20190131KJ) and Graduate Innovation Fund of Jilin 

University (grant no. 101832020CX099). 

REFERENCES 

1. Schwartz, G.T., 2000, Enamel thickness and the helicoidal wear plane in modern human mandibular 

molars, Archives of Oral Biology, 45(5), pp. 401-409. 

2. Sajewicz, E., 2006, On evaluation of wear resistance of tooth enamel and dental materials, Wear, 

260(11-12), pp. 1256-1261. 

3. Baino, F., Verne, E., 2017, Production and Characterization of Glass-Ceramic Materials for Potential 

Use in Dental Applications: Thermal and Mechanical Properties, Microstructure, a nd In Vitro 

Bioactivity, Applied Sciences, 7(12), 1330. 

4. Zhao, X.L., O'Brien S., Shaw, J., Abbott, P., Munroe, P., Habibi, D., Xie, Z.H., 2013, The origin of 

remarkable resilience of human tooth enamel, Applied Physics Letters, 103(24), 241901. 

5. Cui, F.Z., Ge, J., 2007, New observations of the hierarchical structure of human enamel, from nanoscale 

to microscale, Journal of Tissue Engineering and Regenerative Medicine, 1(3), pp. 185-191. 

6. Nanci, A., 2008, Ten Cate's Oral Histology: Development, Structure, and Function, 8th ed, Elsevier, 

Missouri, US. 

7. He, L.H., Swain M.V., 2008, Understanding the mechanical behaviour of human enamel from its 

structural and compositional characteristics, Journal of the Mechanical Behavior of Biomedical 

Materials, 1(1), pp. 18-29. 

8. Gao, X., Qin, W., Wang, P., Wang, L., Weir, M., Reynolds, M.A., Zhao, L., Lin, Z.M., Xu, H., 2019, 

Nano-Structured Demineralized Human Dentin Matrix to Enhance Bone and Dental Repair and 

Regeneration, Applied Sciences, 9(5), 1013. 

9. Lubisich, E.B., Hilton, T.J., Ferrancane, J., Perecdent, N., 2010, Cracked Teeth: A Review of the 

Literature, Journal of Esthetic and Restorative Dentistry, 22(3), pp. 158-167. 

10. Chai, H., Lee, J.J.W., Constantino, P.J., Lucas, P.W., Lawn, B.R., 2009, Remarkable resilience of teeth, 

Proceedings of the National Academy of Sciences of the United States of America, 106(18), pp. 

7289-7293. 

11. Chai, H., Lee, J.J.W., Kwon, J.Y., Lucas, P.W., Lawn, B.R., 2009, A simple model for enamel fracture 

from margin cracks, Acta Biomaterialia, 5(5), pp. 1663-1667. 

12. Lawn, B.R., Lee, J.J.W., 2009, Analysis of fracture and deformation modes in teeth subjected to occlusal 

loading, Acta Biomaterialia, 5(6), pp. 2213-2221. 



14 Y. CAI, L. ZHOU, M. LI 

13. Kato, H., Taguchi, Y., Imai, K., Ruan, Y., Tsai, Y.W., Chen, Y.C., Shida, M., Taguchi, R., Tominaga, K., 

Umeda, M., 2018, The enhancing effects of amelogenin exon 5-encoded peptide from enamel matrix 

derivative on odontoblast-like KN-3 cells, Applied Sciences, 8(10), 1890. 

14. Lara-Carrillo, E., Lovera-Roias, N., Luckie, R.A.M., Robles, L., Garcia-Fabila, M.M., Rosa-Santillana, R., 
Medina-Solis, C.E., 2018, The effects of remineralization via fluoride versus low-level laser IR810 and 

fluoride agents on the mineralization and microhardness of bovine dental enamel, Applied Sciences, 8(1), 

78. 
15. Meredith, N., Sherriff, M., Setchell, D.J., Swanson, S.A.V., 1996, Measurement of the microhardness and 

Young's modulus of human enamel and dentine using an indentation technique , Archives of Oral Biology, 

41(6), pp. 539-545. 

16. Wegest, U.G.K., Bai, H., Saiz, E., Tomsia, A.P., Ritchie, R.O., 2015, Bioinspired structural materials, 

Nature Materials, 14(1), pp. 23-36. 

17. Cuy, J.L., Mann, A.B., Livi, K.J., Teaford, M.F., Weihs, T.P., 2002, Nanoindentation mapping of the 

mechanical properties of human molar tooth enamel, Archives of Oral Biology, 47(4), pp. 281-291. 

18. Braly, A., Darnell, L.A., Mann, A.B., Teaford, M.F., Weihs, T.P., 2007, The effect of prism orientation on 

the indentation testing of human molar enamel, Archives of Oral Biology, 52(9), pp. 856-860. 

19. Park, S., Wang, D.H., Zhang, D.S., Romberg, E., Arola, D., 2008, Mechanical properties of human 

enamel as a function of age and location in the tooth, Journal of Materials Science-Materials in Medicine, 

19(6), pp. 2317-2324. 

20. Habelitz, S., Marshall, S.J., Marshall, G.W., Balooch, M., 2001, Mechanical properties of human dental 

enamel on the nanometre scale, Archives of Oral Biology, 46(2), pp. 173-183. 

21. Balooch, G., Marshall, G.W., Marshall, S.J., Warren, O.L., Asif, S.A.S., Bolooch, M., 2004, Evaluation of 

a new modulus mapping technique to investigate microstructural features of human teeth, Journal of 

Biomechanics, 37(8), pp. 1223-1232. 

22. He, L.H., Fujisawa, N., Swain, M.V., 2006, Elastic modulus and stress-strain response of human enamel 

by nano-indentation, Biomaterials, 27(24), pp. 4388-4398. 

23. Mann, A.B., Dickinson, M.E., 2006, Nanomechanics, chemistry and structure at the enamel surface, 

Monographs in Oral Science, 19, pp. 105-131. 

24. Bajaj, D., Arola, D.D., 2009, On the R-curve behavior of human tooth enamel, Biomaterials, 30(23-24), 

pp. 4037-4046. 

25. Liu, P.F., Hou, S.J., Chu, J.K., Hu, X.Y., Zhou, C.L., Liu, Y.L., Zheng, J.Y., Zhao, A., Yan, L., 2011, 

Finite element analysis of postbuckling and delamination of composite laminates using virtual crack 

closure technique, Composite Structures, 93(6), pp. 1549-1560. 

26. Sun, L., Ma, D.J., Wang, L.Z., Shi, X.Z., Wang, J.L., Chen, W., 2018, Determining indentation fracture 

toughness of ceramics by finite element method using virtual crack closure technique , Engineering 

Fracture Mechanics, 197, pp. 151-159. 

27. Sun, Z.Z., Zhang, X.Y., Zhang, Y.M., 2019, Cracking elements method for simulating complex crack 

growth, Journal of Applied and Computational Mechanics, 5(3), pp. 552-562. 

28. Zmudzki, J., Chladek, G., Krawczyk, C., 2019, Relevance of Tongue Force on Mandibular Denture 

Stabilization during Mastication, Journal of Prosthodontics-Implant Esthetic and Reconstructive 

Dentistry, 28(1), pp. E27-E33. 

29. Zmudzki, J., Chladek, G., Malara, P., 2018, Use of finite element analysis for the assessment of 

biomechanical factors related to pain sensation beneath complete dentures during mastication , Journal of 

Prosthetic Dentistry, 120(6), pp. 934-941. 

30. Zmudzki, J., Chladek, G., Panek, K., Lipinski,P., 2020, Finite element analysis of adolescent mandible 

fracture occurring during accidents, Archives of Metallurgy and Materials, 65(1), pp. 65-72. 

31. Nehdi, M.L., Ali, M.A.E.M., 2019, Experimental and numerical study of engineered cementitious 

composite with strain recovery under impact loading, Applied Sciences, 9(5), pp. 994. 

32. Shabani, S., Cunedioglu, Y., 2019, Free vibration analysis of functionally graded beams with cracks, 

Journal of Applied and Computational Mechanics, 6(4), pp. 908-919. 

33. Zenkour, A.M., Abouelregal, A.E., 2019, Fractional thermoelasticity model of a 2D problem of mode-I 

crack in a fibre-reinforced thermal environment, Journal of Applied and Computational Mechanics, 5(2), 

pp. 269-280. 

34. Cui, W., Zhang, Y.H., Xiao, Z.M., Zhang, Q., 2020, A new magnetic structural algorithm based on virtual 

crack closure technique and magnetic flux leakage testing for circumferential symmetric double -crack 

propagation of X80 oil and gas pipeline weld, Acta Mechanica, 231(3), pp. 1187-1207. 



 Fracture Parameters Evaluation for the Cracked Nonhomogeneous Enamel Based on... 15 

35. Fang, G., Li, L., Gao, X.G., Song, Y.D., 2020, Finite element analysis of the crack deflection in fiber 

reinforced ceramic matrix composites with multilayer interphase using virtual crack closure technique , 

Applied Composite Materials, 27(3), pp. 307-320. 

36. Zhang, W., Jiang, W.C., Yu, Y., Zhou, F., Luo, Y., Song, M., 2019, Fatigue crack simulation of the 316L 

brazed joint using the virtual crack closure technique, International Journal of Pressure Vessels and 

Piping, 173, pp. 20-25. 

37. Zhao, Y.,Sun, X.J., Cao, P., Ling, Y.F., Gao, Z., Zhan, Q.B., Zhou, X.J., Diao, M.S., 2019, Mechanical 

performance and numerical simulation of basalt fiber reinforced concrete (BFRC) using double -K 

Fracture model and virtual crack closure technique (VCCT), Advances in Civil Engineering, 2019, 

5630805. 

38. Chiu, T.C., Lin, H.C., 2009, Analysis of stress intensity factors for three-dimensional interface crack 

problems in electronic packages using the virtual crack closure technique, International Journal of 

Fracture, 156(1), pp. 75-96. 

39. Scheel, J., Ricoeur, A., Krupka, M., 2019, Calculation of stress intensity factors with an analytical 

enrichment of the modified crack closure integral, 25th International Conference on Fracture and 

Structural Integrity, 18, pp. 268-273. 

40. Zhao, X., Mo, Z.L., Guo, Z.Y., Li, J., 2020, A modified three-dimensional virtual crack closure technique 

for calculating stress intensity factors with arbitrarily shaped finite element mesh arrangements across 

the crack front, Theoretical and Applied Fracture Mechanics, 109, 102695. 

41. Parashar, A., Mertiny, P., 2012, Study of mode I fracture of graphene sheets using atomistic based finite 

element modeling and virtual crack closure technique, International Journal of Fracture, 176(1), pp. 

119-126. 

42. Zhang, H., Qiao, P.Z., 2020, Virtual crack closure technique in peridynamic theory, Computer Methods in 

Applied Mechanics and Engineering, 372, 113318. 

43. Heidari-Rarani, M., Sayedain, M., 2019, Sayedain, Finite element modeling strategies for 2D and 3D 

delamination propagation in composite DCB specimens using VCCT, CZM and XFEM approaches, 

Theoretical and Applied Fracture Mechanics, 103, 102246. 

44. Xie, D., Biggers, S.B., 2006, Progressive crack growth analysis using interface element based on the 

virtual crack closure technique, Finite Elements in Analysis and Design, 42(11), pp. 977-984. 

45. Ikeda, T., Sun, C.T., 2001, Stress intensity factor analysis for an interface crack between dissimilar 

isotropic materials under thermal stress, International Journal of Fracture, 111(3), pp. 229-249. 

46. Farkash, E., Banks-Sills, L., 2017, Virtual crack closure technique for an interface crack between two 

transversely isotropic materials, International Journal of Fracture, 205(2), pp. 189-202. 

47. Heidari-Rarani, M., Sayedain, M., 2019, Finite element modeling strategies for 2D and 3D delamination 

propagation in composite DCB specimens using VCCT, CZM and XFEM approaches, Theoretical and 

Applied Fracture Mechanics, 103, 5630805. 

48. Krueger, R., 2004, Virtual crack closure technique: History, approach, and applications, Applied 

Mechanics Reviews, 57(2), pp. 109-143. 

49. Leski, A., 2007, Implementation of the virtual crack closure technique in engineering FE calculations, 

Finite Elements in Analysis and Design, 43(3), pp. 261-268. 

50. Azimi, M., Mirjavadi, S.S., Asli, S.A., Hamouda, A.M.S., 2017, Fracture analysis of a special cracked 

lap shear (CLS) specimen with utilization of virtual crack closure technique (VCCT) by finite element 

methods, Journal of Failure Analysis and Prevention, 17(2), pp. 304-314. 

51. Zeng, W., Liu, G.R., Jiang, C., Dong, X.W., Chen, H.D., Bao, Y., Jiang, Y., 2016, An effective fracture 

analysis method based on the virtual crack closure-integral technique implemented in CS-FEM, Applied 

Mathematical Modelling, 40(5-6), pp. 3783-3800. 

52. Banks-Sills, L., Farkash, E., 2016, A note on the virtual crack closure technique for a bimaterial interface 

crack, International Journal of Fracture, 201(2), pp. 171-180. 

53. Muthu, N., Falzon, B.G., Maiti, S.K., Khoddam, S., 2014, Modified crack closure integral technique for 

extraction of SIFs in meshfree methods, Finite Elements in Analysis and Design, 78, pp. 25-39. 

54. Zhou, L.M., Meng, G.W., Li, x.l., Li, F., 2016, Analysis of dynamic fracture parameters in functionally 

graded material plates with cracks by graded finite element method and virtual crack closure technique , 

Advances in Materials Science and Engineering, 2016, 8085107. 

55. Santare, M.H., Thamburaj, P., Gazonas, G.A., 2003, The use of graded finite elements in the study of 

elastic wave propagation in continuously nonhomogeneous materials, International Journal of Solids and 

Structures, 40(21), pp. 5621-5634. 



16 Y. CAI, L. ZHOU, M. LI 

56. Anlas, G., Santare, M.H., Lambros, J., 2000, Numerical calculation of stress intensity factors in 

functionally graded materials, International Journal of Fracture, 104(2), pp. 131-143. 

57. Hagemann, A., Rohr, K., Stiehl, H.S., 2002, Coupling of fluid and elastic models for biomechanical 

simulations of brain deformations using FEM, Medical Image Analysis, 6(4), pp. 375-388. 

58. Kim, J.H., Paulino, G.H., 2002, Mixed-mode fracture of orthotropic functionally graded materials using 

finite elements and the modified crack closure method, Engineering Fracture Mechanics, 69(14-16), pp. 

1557-1586. 

59. Santare, M.H., Lambros, J., 2000, Use of graded finite elements to model the behavior of 

nonhomogeneous materials, Journal of Applied Mechanics-Transactions of the Asme, 67(4), pp. 819-822. 

60. Li, C.Y., Zou, Z.Z., Duan, Z.P., 2000, Multiple isoparametric finite element method for nonhomogeneous 

media, Mechanics Research Communications, 27(2), pp. 137-142. 

61. Dolbow, J.E., Gosz, M., 2002, On the computation of mixed-mode stress intensity factors in functionally 

graded materials, International Journal of Solids and Structures, 39(9),  pp. 2557-2574. 

62. Qian, Q., Xie, D., 2007, Analysis of mixed-mode dynamic crack propagation by interface element based 

on virtual crack closure technique, Engineering Fracture Mechanics, 74(5), pp. 807-814. 

63. Yang, Z.J., Yao, F., Huang, Y.J., 2020, Development of ABAQUS UEL/VUEL subroutines for scaled 

boundary finite element method for general static and dynamic stress analyses, Engineering Analysis with 

Boundary Elements, 114, pp. 58-73.