Preventing stress singularities in peri-implant bone – a finite element analysis using a graded bone model

Abstract In finite element analysis bone is often treated as two-layered material that has a discontinuity between the cortical and cancellous bone, which leads to a singularity and incorrect stresses. The goal of this study was to eliminate this singularity and to create a more realistic representation of bone which also considers the transition zone between cortical and cancellous bone as observed in natural bone. This was achieved by modelling bone as a graded material and inserting node-specific values for Young’s modulus in the finite element simulation, whereas the transition zone thickness was derived from a CT scan. The modelling was performed semi-automatically, and the maximum principal stresses of the new approach were compared to those of a conventional approach. The new approach was found to effectively avoid singularities and provides more accurate predictions of stress in areas of the bone transition zone. As the approach is automatable and causes rather small overhead it is recommended for use in future work, when the problem at hand requires evaluating stresses close to the former singularity.


Introduction
Nowadays implants are commonly used to treat patients with one or more missing teeth.Although implant treatment has been practised for several years, complications after a successful surgery can still occur due to various biological or biomechanical effects.The main reasons for implant failure are a lack of osseointegration, peri-implantitis, overloading and implant fracture (Manor et al. 2009).While those issues may be a consequence of poor oral hygiene, some complications are associated with biomechanical issues, especially incorrect loading of the bone (Misch 2015a).During mastication, the resulting force is transferred directly from the implant to the bone.Compared to a natural tooth, an implant has a different contact surface with the bone and different material properties.In addition, it lacks the periodontal ligament.This results in altered forces and stresses in the bone (Misch 2015a).According to a bone remodelling theory (Frost 2003), loading outside the adapted window leads to bone loss and possibly an increased probability of implant failure over time.Consequently, finding an optimal biomechanical design is of great interest to dental researchers.
The FEM allows the determination of strains and stresses in both the implant and the surrounding bone and has become a recognised tool for conducting in silico studies (Lisiak-Myszke et al. 2020;Al-Kordy and Al-Saadi 2022;Alberto et al. 2022;Babaei et al. 2022).Despite its growing reputation, issues such as model validation, mesh quality and verification are not always adequately addressed (Viceconti et al. 2005;Burkhart et al. 2013).Though, addressing those issues is crucial to obtain an accurate numerical solution as well as proper model approximation (Westermann 2021).
In literature, finite element models usually consists of a varying detailed implant system and a section of the mandibular or maxillary jaw containing the implant on which a load is applied (Oswal et al. 2016;Asgharzadeh Shirazi et al. 2017;Jarrahi et al. 2018;L opez-Campos et al. 2018;Udomsawat et al. 2019;Didier et al. 2020).To distinguish the elastic behaviour of cortical and cancellous bone, a different material value for Young's modulus is applied in the simulation tool to both bone tissues.However, bone consists of several phases and has a heterogeneous and anisotropic structure (Doblar e et al. 2004).
Clinical studies show the existence of a transitional zone between cortical and cancellous bone, suggesting that there are no abrupt changes between the bone tissues (Keshawarz and Recker 1984;Zebaze et al. 2013).Nevertheless, the material properties are considered isotropic in many finite element studies, although this assumption does not correspond to reality and thus the results may only be evaluated qualitatively (Trivedi 2014).
Furthermore, this idealization leads to an abrupt material transition between the tissues causing unwanted and sometimes problematic effects like stress singularities, as noted in several publications (Petrie and Williams 2005;Schwitalla et al. 2015;Brune et al. 2019).A stress singularity occurs at a location, where stress does not converge, but tends to infinity with mesh refinement.There are several reasons for that behaviour, such as point loads, corners with no radius or dissimilar material joints (Wood et al. 2015).The origin of such singularities lies in the simplification during model building.For example, infinite sharp edges are not possible in reality or are prevented by plastic deformation.Therefore, an option to prevent singularities is to enhance the model by decreasing the grade of idealisation.
Another way to deal with singularities is to ignore the singularity, if it lies outside the area of interest.This approach can work, since according to the Saint Venant's Principle, the further away the stresses are from the singularity, the less they are affected by it (Mises 1945).However, this approach should be avoided, when stresses close to the singularity are examined.
In the present study the discontinuous transition between cortical and cancellous bone that is often chosen in bone models was replaced with an automatically generated graded transition between cortical and cancellous bone.The aim was to compare the graded and conventional approach regarding convergence behaviour, the resulting stress distribution and to highlight key aspects of the automation.Taking the heterogeneity and anisotropy into consideration when creating the finite element model could also result in a better representation of bone.

