Comparison of two multiaxial fatigue models applied to dental implants

This paper presents two multiaxial fatigue life prediction models applied to a commercial dental implant. One model is called Variable Initiation Length Model and takes into account both the crack initiation and propagation phases. The second model combines the Theory of Critical Distance with a critical plane damage model to characterise the initiation and initial propagation of micro/meso cracks in the material. This paper discusses which material properties are necessary for the implementation of these models and how to obtain them in the laboratory from simple test specimens. It also describes the FE models developed for the stress/strain and stress intensity factor characterisation in the implant. The results of applying both life prediction models are compared with experimental results arising from the application of ISO-14801 standard to a commercial dental implant


INTRODUCTION
esigning real mechanical components subjected to cyclic loading is a traditionally complex matter due to the numerous variables affecting fatigue life.In addition to the material properties, the stress concentration, the component size and its surface characteristics as well as the multiaxiality of the stress state induced, influence its fatigue response.The use of prediction models allows studying this response but, taking into account the randomness of the process and that the effect of these variables on the fatigue life of different materials is not perfectly known, experimental analysis is often needed.Numerous methods have been proposed to carry out the fatigue design of components in presence of stress risers and multiaxial stress states.Some of these methods consider only the crack initiation phase [1], assuming that the length of the propagation phase is negligible compared to that of initiation.This approach usually provides good results for large notches and low stress levels.However, it is difficult to know a priori if, for a given combination of geometry, materials and load state, the above assumption is correct.In general, these methods tend to produce more or less conservative results, depending on the actual importance of the propagation phase with respect to the total life.A common aspect to these methods is the definition of the point or depth from the surface where the stresses have to be evaluated in order to apply the corresponding initiation criteria.Other methods only consider the propagation phase assuming that the duration of the initiation is small, either by the prior existence of small defects (microcracks) in the material [2], or by the assumption that these microcracks are quickly initiated during the first cycles, due to the high stress levels produced by a stress riser.This stress concentration can be D due to the geometry or to defects in the material such as inclusions [3].In these methods, if the initial crack is large enough, the propagation is analysed by LEFM by applying a propagation law for long cracks.If the crack is considered to be initiated from a microstructural defect, various models take into account the behaviour of short cracks [4][5][6][7].However, models regarding short crack growth are not sufficiently developed, and some of them do not work properly under certain conditions or with certain materials, so that there is currently no method that can be applied successfully in any circumstance.To avoid these problems, the defect is considered to be long enough so it is correct to apply the LEFM or some modification of it.The problem here is that predictions may be too conservative, with the addition that it is not possible to quantify the error.Another group of methods considers the fatigue process as the combination of an initiation and a propagation phase, analysing life as the sum of the durations of both phases [8,9].In most of these methods, the duration of the initiation phase is determined by local strain methods, through the ε-N curves, and the propagation based on fracture mechanics methods.In this paper, to characterise the fatigue life of a dental implant, two life prediction models for notched elements subjected to multiaxial stress states are used.As it will be shown, one of them analyses separately the initiation and propagation phases, and combines these results to obtain the total life of the specimen.The other life prediction model used is based on the Theory of the Critical Distance (TCD) and considers that, in presence of notches, the crack propagation phase is negligible compared to the total life of the specimen.Therefore, it analyses the initiation phase and initial propagation of micro/mezzo cracks using a combination of the TCD and a critical plane damage model.

THEORETICAL MODELS
oth life prediction models used in this paper are described in this section.

Variable initiation length model (VIL)
This paper uses a model for the prediction of life already proposed by the authors [10].It bears the characteristic of combining the initiation and propagation stages without having to predefine the crack length at which the initiation ends and propagation begins.However, each phase is analysed separately.The initiation phase is analysed by determining the number of cycles required to generate a crack.This number is calculated from the stresses along the path followed by the crack and from a fatigue curve ε-N, which will be detailed later.The result is a curve, a -Ni, representing the cycles required to cause a crack of length a.In the propagation phase, the number of cycles needed to propagate a crack from any length a up to failure, is calculated using fracture mechanics.To do this, the growth law is integrated from each crack length, a, until failure, yielding the curve (a -Np).The sum of these two curves would provide the total life depending on what value is taken for the crack length separating the initiation and propagation phases.It has been shown [10] that the initiation process dominates near surface, and the propagation process does so farther from it, so that the link between the two is found in the minimum of the total life curve described above.For this reason, and because it is the most conservative value, the minimum of the curve is taken as the total life of the specimen.

