FEM and BEM Stress Analysis of Mandibular Bone Surrounding a Dental Implant

In the present work the structural behaviour of a mandible with a dental implant, considering a unilateral occlusion, is numerically analysed by means of the Finite Element Method (FEM) and the Boundary Element Method (BEM). The mandible, whose CAD model was obtained by computer tomography scans, is considered as completely edentulous and only modelled in the zone surrounding the implant. The material behaviour of bone is assumed as isotropic linear elastic or, alternatively, as orthotropic linear elastic. With reference to the degree of osteo-integration between the implant and the mandibular bone, a partial osteo-integration is considered; consequently a nonlinear contact analysis is performed, with allowance for friction at the interface between implant and bone. A model of a commercial dental implant is digitised by means of optical 3D scanning process and fully reconstructed in all its geometrical features. Special attention is drawn to the mathematical reconstruction of the CAD model in order to facilitate the meshing process in the BEM environment and reduce the geometrical imperfections generated during the CAD to CAE translation process. The results of FEM and BEM analyses in terms of stress distribution on the mandible are compared in order to benchmark the two methodologies against accuracy and pre-processing efforts.


INTRODUCTION
Endosteal dental implants can cause resorption in the surrounding bone, leading to gradual loosening and ultimately to a complete loss of the implant; in particular a direct correlation was found between overstressed regions and bone resorption [1].Stress distribution in the bone strongly depends on the implant shape so it becomes of uttermost importance to test different commercial implants in order to devise the configuration providing the lowest possible stress concentration in the bone, thereby reducing the resorption risk.
In [2] authors studied the stress distribution on endosseous dental implant and surrounding bone, by using both Finite Element Method (FEM) and Boundary Element Method (BEM) in a three dimensional modelling approach.A FEM-based study in [3] has pointed out the impact of the articular disc stiffness and of the temporo-mandibular joint friction coefficient on the mandible stress peaks and occlusal forces.A later comparison also with BEM was then carried out in [4].In [5,6] also the occlusal stress transmitted to the inferior alveolar nerve was analysed.In [7] the biomechanical behaviour, in terms of stress concentration and distribution, of different commercial dental implants with different thread profiles was studied with FEM, analysing conditions of both perfect and partial osteointegration.Similarly, by using FEM in [8] a statistical approach was described to evaluate the geometrical parameters of the implant that significantly affect the induced stresses and damage in the bone.
All these studies require the modelling of the dental implant and of the surrounding bone topology.In some cases the analysis of simplified geometrical models is sufficient for a preliminary analysis, but for a closer investigation, such that adopted in the present study, detailed geometries of the real models are necessary.The virtual model of the dental implant can be generated by using high density optical scanners providing a dense point data set, which is used to create a tessellated/polygonal model that can be converted into a CAD model by fitting patch surfaces following the common reverse engineering (RE) procedures.
On the other hand, the process of virtual modelling of anatomical parts is slightly different: the data set usually comes from 2D medical imaging like CT scan; it is then converted into a tessellated 3D model by volumetric segmentation algorithms.The polygonal model can be edited with RE techniques or, taking into account the specific topology of the part, CAD-like modelling functions (extrusion, sweep, loft, blend) can be used to create the final 3D CAD model [11].
The present work starts from an accurate modelling of a commercial dental implant as well as the human edentulous mandible and is intended for numerically analysing the stress gradient induced by the loaded implant, using both FEM and BEM approaches under different load conditions.The accurate 3D CAD models of both mandible and implant is performed to tackle the complexity of models closer to real geometries and point out pros and cons of the numerical analyses carried out with Comsol Multiphysics (FEM) and BEASY (BEM) commercial codes.Specific attention, as detailed in the next section, is drawn to BEM model generation, due to some restrictions of the adopted commercial code (BEASY) when involving the CAD-CAE data exchange and meshing procedure.
Numerical simulations involve axial and inclined loads, not-perfect osteo-integration, and isotropic and transversely isotropic material models.Also the effects of the friction have been studied.
The cross validation between the two methodologies is mainly based on the comparison of the calculated pressure distributions at the bone implant interface.This is a critical parameter when considering the problem of the micromotion and fretting damages at the dental implant/bone interface, where an accurate evaluation of contact pressure is of the uttermost importance for the correct assessment of the initial success of osteo-integration and the life time of dental implant [9,10].