Materials and methods
In this study, a finite element analysis was conducted with two different generic implant systems and two different bone modelling approaches: a graded and a conventional transition between cortical and cancellous bone.The maximum principal stresses in the peri-implant bone were computed for both approaches and compared.Additionally a mesh convergence study was carried out to provide information on convergence behaviour.

Geometry and material properties
The geometric model for this finite element analysis consisted of an implant system, which was inserted into a left mandibular jaw section.The implant system was comprised of an implant, abutment, abutment screw and a crown and was modelled using Rhinoceros 3D (Version 7 SR15, Robert McNeel & Associates, Seattle, WA, USA) and Autodesk Inventor Professional (Version 2022.3, Autodesk, Mill Valley, CA, USA).Two types of bone level implants were modelled, one with a generic v-thread and one with no thread, both designs with a conical tip (Figure 1).The bone geometry is based on previous work (Rand et al. 2017;Brune et al. 2019) and consisted of the cancellous bone and a 1:5 mm thick layer of cortical bone on the outside (Figure 1c).Additional to the current bone geometry, a new graded bone was modelled with the same shape.This approach is described in more detail in Section 2.2.Each implant was inserted into both bone models, which resulted in a total of four different simulation models.The material properties of the implant and the conventional bone model were modelled as isotropic and homogenous and are shown in Table 1.All contacts were modelled as bonded.

Graded modelling approach
The aim of this approach was to model a transition zone between the cortical and the cancellous bone, as shown in Figure 2.This was achieved by applying a certain value for the Young's modulus to every node of the bone to create a transition zone between the two bone tissues.As a result, the bone structure became inhomogeneous and the stiffness anisotropic.In order to determine the Young's modulus, the distance to the surface was calculated for each node.This distance was used to identify the zone (cortical zone, cancellous zone, transition zone) in which the node was located.Then, a value for the Young's modulus was applied accordingly.In order to achieve this, the software Ansys (Version 2022R1, Ansys Inc., Canonsburg, PA, USA) and Matlab (Version R2021b Update 3, MathWorks, Natick, MA, USA) were used.The whole process is depicted in Figure 3.
First, all required geometries, in this case the bone and the implant, were modelled in Rhino and Inventor.The entire process can be divided into two subprocesses.In the first subprocess, the bone geometry was used to define the position of the outer hull, which was later used in the calculation of the distance.For this purpose, a self-written Matlab script was executed.In this script, the bone was imported in the .stlfile format.The built-in Matlab function pdegplot with activated FaceLabels-setting was used to visualise the previously loaded file.In this representation, all surfaces are labelled with a unique number.The user has then to determine the cortical surfaces and enter the identifying numbers of all desired surfaces in a vector.In the next step, the generateMesh function was automatically called, which creates a tetrahedral body mesh of the 3-D bone geometry.In the last step of this subprocess, points were extracted from the mesh nodes on the previously selected faces.These points represents the outer hull of the bone.
In the other subprocess, the whole model including the bone and the implant were imported into Ansys in the .stpfile format.Here, the model was set up as usual and a finite element mesh was created.The mesh of the bone was then exported to Matlab as a text file.This text file contains the node ids and their respective coordinates.For each node of the finite element mesh, the distance to the closest surface point was then computed.Depending on this distance, a value for the Young's modulus was assigned to every node in the following way (equations 1-3): With x being the perpendicular distance from the cortical surface, l cort the width of the cortical bone, l transition the width of the transition zone, E cort the defined Young's modulus for the cortical bone and E canc the defined Young's modulus for the cancellous bone.The node ids with the corresponding coordinates and the values for the Young's modulus were then reimported into Ansys and applied to the simulation.The values for E cort and E canc were the same as for the conventional modelling approach (Table 1).In this study the width of the cortical layer and the transition zone was l cort ¼ 1:355 mm and l transition ¼ 0:29 mm respectively.The determination of both values is further described in Section 2.3.
Since the surface defined for distance calculation was approximated by points, deviations from the real distance inevitably occur.For this reason, the generateMesh function mentioned above allows to pass a parameter for the mesh size, where lower values results in more accurate results but longer computational time.The maximum expected deviation error for the selected mesh and the selected parameter is presented to the user, so that he can change the mesh settings accordingly.In this case, a mesh element size was chosen between Meshsize min ¼ 0:0896 mm and Meshsize max ¼ 0:1792 mm, which resulted in a maximum deviation error of 0:44 % for nodes close to the layer edges and in the transition zone.This error indicates the maximum deviation of the calculated distance from the real distance for the mentioned nodes.
For utilizing the graded approach in Ansys, an artificial thermal load was applied to the bone body like previously described by Ashirbekop et al. (Ashirbekob et al. 2020).By disabling all thermal strain effects and setting a temperature dependent elasticity modulus, inhomogeneous material properties were simulated.The artificial temperature in this case ranged from t ¼ 0 (equivalent to a Young's modulus of E ¼ 1:37 GPa) to t ¼ 100 (equivalent to a Young's modulus of E ¼ 13:7 GPa).Temperature values between the two limits were interpolated linearly, which had the corresponding effect on the Young's modulus.In Ansys, the temperature values were applied to the mesh nodes of the bone.Since all thermal effects were turned off, this load only affected the defined elastic properties in the desired way.