Variable initiation length model (VIL). Initiation phase
The analysis of the initiation phase uses the idea of McClung et al. [11] for notches, where they obtain an initiation curve by subtracting the propagation cycles from the total number of cycles until fracture.
The first step consists in obtaining a fatigue curve,   | t a N , in smooth unnotched specimens, which provides the number of cycles required to generate a crack of length at , as a function of the applied strain.For each level of strain, εj, the number of cycles of this curve,  j t a N , is obtained by subtracting the number of cycles to grow the crack from a t to failure from the total number of cycles.Each curve,   | t a N , for different values of at will be called initiation curve.In case of application of the model in a simple fatigue test, the number of cycles required to generate a crack of length, a t , could be calculated using the appropriate curve,   | t a N .In the case of a component with a multiaxial state and stress gradient, the same process may be applied, but with some modifications.Firstly, this requires a multiaxial fatigue criterion; in this case, Fatemi-Socie [12] will be used.Subsequently, the Fatemi-Socie parameter (FS) is calculated for each strain B level in the initiation curves,   | t a i N , obtained previously.With this, the new curves are constructed,  | t a i FS N .On the other hand, when there is a notch, stress decreases rapidly with depth, from a maximum at the surface.Therefore, the estimated initiation life will be one or another, depending on where the damage parameter used is assessed.The option considered most appropriate is to calculate the average FS between the surface and the crack length at, and with it, to enter the curve and obtain the number of cycles required to generate a crack of length at.This option implies the hypothesis that an equal value for the average damage parameter in the area will produce the same number of cycles to initiate the crack of that length.

Variable initiation length model (VIL). Propagation phase
Fracture mechanics is applied for the propagation phase, taking as initial length a generic length, a.The growth law used also attempts to model the growth of small cracks, since the defined initiation length can be in the order of microns.The way to do this is by introducing a modified growth threshold as a function of crack length [10].
where ΔK th is the long crack fatigue crack growth threshold, f is a parameter generally taken as equal to 2.5 [12], d is the typical distance to the first microstructural barrier (in this case half the grain size) and a0 is the so-called El Haddad constant [13], defined in the expression: where σFL and Kth are the amplitude of the stress at the fatigue limit and the amplitude of the SIF at the threshold, respectively, for R = −1.In this case, these values also correspond to the positive part of the stress and SIF cycle.

Variable initiation length model (VIL). Combination of initiation and propagation
Once the two previously mentioned curves have been obtained, (a -Np y a -Ni) they are both added, rendering a curve that represents the total life as a function of the value taken for initiation length.Fig. 1 depicts an example of those curves, obtained for a test analysis with the implant under study subjected to cyclic loads varying between 220 N and 22 N.The minimum of the curve is taken as the fatigue life, and the crack length for which the minimum occurs is taken as the initiation length.In the test represented in Fig. 1, the initiation length ai is close to 85 microns.It can be seen that, in this case, the initiation life is very small compared to total life.This is due to the high stress concentration which makes the cracks to initiate rapidly.This model can be compared with others where the length from which propagation is taken is defined a priori.This would be equivalent to entering the graph in Fig. 1 with a predetermined crack length a, obtaining an initiation and a propagation life.The advantage of the VIL model is that it is more conservative and there is no need to make a decision regarding when one phase ends and the other begins.

Multiaxial critical distance model
In this section a multiaxial life prediction model [15] for mechanical notched components based on the Theory of Critical Distance (TCD) is presented.The TCD is used to estimate the fatigue damage in presence of stress concentrators when the linear-elastic stress field in the material is known.Its multiaxial nature will be taken into account by the so-called Modified Wöhler Curve Method (MWCM).