CAD MODELLING
CAD model of the mandible, related to an edentulous patient, is generated starting from CT scan images (Fig. 1), in order to build up a 3D tessellated model.Then this model is edited in SolidWorks 2012 by using the ScanTo3D add-on [11]: some cross-section curves, taken with different orientations according to the topology of the mandible, are used to model the implant area by using a CAD loft surface (Fig. 2).The cortical and spongy bones are reconstructed with the same approach as shown in Fig. (3).This procedure, also successfully applied in [12], offers the advantage, over other reconstruction methods (in particular those based on patch fitting), of providing few smooth surface patches (just one in our case), amenable to easy handling and meshing in a FEM/BEM environment.The level of accuracy depends on the number of cross-section curves.
In parallel, a model of a commercial implant (DentSply FRIADENT ® ) is digitised by means of a high density structured-light 3D scanner (Comet 5 by Steinbichler Optotechnik GmbH), and fully reconstructed in all its geometrical features (Fig. 4) using Rhino3D v.4, which offers several tools to edit meshes and generate surface models.The implant model is then imported as solid geometry into SolidWorks 2012 and assembled with the mandible model (Fig. 5).The choice of the latter CAD system is based on its parametric features that allow to easily control the implant position.Moreover, it is well integrated into Rhino3D and into the adopted FEM code (Comsol Multiphysics), as described below.Fig. ( 5) also shows the subdivision of the whole model into three portions to facilitate the sub-modelling analyses made in BEASY.A Boolean operation is adopted in order to generate proper interfacial boundaries between mandible and dental implant.
Comsol Multiphysics 4.4 has a direct link to SolidWorks CAD system so that CAD models can be directly read without losing geometrical and topological information.The so-generated model is meshed in Comsol tightly controlling the mesh in the surrounds of the implant threads.
IGES neutral format, instead, is used to import CAD models into BEASY package.A successful translation into BEASY is obtained only if appropriate tolerances and topology settings are chosen: CAD surface patches are firstly converted into separate trimmed surfaces (144 IGES entities), and then tolerance settings are chosen as described in [13].Nevertheless, this is not sufficient to get an acceptable model mesh in BEASY: as shown in Fig. (6), patches with more than 4 edges are roughly meshed and this may affect the accuracy of the analysis.A regular mesh is only possible with 3-or 4-edge patches.Due to this limitation the CAD model is firstly checked and processed in Rhino3D (for its wide set of specific commands) by manually editing edges and splitting/merging patches so as to generate regular patch shapes.The so created final model is then exported via IGES into BEASY code.

MATERIALS AND METHODS
Material properties greatly influence the stress-strain distribution in a bone-implant structure [5,14].If, on the one hand, for the implant a linear elastic and isotropic behaviour can be assumed (dental implants are usually made of Titanium alloys), on the other hand bone tissue is a more sophisticated structure, whose mechanical properties depend on several factors: age, gender, pathology, anatomical site, liquid content and so on [15].With respect to the maxillarymandible site, several contributions pointed out the bone mechanical properties [16].Many of these contributions assumed both cortical and spongy bones as linear and isotropic.As demonstrated also in [17], this assumption corresponds to the approximation of a type II bone density.Other contributions, instead, considered a more realistic anisotropic law [2][3][4][5][6], also trying to capture the change of bone density in all space directions.A review of the available literature on computational modelling in the area of bone biomechanics, fracture and healing is available in [18].Since the aim of the present research is to compare numerical results coming from FEM and BEM numerical approaches, pointing out advantages and drawbacks in terms of stress-strain accuracy, both cases of isotropic and transversely isotropic (Fig. 7) behaviour for both cortical and spongy bones are considered.
Tables 1 and 2 report the adopted isotropic and transversely isotropic material mechanical properties, respectively.
Mastication involves repeated cyclic forces that are transferred to the bone tissue by the implant.Literature shows that mastication loads can significantly vary from one area of the mouth to another.As reported in [15], average forces of nearly 800 N for male young adults and 600 N for female young adults have been recorded in the molar region.In the premolar region, occlusal forces range from 200 to 600 N. Forces lower than 200 N have been measured in the incisal region.Such a variation may be related to many factors, such as muscle size, bone shape, gender, age, degree of edentulism.
In the present paper, we work with an endosseous implant located in the premolar-incisal regions.The occlusal load is assumed equal to 350 N, statically applied on the upper surface of the implant and acting along its longitudinal axis or with a given inclination (α 1 =16°), configuring in such a way an added lateral mastication force (Fig. 7).Only half of the original reconstructed mandible is considered for numerical analyses, clamped on the two cutting surfaces as shown in Fig. (8).In all the analyses contacts are introduced at the interface between implant and cortical bone (Fig. 9), with and without friction as detailed in the following.Where friction is applied it is assumed µ=0.42.
Table 3 summarises the five alternatives hypotheses adopted in the numerical FEM and BEM analyses.
The compenetration between implant and cortical bone at the ring "planar supporting" area (Fig. 10) can be allowed, using internal spring of negligible stiffness (Hyp. 2 to 5), or not allowed using contact pairs (Hyp.1).The former configures a more realistic condition in which the occlusal load is mainly absorbed at the implant thread level rather than in a very restricted area: the ring "planar supporting" area would absorb nearly 60% of the applied load if compenetration is not allowed.Moreover, the sliding between implant and bone is modelled with or without friction (Hyp. 2 to 5): where applied, the friction coefficient µ is set to 0.42 [2].