Choosing the width of the transition zone
In order to chose a value for the transition zone between the cortical and cancellous bone, the software Stradview (Version 7.1, Cambridge University, Cambridge, UK) was used, which is able to measure the endocortical region from CT scans.The endocortical region exists between the compact cortical bone and the porous trabecular bone (Zebaze et al. 2013;Pearson and Treece 2018).In this study, the concept of the endocortical region was used for the transition zone.As a simplification, it is assumed, that the cortical bone and the cancellous bone have a constant density and the transition zone describes a linear decrease of density from the cortical region to the cancellous region (compare to Figure 2).
Here, the INCISIX dental CT scan from the OsiriX repository (DICOM Image Library 2022), which includes the maxilla and the mandible, was used to determine the width of the transition zone.In the first step, the mandibular jawbone was segmented.After segmentation, the area around the left first lower molar was selected by drawing a closed contour through manually placed points.Based on this contour, a new object was created and the width of the endocortical region was calculated by measuring the radio density in a straight line, starting from a defined point on the outer cortical surface and moving perpendicular to the surface inside the bone.This process was automatically performed for all points within the created surface by calling the Endo CBM v2 function.
The next step was to choose a value for the transition zone used in this study from all of this measured points.For this purpose, the mode was calculated from the measured data.The mode is the value that occurs most frequently in data set and is more insensitive to measurement errors or artefacts then the mean value.As the data is continuous, a kernel density estimation approach was chosen to calculate the mode.In this study the Matlab function ksdensity was incorporated in the automation script to return a probability density estimate with an automatically chosen bandwidth of 0.2189.The mode of the density function delivers a value of 0:291 mm and was chosen as the width for the transition zone.In order to calculate the cortical width of the graded approach, the 1:5 mm thick cortical layer of the conventional approach served as the starting point.It was assumed, that half of the transition zone covers the former cortical region.This resulted in a thickness of the cortical layer of l cort ¼ 1:355 mm:

Boundary conditions
A load was directly applied to the occlusal surface of the crown (red surface in Figure 4), with 100 N in axial direction of the implant and 20 N in mesial direction.A frictionless support was added to the lateral interfaces of the bone (blue surfaces in Figure 4), preventing any movement in normal direction.The lower face (black surface in Figure 4) was fixed, which prevented any movement.