Multiaxial critical distance model. Theory of the Critical Distances in notches (TCD)
The Point Method is a well known version of the TCD [16] and establishes that the fatigue limit of a notched component subjected to uniaxial stress is reached when the effective notch stress, evaluated at a distance L/2 from the bottom of the notch, equals the fatigue limit of the unnotched material: where L is the critical distance, considered as a material property [17][18][19], and Δσ 0 is the fatigue limit.However, if finite life in a notched component is to be investigated instead of its fatigue limit [20], the authors state the hypothesis that the critical distance is dependent on the life of the component by a power law in the form: where A and B are constants depending on the material and on the load (through R), which may be obtained experimentally using two calibration curves.These can be two fatigue curves for the notched and unnotched material, obtained under the same loading conditions (same R).Indeed, as shown in Fig. 2, for any given fatigue life N f,k , the range of the maximum nominal principal stress which produces that life for an uniaxial notched specimen Δσgross (referring to gross section) and for an unnotched one, Δσ1,k, can be obtained.The critical distance associated to that life, Nf,k, will be the distance at which the maximum principal stress takes the value Δσ 1,k in the notched specimen.This distance, L, can be determined from a FE analysis of this notched specimen.
Defining two fatigue lives Nf,1 and Nf,2, the critical distances L1 and L2 can be obtained.By fitting Eq. ( 4) to theses two lengths, it is easy to obtain the constants A and B in this power law: If the maximum principal stress vs. distance from the notch is obtained in a notched component, for example, by a FE model, a function representing the evolution of Δσ with r, Δσ=f(r), can be calculated.Combining this function with the fatigue curve , the fatigue life as a function of r, N f (r), can be determined: On the other hand, the fatigue life Nf is related to the critical distance according to Eq. ( 4).Therefore, combining Eq. ( 4) and ( 5), and bearing in mind that according to TCD r = L/2; it is possible to obtain an expression depending only on r, h(r)=0.This expression is non-linear but can be easily solved by iterative methods.Solving this equation, the value of the distance to evaluate the stresses, r, is obtained.Then, the Eq. ( 4) or (5) gives the fatigue life.
It is important to highlight that the approach outlined in the preceding paragraphs is valid for estimating the fatigue life of notched components subjected to uniaxial cyclic loading.However, the actual mechanical and structural elements are usually subjected to loads that generate multiaxial stress states.In these cases it is useful to combine the TCD with some critical plane model to take into account that multiaxiality.

Multiaxial critical distance model. Modified Wöhler Curve Method (MWCM)
The Modified Wöhler Curve Method [15] uses a critical plane damage model so the effect of the different components of the stress on the crack initiation is quantified.In this model, the critical plane is the one where the amplitude of the tangential stress reaches its maximum value.This critical plane model can be used to define the multiaxiality index ρ eff , a parameter that can be calculated, knowing the stress state at each point, as: where σn,m and σn;a are, respectively, the mean value and the amplitude of the stress normal to the critical plane and τa is the amplitude of the tangential stress in this plane.The parameter m is the mean stress sensitivity index [21], a material property that can be determined experimentally.In this paper we have assumed m=1, which is the worst case scenario.This indicates that the material has a high sensitivity to the presence of superimposed static loads.The multiaxiality index ρ eff allows modifying the Wohler curve (curve S-N) of the material in terms of the shear stress amplitude (τa-N curve).To do this, experimental studies [22,23] suggest that the slope of the curve, kτ, as well as the fatigue limit, τA,Ref may be linearly related to the multiaxiality index ρeff: where a, b, α and β are constants which depend on the material and can be determined experimentally from two calibration curves, for example an uniaxial fatigue curve and a pure torsion fatigue curve, cases where the multiaxiality index ρ eff takes the values 0 and 1 respectively.Once kτ(ρeff) and τA,Ref(ρeff) are known, the specimen fatigue life can be obtained as: ,R ( ) Thus, once the stress state at a depth equal to the critical distance is known, and using τa, the number of cycles to failure N f is obtained directly.This idea is illustrated in Fig. 3.