FE MODEL
As previously mentioned, the CAD model shown in Fig. ( 8) is exported from SolidWorks to Comsol by using the LiveLink connection by Comsol.In the FEM environment the tetrahedral mesh is generated (Fig. 11) with a high refinement with linear elements in the middle portion of the model around the implant area and, in particular, along the threaded region between the dental implant and the cortical bone, in order to obtain an accurate analysis in the contact zone.The final model consists of about 730'000 degrees of freedom (dofs).Numerical simulations are performed by considering a not-perfect osteo-integration.This is modelled by defining contact pairs at the implant-cortical interface.We assume that the implant has a lack of integration just in correspondence of those threads interacting with the cortical bone, where the stress reach the highest values.Thus, the implant-spongy interface is modelled with a node-to-node identity so having a congruent mesh within Comsol.The choice of linear elements is done in order to reduce the number of dofs and to speed up the solution convergence.

BE MODEL
BE model (Fig. 12) is based on a mixed linear-quadratic mesh with 4'300 elements: the higher polynomial order is provided in correspondence of the interface between the cortical bone and the implant, where highest stress gradients are expected.In the BEM environment we adopt the approach based on the extraction of a sub-model involving the volume surrounding the implant, i.e. the middle portion of the whole model (Fig. 12c).As for the FEM model, in case of not perfect osteointegration between the implant and the cortical bone, the analysis becomes nonlinear with gap elements at the interface; for the remaining part of the implant, in contact with the spongy bone, a perfect osteo-integration is supposed and consequently the continuity is modelled for tractions and displacements.
When considering a pure axial load, a mesh made of only linear elements (also in the critical area of implant/cortical bone interface) is sufficient to get accurate results, whereas, with an inclined load the use of quadratic elements in the critical areas is recommended.

COMPARISON BETWEEN FEM-BEM RESULTS
For each hypothesis listed in Table 3 the results of FEM-BEM analyses are reported and discussed in the following.
FEM and BEM models are based on about 730'000 and 39'000 dofs, respectively.Fig. (13) shows the resultant displacements provided by FEM and BEM simulations under the hypothesis "hyp.1".
Normal tractions reach very high values in correspondence of the edge of the contact area between implant and bone (the ring "planar supporting" area), also because of a zero radius fillet with consequent strong notch effects (Fig. 14); as a result the only first thread turns out to be loaded by non negligible pressure values.
With reference to the geometric points highlighted in Fig. (15), a quantitative comparison between FEM and BEM results is provided in Fig. ( 16), a satisfactory consistency with reference to both displacements and pressures.In particular, the maximum difference between FEM and BEM displacements is lower than 2% (Fig. 17).
Fig. (18) shows FEM and BEM contour plots of normal tractions at the implant-cortical bone interface, with axial load, bone isotropic properties and no friction according to the hypothesis "hyp.2".Differently from the previous case, the load is mainly absorbed by the whole cortical bone thread, thus providing a more realistic condition [20] in comparison to "hyp.1".The level of consistency between FEM and BEM results is satisfactory.friction: a slight pressure decrease is observed in comparison with the case at µ=0.Stress concentration areas are localised around implant neck at the cortical bone interface.A comparable behaviour is reported in the work performed by Citarella et al. [2], in terms of stress field at cortical bone-implant interface.Fig. (20) shows FEM and BEM contour plots of normal tractions at the implant-cortical bone interface under the hypothesis "hyp.3", with inclined load, isotropic material properties and no friction.In this case quadratic elements are used in the BEM simulations in order to reach a satisfactory agreement with FEM results.In both cases it is possible to remark the load concentration on one side of the implant, even with a relatively small lateral load component.