Evaluation
The maximum principal stresses in the peri-implant bone were used to compare the conventional and graded approach.A mesh convergence study was conducted, beginning with an overall mesh element size of 0:5 mm for all bodies to see if divergence arose in the conventional models and if the graded modelling approach can prevent singularities.For this purpose, the mesh elements of the intersection surfaces between the implant and the bone were refined nine times with a refinement factor of 1.3, resulting in 10 different refinement steps for each simulation with a final mesh size of 0:047 mm (Table 2).The mesh was composed of quadratic tetrahedral elements.The highest maximum principal stresses of the nodes at the bone transition were recorded in order to enable the evaluation of the modelling approach.
An assessment of the mesh quality metrics were performed according to Burkhart et al. (Burkhart et al. 2013).The aspect ratios, the deviations from the ideal angles and the element jacobians were recorded.The requirements for 95 % of all elements were the following: An aspect ratio < 3, a maximum angle deviation of 6 70 and a jacobian ratio > 0:7: All requirements were fulfilled for the refined meshes.Furthermore it was verified, that the area of interest (the peri-implant bone near the transition) contains no failed elements.

Results
Figure 5 shows the distribution of the maximum principal stress in the peri-implant bone near the implant.The point of view is depicted in Figure 5e.For the conventional approach with a non-threaded implant (Figure 5a), the highest maximum principal stress can be found at the edge of the cortical bone near the implant and the cancellous bone.The stresses in the cortical bone were higher than in the cancellous bone.Additionally the stresses change abruptly between the two bone tissues.For the graded approach with a non-threaded implant (Figure 5b), the stresses are higher in the cortical region too.In contrast to the conventional approach, the peak stresses were overall lower and stress transitions more smoothly.
For the conventional approach with a threaded implant (Figure 5c), the highest maximum principal stress can be found at the edge of the cortical bone near the implant and the cancellous bone on the lower part of the first thread flank.The stresses in the cortical bone were higher than in the cancellous bone.Additionally, the stresses change abruptly between the two bone tissues.For the graded approach with a threaded implant (Figure 5d), the stresses are higher in the cortical region too.In contrast to the conventional approach, the peak stresses were overall lower and stress transitions more smoothly.In the transition zone, the stresses decreased continuously towards the cancellous bone.For both graded modelling approaches, the stress distribution outside the transition zone looks similar compared to the respective conventional modelling approach.
The graphs in Figure 6 show the results of the mesh convergence study.The highest maximum principal stresses in the peri-implant bone near the transition for the non-threaded implant (Figure 6a) are higher for the conventional modelling approach, due to the singularity.Furthermore, they did not converge but rose constantly with a refining mesh.The results for the graded modelling approach converged to a value of about r ¼ 4:8 MPa: Similar results were observed for the Table 2.The number of mesh elements used.Two meshes of the threaded implant with the conventional approach could not be acquired due to meshing errors.threaded implant (Figure 6b).Here, the peak stresses for the conventional modelling approach were higher than the peak stresses for graded modelling approach.
Additionally, the values did not converge either.The results for the second and the last step of the mesh refinement could not be acquired due to meshing errors.The graded modelling approach for the threaded implant shows a convergence to a value of about r ¼ 9:3 MPa (Figure 6b).The peak stresses for both modelling approaches were higher than the respective peak stresses for the non-threaded implant.