Multiaxial critical distance model. TCD and MWCM combination: Life estimation in multiaxial notches
The combined application of the methodologies discussed above can be accomplished by the following steps: 1. First, it will be necessary to identify the crack initiation point, which is the point of maximum stress concentration caused by the notch.2. From this point and following a straight line perpendicular to the surface of the component, the potential crack path is established, and described by the coordinate r. 3. The stress state must be determined along the potential path of the crack.This is usually achieved by a linear elastic analysis based on a FE model of the system under study.4. From the stress state, the evolution of the critical plane along the crack path can be determined, as well as the different components of the stresses tangential and perpendicular to that plane, τa(r), σn,m(r) and σn,a(r ). 5. The multiaxiality index ρeff(r) along the crack path can be obtained using the stress components calculated previously.6.For every value of ρeff a modified Wöhler curve (τa-N) can be calculated using the expressions (7) for each point, r, along the crack path.7. Once these fatigue curves, τ a -N, and the tangential stress amplitude τ a (r) are known for every point in the crack path, r, the number of cycles to failure Nf(r) associated to each of those points is obtained.8.According to the Point Method, to estimate the fatigue life in a notch, stresses have to be evaluated at a distance of L/2 from the bottom of the notch, i.e., r = L/2, where L is given by Eq. ( 4).Therefore, substituting into this expression the relationship obtained in the previous point, Nf(r), the value r where the fatigue life will be evaluated is obtained (see Eq. 5).9. Once r is known, the life associated to that r is the fatigue life.

NUMERICAL MODELS
s mentioned above, the prediction models here described, were used to estimate the fatigue life of a dental implant when tested under the conditions imposed by ISO 14801 [24].The FE models made faithfully reproduce the geometry of the elements involved in this type of test, as well as the boundary conditions they are subjected to, which are shown in Fig. 4.This figure also shows the area where the crack will develop.The FE model used to characterise the stress state in the implant is shown in Fig. 5.In this model, 10-node tetrahedral elements (SOLID 187 ANSYS) and 20-node hexahedral elements (SOLID 186 ANSYS) were used.The contact between the different system components was modelled using a "bounded" type contact.Conditions of null displacement were applied to the external surface of the fixed support.As shown on the right side of the figure, the implant body has been divided into different volumes to have a better control over the mesh.Thus, in the area where the crack will develop, stress convergence lower than a 4% has been achieved, with an element size of 6 microns.The material behaviour of the model close to the crack initiation zone has been defined as elastoplastic, the rest being elastic.This is useful because the TCD model uses elastic stresses but the VIL takes plasticity into account to study the A crack initiation.At the boundary between the two different behaviour zones, the stress level was within the elastic regime, so that there were no discontinuities in stresses.Plasticity was modelled using kinematic hardening, with properties obtained from tests conducted in the laboratory.The model shown is taken as a starting point for a second one, described in more detail in [26], used to calculate the SIF and therefore study the propagation phase in the VIL model.In order to calculate the SIF, the crack is assumed to initiate at the bottom of the external thread, due to the stress concentration existing in this zone, Fig. 6.As can be seen in Fig. 8, it is considered that the crack extends along the bottom of the thread as it propagates trough the interior of the implant, along a helical surface, whose axis coincides with the axis of the implant.The central point of the crack front is supposed to follow a straight line, which obviously belongs to the propagation plane, perpendicular to the axis of the implant.When the crack is small, shorter than 100 microns, it is modeled through ANSYS WORKBENCH assuming a semielliptical flat crack inclined 11.4º, since the crack grows along the valley of thread.For larger cracks, where it can no longer be assumed that it is flat, its geometry is introduced in the solid model and defined by orthogonally projecting an ellipse onto the helical propagation surface, Fig. 8. Therefore, the parameters defining the ellipse (minor axis a and major axis b) also characterise the real crack.The center of the ellipse is located at the initiation point and its minor axis is collinear to the potential propagation path of the central point of the crack.The major axis of the ellipse is tangent to the propagation plane.For these larger cracks, submodelling was used since the number of elements was too large.A crack, with the previously described characteristics, is inserted in the FE model of the implant system and the J-Integral method is used to estimate the SIF at the crack front, Fig. 9.The mesh quality as well as the number of integration contours in the J-Integral calculation chosen ensure a 5% convergence during the SIF calculation.The material behaviour is supposed to be linear elastic.Several simulations have been performed in order to estimate the SIF evolution as a function of the crack length a, as the crack extends through the implant thickness.An initial semicircular crack of length 5 microns has been assumed.Later the advance of the crack at the surface (S) and at the deepest point (D), is considered to be controlled by the Paris law and the SIF at each point.Under these conditions, the following relation is obtained Therefore, based on the SIF obtained at a certain crack length, a given increment of  D a in the crack size at the crack front is imposed.To determine the crack aspect ratio at this new length, and hence the SIF, the following relation is applied With this value,  S a , the new crack length at the surface point is obtained and a new FE problem can be solved to obtain the SIF.This process is repeated until the final crack length is reached.Under the previously exposed propagation conditions, as can be seen in Fig. 10, the aspect ratio decreases strongly at the beginning of the crack growth.This shows that the almost circular crack grows faster at the surface than inside at the central point.This is due mainly to the effect of the stress concentration induced by the external thread of the implant body.The stress concentration always affects the crack at the surface, but its effect tends to disappear at the crack center as the crack grows.It can also be observed that the aspect ratio tends to an almost constant value as the crack length increases.Fig. 11 shows the SIF vs. crack length at the deepest point of the crack.The high stress gradient close to the surface makes the SIF increase rapidly at the beginning and at a much lower pace afterwards.Therefore, more simulations had to be performed at small crack lengths to get a better estimation of the SIF evolution with the crack length.
The SIF calculated is in mode I in the local reference system of the crack.Therefore, at the centre of the crack there really is a mixed mode growth of mode I and mode III because the crack is inclined an angle of 11.4º.Nevertheless, this angle is so small that life estimation calculations have been performed using only mode I because the contribution of mode III can been neglected, as can be observed in Fig. 12.