CONCLUSION
This paper shows a cross comparison between two numerical methodologies (FEM and BEM) as implemented in two different commercial codes applied to a biomechanical problem.
Different conditions are adopted to handle dental implant issues such as the isotropic and transversely isotropic material behaviour for cortical and spongy bones as well as the presence or absence of friction on the interface between implant and cortical bone.Moreover, the modelled complex geometries are very close to the real ones as accurate reconstruction techniques are used.
Numerical results demonstrate an overall good agreement between the two methodologies even though they are related to the specific software adopted in the simulations: Comsol Multiphysics for FEM and BEASY for BEM.With respect to these codes some remarks can be outlined from this study.
The CAD-CAE integration level is very high in Comsol environment as a direct link with common CAD software is available and well tested for several native and neutral formats of data exchange.Thus, complex geometric models can be successfully imported and, in case of data exchange   problems, healing tools are available to repair the model.The internal automatic mesher performs satisfactorily with little user intervention.Nevertheless, a very dense mesh is required to get good results, especially in case of nonlinear contact analysis without or with friction.In the latter case, in the previous versions of the Comsol software it was necessary appropriately setting the internal parameters for the penalty factor, while with the used version 4.4 it is just enough to adopt the suggested preset parameters and the solution converges.
On the other hand, using BEASY code, it is difficult to handle complex geometries and particular attention must be devoted to the right settings of the parameters controlling the translation of native CAD geometries.Moreover, very few repairing tools are available to heal topology issues, so that a trial and error approach is often necessary.To obtain a good quality mesh the model topology must be controlled at a CAD level: an accurate tailored modelling task is necessary to generate in the CAD environment the single patches with only 3-or 4-edges.For this task Rhino3D is used: it provides several specific and powerful editing tools, not available in many other CAD environments.This is a very timeconsuming task.Nevertheless, when facing non-linear simulations, involving contact pairs, BEM approach gives results comparable with those obtained with FEM on the boundary pairs, without drastically reducing the mean mesh size.On the contrary, the FEM solution is quite sensitive to the local mesh size around the contact pairs: to improve the accuracy, a very fine local mesh is required.
The comparison between the two models is judged acceptable, although it is the authors' opinion that for a complete validation further analyses should be carried on.

Fig. ( 3
Fig. (3).CAD model of the mandible with highlight of cortical and spongy bone.

Fig. ( 6 ).
Fig. (6).Detail of the distorted mesh originally generated in BEASY due to the multi-edge patches coming from the CAD model.

Fig. ( 9 ).
Fig. (9).Highlight of the interface area between implant and cortical bone where contacts are added.

Fig. ( 10
Fig. (10).The occlusal load is mainly absorbed by the planar supporting provided by a very limited ring external cortical bone surface.

Fig. ( 19
Fig. (19) shows FEM and BEM contour plot of normal tractions at the implant-cortical bone interface in case of

Fig
Fig. (22) shows FEM and BEM contour plots of normal tractions at the implant-cortical bone interface under the

Fig. ( 25
Fig. (25) shows the FEM and BEM contour plots of normal tractions at the implant-cortical bone interface, in case of transversely isotropic properties and inclined load according to the hypothesis "hyp.5".Again, quadratic elements are necessary in the BEM simulations in order to get a satisfactory consistency with FEM results.The comparisons of the contour plots of normal tractions and maximum principal stresses at the cortical bone in case of friction are depicted in Figs.(26, 27), respectively, displaying, as in the previous cases, a sound agreement.