Discussion
The aim of this study was to investigate the effect of a graded modelling approach of the cortical and cancellous bone on the stresses in the peri-implant bone.
For this purpose, a finite element analysis and a mesh convergence study were conducted.The results of the mesh convergence study of the conventional modelling approach confirm the hypothesis that stress singularities can occur between contacting bodies with different material properties.Moreover, the abrupt decrease of stress at the boundary region from the cortical to the cancellous bone seems to be connected to this phenomenon.Some studies observed sudden transitions in stress as well (Petrie and Williams 2005;Schwitalla et al. 2015;Brune et al. 2019), while in other dental research articles it is unfortunately often not mentioned whether stress singularities are present and whether proper mesh verification methods like a mesh convergence study were conducted to find them.
With the graded modelling approach, it was possible to achieve a successful mesh convergence.Using this method seems appropriate when looking for an accurate solution and to prove the independence of the result from the mesh size.In numerous studies, FEM is used to investigate the effects of various parameters of dental implants on bone (Santiago et al. 2013;Oswal et al. 2016;Asgharzadeh Shirazi et al. 2017;Jarrahi et al. 2018;L opez-Campos et al. 2018;Udomsawat et al. 2019;Didier et al. 2020;Farhan-Ul-Haq et al. 2022;Onone Gialain et al. 2022).In these cases, the present method seems oversized at first glance, since qualitative statements about the resulting forces in the bone are often sufficient.However, when dealing with the bone mechanostat theory according to Frost (Frost 2003), it is essential to know exactly which stresses/strains are present in the bone to accurately predict bone remodelling.Looking at the stresses of the conventional approach with a nonthreaded implant as an example (Figure 6), it is impossible to make an exact statement about the occurring stresses, because the result varies in this case from 8:3 MPa up to 12:9 MPa: With the graded approach, on the other hand, the solution converges to a value, so that in this case it is possible to make a statement about the stresses that actually occur.
Another important finding from this study shows the influence of implant geometry on stresses in bone.
Looking at Figure 6, there are several things to notice.The stresses for the conventional approach diverge more linearly with the non-threaded implant than the stresses with the threaded implant.Furthermore, the graded approach has a huge impact on the highest maximum principal stresses on the peri-implant bone when using both implants.For the non-threaded implant, the highest maximum principal stresses decreases between 37:25 % and 53:40 % with the graded approach.When using the threaded implant, the highest maximum principal stresses decreases between 39:10 % and 47:53 %: The overall higher maximum principal stresses in the peri-implant bone of the threaded implant show that the implant geometry also has a significant influence on the stresses in the bone, which is in line with other studies (Misch 2015b;Oswal et al. 2016).Principal stresses are considered appropriate measures for brittle materials, such as bone, and are frequently utilized in biomechanical finite element studies (Yamanishi et al. 2012;Ramos Verri et al. 2015).In general, it is useful to analyse both the minimum and maximum principal stresses.In this study, we have focused on the maximum principal stresses because they are relevant in the area where the singularity occurs.
Non-threaded implants have a limited clinical benefit, but are often used as a comparison baseline in implant thread studies (Chou et al. 2008;Rismanchian et al. 2010).In this study a nonthreaded implant was used to show the robustness of the graded approach and to demonstrate that stress singularities can occur regardless of the geometry of the implant-bone interface.We expect similar results for other thread designs.
Another aspect of this study was to create the graded model with a high degree of automation as a basis for future work.It is possible to create a transition layer based on CAD files in an automated way, using parameters to control the transition and cortical thickness and the values for the Young's modulus.This not only saves modelling time, but also ensures the robustness of the model to a certain extent.A contact surface between cortical and cancellous bone is no longer required.This led to a considerably improved mesh creation and a lower computation time in Ansys.For the threaded implant, the computation times for the graded approach have decreased by 7:3 % on average, for the non-threaded implant the computation times have decreased by 12:5 % on average.In addition, no meshing errors occurred.Furthermore, this approach can also be used to change the thickness of the cortical layer without the need for using CAD modelling again.However, the automation has still potential for further research.For example, the transition and cortical thickness in this approach are identical at every location and the geometry of the bone must also be created manually.Based on CT data, both points can in principle be addressed more accurately with an increased degree of automation.
Acquiring the width of the cortical layer and transition zone by measuring CT scans with a device capable of measuring radiodensity seems to be an appropriate step.To the authors' knowledge, the endocortical region measured from CT scans was not used for a graded modelling approach in a dental finite element analysis up to now.Density measurements of the endocortical region have been used in previous studies mainly to detect changes in the bone (Lang et al. 2004;Treece and Gee 2018).
In literature, there are also other approaches that consider the graded bone structure (Rieger et al. 2018;Sheidaei et al. 2019).With data from MicroCT scans, micro finite element models can be generated that directly represent the trabecular structure of the bone.This method is suitable for determining the biomechanical properties of bones in a patient-specific manner.However, a form of data acquisition such as through a CT scan is mandatory.The current approach from this study has a more general application.
In this study, the transition zone was modelled through a graded approach.The transition zone proven in clinical studies (Keshawarz and Recker 1984;Zebaze et al. 2013) was also visible in the CT scan.It is expected that the graded approach results in a more realistic representation of the bone.In preliminary work, it was also investigated what happens when the transition zone thickness is set to 0 mm in order to investigate whether the contact situation between the bone tissues has an influence on the singularity.However, the singularity occurs even when the transition zone thickness is set to 0 mm: We conclude that the singularity is exclusively caused by the sudden increase in Young's modulus.
The occlusal surface of the crown was loaded with a force of 100 N in axial direction and 20 N in mesial direction.These values are consistent with ranges found in the literature for typical food types (Schindler et al. 1998).However, occlusal load during bruxism or maximum voluntary masticary forces can be much higher (Van Der Bilt 2011).
There are some limitations in this study.Besides the bone of the graded approach, all other materials were treated as isotropic and homogenous.The interfaces between the parts were modelled as bonded, which prevents sliding and friction between the parts.The abutment tightening torque was not modelled in our study.While it can be expected that it can have an effect on the stress magnitude (J€ orn et al. 2014), the effect of the singularity occurs in both modelling variants and our approach works as well when the preload on the abutment screw is considered.The typical microstructure of the implant surface was not taken into account and replaced by a smooth, idealized surface.This simplifies for instance the contacts between bone and implant, or the contacts between implant, implant screw and abutment.Furthermore, the forces were applied directly to the occlusal surface and thus the contact situations between the implantsupported crown and an antagonist are greatly simplified compared to previous studies (Rand et al. 2017;Brune et al. 2019).Those changes were done to reduce the complexity of the model.It is expected that the findings of this study apply to more sophisticated contact situations as well.While the represented modelling approach has no direct clinical implications, the improvement in predicting the stresses in the peri-implant bone more accurately can serve as basis for follow-up studies, such as determining the best implant to use in a certain situation.