LIFE ESTIMATION AND COMPARISON OF RESULTS
inally, a comparison of the lives estimated by applying the TCD and VIL models, together with the results obtained experimentally is shown in Fig. 13.From a general point of view it can be said that both models have a reasonable behaviour.As can be seen in Fig. 13, the VIL model presents a better performance in estimating the fatigue life of the dental implant system, being able to correctly reproduce the slope of the fatigue curve and the presence of fatigue limit, although in this case the model seems to underestimate this limit.The latter is a conservative result, while it has to be experimentally confirmed that the actual fatigue limit is located where the only trial that has been able to withstand 5 million cycles suggests.It is important to highlight that the entire process of crack initiation and propagation takes place within about half a millimetre, and that the model is able to reflect this in an acceptable manner.On the other hand, the TCD model overestimates the life and not collects properly the slope of the fatigue curve.Fig. 14 shows a comparison between the initiation length obtained using the VIL model and the critical distance calculated through the TCD.Although the values are similar, the meaning of these parameters is very different.The application of both models to other situations should show if this similarity is a coincidence or if there is a relation between them.As a result of the VIL model application, it has also been observed that in the short/medium fatigue regime (lower than 10 5 cycles), the duration of the crack initiation phase is negligible compared to the total fatigue life.This is due to the severe notch existing at the thread root of the implant body, which causes the crack to initiate within the first loading cycles.

CONCLUSIONS
he results of this study show that the VIL life prediction model, also used in various other situations, is very versatile and robust, adapting to different circumstances.It is important to emphasize that this model fits the life prediction of implants, elements which, given their small size, can present problems of scale.Once fatigue and fracture properties of the material are known, using this model allows estimating life without presupposing the relative importance of the stages of initiation and crack propagation.The applicability of each model can be compared.Regarding the material properties needed, the VIL model requires an elastic-plastic FE analysis to determine the stresses and the strains.In this FE model, for example, a kinematic hardening curve is needed to describe the plastic behaviour of the material.Also, the crack growth properties (Paris parameters) of F T the material have to be known for SIF calculation.Finally, the uniaxial fatigue curve of the material is needed to obtain the initiation curves.All of these material properties can be obtained in the laboratory by testing simple geometry specimens, or found easily in the literature for the most commonly used materials.In the case of the TCD model, tests are needed to obtain the material fatigue curves under traction and torsion conditions and in the presence of a notch.These are not complicated tests, but the curves are not so easily found in the literature.Regarding the calculations, the VIL model requires a longer time, first to obtain the initiation curves, but majorly to characterise the stress intensity factor.