Conclusions
The following conclusions can be drawn from this study: 1.The deployed modelling technique can be used to prevent stress singularities and can predict occurring stresses more accurately.
2. A graded bone model leads to a smoother stress distribution for threaded and non-threaded dental implants compared to the conventional modelling approach.3. The graded approach reduces the complexity of the model and the computational time of the simulation compared to the conventional approach

Figure 1 .
Figure 1.The geometry that was used in this study.The units of all shown measures are indicated in millimetres.

Figure 2 .
Figure 2. Forming a new transition zone between the cortical and cancellous bone: The graph on the left shows a step decrease of Young's modulus, which occurs in the conventional approach when modelling two separate bodies.The graph on the right shows the graded course of Young's modulus by including a transition zone, which was acquired by applying a position dependent value for Young's modulus.

Figure 3 .
Figure 3. Flowchart which shows the process of assigning the Young's modulus to the bone.

Figure 4 .
Figure 4.The boundary conditions of the model from two different viewpoints.A force was directly applied to the red surface.The faces in blue depicts frictionless support, the black face depicts fixed support.

Figure 5 .
Figure 5. Distribution of the maximum principal stresses in the peri-implant bone with a mesh size of 0:104 mm: (a) Shows the conventional modelling approach with a non-threaded implant and (b) the graded modelling approach with a non-threaded implant.(c) Shows the conventional modelling approach with a threaded implant and (d) the graded modelling approach with a threaded implant.The white ellipses indicate the locations with the highest maximum principal stresses.The black dashed lines indicate the position of the transition zone.(e) Shows the point of view for the simulation results.

Figure 6 .
Figure 6.Highest maximum principal stresses in the peri-implant bone near the transition for different mesh sizes.The areas where the highest stresses occurred are marked with white ellipses in Figure 5a-d.

Table 1 .
The relevant material properties of the components used in this study.