Figure 1 :
Figure 1: Application of the prediction model in the test with F = 220 N.

Figure 2 :
Figure 2: Critical distance determination using two fatigue calibration curves.

Figure 4 :
Figure 4: Dental implant system geometry and testing conditions.

Fig. 6
Fig.6shows the von Mises stress distribution in the threaded zone in a section of the implant.This figure shows the high stress concentration at the outer thread root, where the crack initiates.

Figure 5 :
Figure 5: FE model of the implant.Figure 6: Von Mises stress distribution in the implant body.

Figure 6 :
Figure 5: FE model of the implant.Figure 6: Von Mises stress distribution in the implant body.

Fig. 7
Fig.7shows the linear and non-linear solutions of the model at the maximum and minimum of a cycle for a test with Fmax = 220 N and R = 0.1.The magnitude represented is the evolution of the stress normal to the crack along the crack's path.A plastification at the bottom of the external thread can be observed, as well as the influence of the stress raiser which reaches a depth of about 100 microns.The model shown is taken as a starting point for a second one, described in more detail in[26], used to calculate the SIF and therefore study the propagation phase in the VIL model.In order to calculate the SIF, the crack is assumed to initiate at the bottom of the external thread, due to the stress concentration existing in this zone, Fig.6.As can be seen in Fig.8, it is considered that the crack extends along the bottom of the thread as it propagates trough the interior of the implant, along a helical surface, whose axis coincides with the axis of the implant.The central point of the crack front is supposed to follow a straight line, which obviously belongs to the propagation plane, perpendicular to the axis of the implant.

Figure 7 :
Figure 7: Normal stress to the crack at the maximum and minimum of the cycle (F=220N, R=0.1).Figure8: Scheme of the crack geometry in the implant.

Figure 8 :
Figure 7: Normal stress to the crack at the maximum and minimum of the cycle (F=220N, R=0.1).Figure8: Scheme of the crack geometry in the implant.

Figure 9 :
Figure 9: A section of the implant showing one the cracks modeled.

Figure 10 :
Figure 10: Evolution of the crack aspect ratio when it grows.Figure 11: Range of the stress intensity factor vs. crack length in the test with F = 220N.

Figure 11 :
Figure 10: Evolution of the crack aspect ratio when it grows.Figure 11: Range of the stress intensity factor vs. crack length in the test with F = 220N.

Figure 12 :
Figure 12: Mode I, II and III stress intensity factors along the crack front.

Figure 13 :
Figure 13: Fatigue tests in implants and theoretical estimates.Figure 14: Initiation length and critical distance vs. estimated life.

Figure 14 :
Figure 13: Fatigue tests in implants and theoretical estimates.Figure 14: Initiation length and critical distance vs. estimated life.
NOMENCLATURE a = crack length, minor semiaxis of the crack af = crack length at failure a0 = El Haddad constant b = major semiaxis of the crack f = parameter in a approximation to Kitagawa-Takahashi diagram d = distance from the surface to the first microstructural barrier L = critical distance n = exponent in Paris law da/dN = crack growth rate C = coefficient in Paris law E = Young modulus F = force applied to the implant j S f N = number of cycles to failure in a plain fatigue test with stress Sj i j a S N = number of cycles to generate a crack of length a i in a plain fatigue test with stress S j Ra = surface rougness K = stress intensity factor range KD = stress intensity factor range at the deepest point of the crack KS = stress intensity factor range at the surface K th∞ = fatigue crack growth threshold for long cracks σ FL = fatigue limit σy = yield strength σ u = tensile strength