Analytical Study of Stress Distributions around Screws in Flat Mandibular Bone under In-Plane Loading

A known complication for mechanically loaded bone implants is the instability due to screw loosening, resulting in infection and the non-union of fractures. To investigate and eventually prevent such bone degradation, it is useful to know the stress state in the bone around the screw. Considering only in-plane loadings and simplifying the mandibular bone into an orthotropic laminated plate, the analysis was reduced to a two-dimensional pin-loaded plate problem. An analytic model, based on the complex stress analysis, was introduced to the bone biomechanics field to obtain the stress distributions around the screw hole in the bone. The dimensionless normalized stresses were found to be relatively insensitive to the locations of the screw hole over the mandible. Parametric analyses were carried out regarding the friction coefficient and load direction. It was found that the load direction had a negligible influence. On the contrary, the friction coefficient had a significant effect on the stress distributions. Whether the screw was well bonded or not thus played an important role. The proposed analytic model could potentially be used to study bone failure together with stress-based failure criteria.


Introduction
Fractures of human bones, e.g., the mandible, require adequate fixation in a correct occlusal position in order to allow the absorption of masticatory forces. Adopting plates fixed with screws is the standard procedure for the treatment of mandibular fractures currently, mainly because of its advantages of minimal tissue injury and high fixation stability [1]. However, bone loss is observed in patients after fracture treatment, especially when load bearing plates are used or bone grafting is involved. Resorption of bone surrounding the screws will consequently result in screw loosening, which in turn may cause instability, bone infection, further loosening, and occasionally osteomyelitis [2][3][4]. This leads to the need for surgical revision, removing the fixation material and new adequate fixation with plates and screws [5]. The failure is most often due to the incorrect analyses of surgeons about the biomechanical stability requirement in a specific fracture location and pattern and/or the wrong choice of material and an inadequate anatomical location for osteosynthesis not providing mechanical stability. From a mechanical perspective, the stress shielding and the high peak stresses at the interface are believed to be determinant for bone loss [6][7][8][9][10]. Thus, better understanding of stress distributions around screw holes is of great importance for preventing implant failures.
Several studies have been performed in the past few decades to investigate the biomechanics of bones around artificial implants. Most of the research is based on finite element simulations [11][12][13]. They provided detailed and precise stress distributions in the bone at the cost of a considerably large amount of time and computations. An analytic model to capture trends and parametric influences would be useful in the general understanding and model to capture trends and parametric influences would be useful in the general understanding and for the design approach. The human mandible consists of two cortical layers and an intermediate layer of cancellous bone. For simplification, mandibular bones were idealized and treated as thin laminated plates in this work. This is acceptable since, in this work, only the in-plane loads were considered. The biomechanical problem could thus be converted into the well-solved pin-loaded plate problem, which has been researched substantially [14][15][16][17][18]. Among the theoretical investigations, most are analytical studies based on the complex stress function method [16]. De Jong adopted the general solution of the complex stress function to analyze stress distributions around a hole in a linear elastic plate loaded by a rigid pin [14]. The study was limited to frictionless and tight-fit connection. Hyer and Klang combined the complex stress function method with a collocation method in a general stress analysis, in which the deformation of the pin, clearance, and friction between the pin and the plate were all accounted for [17]. Zhang and Ueng obtained a compact analytic solution for the stress functions, which satisfies the major displacements on the hole loaded by a rigid, perfectly fitting pin [18].
The aim of the current work was to revive and introduce the complex stress analysis method into the bone biomechanics field for predicting the stress distribution around a screw-loaded hole in human mandibular bone. Parametric studies of the stress distribution in the bone regarding the relevant design parameters were easily conducted by utilizing an analytic model derived with the complex stress analysis method. Since the stresses determined were averaged with respect to the bone thickness, the stress distribution was a coarse approximation. Stress approximations in the cortical and cancellous layers can, however, be obtained from the elastic properties of these layers. The idea of the compact analytic model is to provide rough, but timely guidelines in the design process of the screws and the locations of fixation in the bones. For more accurate studies for cases where an in-plane pin loading and plate analogy cannot be assumed, an analytical approach is most likely not possible. In that case, one must resort to a full three-dimensional finite element analysis. Notably, such analytic models could easily be expanded and applied to other types of flat bone, such as the skull, in human bodies.

Stress Analysis
Prior to the analytical analysis, the following simplifications were made and are listed regarding the screw and mandibular bone, schematically illustrated in Figure 1. (1). The screws are supposed to be symmetrically loaded herein. As reported, throughthickness bicortical screw placement is often used for the rigid fixation of mandibular fractures [19]. Thus, an in-plane loaded rigid pin in a through hole was considered reasonable here. The effects of transverse pull-out loads and contributions from large moments deserve special attention and necessitate a full three-dimensional finite element study. (1). The screws are supposed to be symmetrically loaded herein. As reported, throughthickness bicortical screw placement is often used for the rigid fixation of mandibular fractures [19]. Thus, an in-plane loaded rigid pin in a through hole was considered reasonable here. The effects of transverse pull-out loads and contributions from large moments deserve special attention and necessitate a full three-dimensional finite element study. (2). The screw-load direction is arbitrary with respect to the principal material direction, as the loading direction has an influence on the strength of pin-loaded anisotropic plates [20,21].
(3). The screw was assumed to be infinitely rigid. According to Hyer et al. [22], the stiffness of the pin has a minor effect on stress distribution around the hole. (4). Threads on the screw and screw hole were neglected. (5). The location of the screw hole in the bone was assumed to be far away from the edge of the bone. The bone was considered to be of uniform thickness and flat in the vicinity of the screw. (6). The contact between the screw and the bone was considered, and the friction coefficient was assumed to be constant over the contact region.

Elastic Properties of Human Bones
The elastic properties of bones are of fundamental interest in functional morphology, which describe how anatomical structures resist or deform by external loads [23]. In this context, the mandible was assumed to be homogenous and elastically orthotropic [24]. The elastic properties of mandibular bone have been experimentally determined [25][26][27], and the values for different regions of the mandibular cortical bone are presented in Table 1. The values were taken from [28], which was based on the experimental study in [29]. The material directions were in accordance with local Cartesian coordinate systems, in which Direction 3 is normal to the bone plane. Table 1. Orthotropic elastic properties of cortical bone [28]. The cancellous bone in the mandible was assumed to be elastically isotropic in this analysis, with an elastic modulus of 1.37 GPa and a Poisson ratio of 0.3 [30]. According to its anatomic structure, the mandible is represented by a laminated plate comprising three layers [31], as shown in Figure 2. The effective elastic properties were calculated through classical laminate theory (e.g., [32]):

Positions in Mandible
where H is the thickness of the bone; A ij are the components of the stiffness matrix, defined as in which C ij k are the stiffness matrix components of the kth layer and t k is the thickness of the kth layer.

Stress Analysis Based on Complex Stress Functions
The stress distribution in an orthotropic plate with elliptic or circular openings has been extensively explored in the literature. Analytical solutions of the stress distribution in orthotropic homogeneous plates with circular openings that are filled by an elastic material with different elastic properties have been derived by Lekhnitskii [33]. As has been stipulated, when the opening (far from the edge) is relatively small in comparison with the plate size, the plate could be simplified as infinite, and the effect of the external edge is neglected. The general solutions of the complex stress functions and stress components in this analysis followed those in [18,33]. Only the formulas are reproduced and explained here with consistent notation, given by Equations (3)-(9).

General Solution of Stress Functions
A Cartesian coordinate system -was set up, such that the displacement of the pin was in the direction of the positive -axis, as shown in Figure 3. The radii of the screw and screw hole are identical and are denoted by . The resultant force from the pin to the plate is known. The principal material directions of the plate are such that the intersection angle between the principal material direction and coordinate axis is . The Young's moduli along the principal material directions and are and , respectively. The shear modulus is , and the Poisson ratio characterizing the lateral shrinkage in direction caused by extension in direction is . Plane stresses were derived by solving the governing biharmonic differential equation in the complex domain [33]: where are complex stress functions and the first-order derivatives of functions are with respect to = + . The variables are the complex roots of the characteristic equation in the coordinate system -, while for principal coordinate system -. According to Lekhnitskii [33], if the resultant force components are known and displacements * , * at the plate hole are given by

Stress Analysis Based on Complex Stress Functions
The stress distribution in an orthotropic plate with elliptic or circular openings has been extensively explored in the literature. Analytical solutions of the stress distribution in orthotropic homogeneous plates with circular openings that are filled by an elastic material with different elastic properties have been derived by Lekhnitskii [33]. As has been stipulated, when the opening (far from the edge) is relatively small in comparison with the plate size, the plate could be simplified as infinite, and the effect of the external edge is neglected. The general solutions of the complex stress functions and stress components in this analysis followed those in [18,33]. Only the formulas are reproduced and explained here with consistent notation, given by Equations (3)-(9).

General Solution of Stress Functions
A Cartesian coordinate system x-y was set up, such that the displacement of the pin was in the direction of the positive x-axis, as shown in Figure 3. The radii of the screw and screw hole are identical and are denoted by r. The resultant force P from the pin to the plate is known. The principal material directions of the plate are such that the intersection angle between the principal material direction x 0 and coordinate axis x is ϕ. The Young's moduli along the principal material directions x 0 and y 0 are E x 0 and E y 0 , respectively. The shear modulus is G x 0 y 0 , and the Poisson ratio characterizing the lateral shrinkage in direction y 0 caused by extension in direction x 0 is v x 0 y 0 . Plane stresses were derived by solving the governing biharmonic differential equation in the complex domain [33]: where φ k are complex stress functions and the first-order derivatives of functions φ k are with respect to z k = x + µ k y. The variables µ k are the complex roots of the characteristic equation in the coordinate system x-y, while µ k 0 for principal coordinate system x 0 -y 0 . According to Lekhnitskii [33], if the resultant force components are known and displacements u * , v * at the plate hole are given by The general solution of stress functions and could be expressed as Here, the arbitrary constants were neglected since only derivatives of ( ) and ( ) were of interest for stresses; is the imaginary unit, and is the angle with the positive -direction, as shown in Figure 3; the bar sign "¯" denotes the complex conjugate; and are functions of complex variables and , respectively; , , , , , , and are constants determined by the material properties. All aforementioned parameters are defined in [18,33].

Boundary Conditions
The unique solutions of the stress components for the pin-loaded plate problem depend plainly on the unknown parameters and , which should fulfill specific boundary conditions of the problem. In the present work, we made use of the boundary conditions proposed previously by other researchers. According to Zhang and Ueng, the major displacement boundary conditions of this problem may be formulated as [18] where and are the displacements of the points on the plate hole circumference in the -and -axis. , / , and / are the displacements of points , , and , shown in Figure 3, respectively. Furthermore, denoting the friction coefficient as , certain stress boundary conditions on the plate-hole circumference are also introduced [18]: The general solution of stress functions φ 1 and φ 2 could be expressed as Here, the arbitrary constants were neglected since only derivatives of φ 1 (z 1 ) and φ 2 (z 2 ) were of interest for stresses; i is the imaginary unit, and θ is the angle with the positive x-direction, as shown in Figure 3; the bar sign " " denotes the complex conjugate; ξ 1 and ξ 2 are functions of complex variables z 1 and z 2 , respectively; A, B, D, p 1 , p 2 , q 1 , and q 2 are constants determined by the material properties. All aforementioned parameters are defined in [18,33].

Boundary Conditions
The unique solutions of the stress components for the pin-loaded plate problem depend plainly on the unknown parameters α m and β m , which should fulfill specific boundary conditions of the problem. In the present work, we made use of the boundary conditions proposed previously by other researchers. According to Zhang and Ueng, the major displacement boundary conditions of this problem may be formulated as [18] where u and v are the displacements of the points on the plate hole circumference in the xand y-axis. u 0 , u 0 /c 1 , and u 0 /c 2 are the displacements of points M, M 1 , and M 2 , shown in Figure 3, respectively. Furthermore, denoting the friction coefficient as f , certain stress boundary conditions on the plate-hole circumference are also introduced [18]: Bioengineering 2023, 10, 786 6 of 15

Exact Solution of Stresses
To determine the unknowns α m and β m , the expressions for displacements u and v are represented in the following form [18]: where U i and V i are unknown constants. Substituting the expressions for u and v into Equation (6), the expressions of U i and V i in terms of u 0 , c 1 , and c 2 are obtained. Comparing the corresponding terms from Equations (4) and (8), the unknowns α m and β m can then be found. Thus, the stress components σ xx , σ yy , and τ xy can be obtained from Equations (3) and (5). Applying the stress transformation between the polar and Cartesian coordinate systems, compact solutions for stress components in the polar coordinate system can be expressed as [18] where S, C σ 1 , S σ 1 , C τ 1 , S τ 1 , etc. are expressed in detail in Appendix A; P x = Pcos θ 1 is the x-component of screw-load P, in which θ 1 is the angle between the screw-load direction and the positive x-direction. Unknowns u 0 , c 1 , and c 2 are determined by substituting the stress components into the stress boundary condition, e.g., Equation (7).
For a finite-width plate, the stress concentration factor can be calculated by the use of the finite-width correction factor [34,35], which is a scale factor defined as the stress concentration factor of an infinite-to a finite-width plate. Thus, the stress concentration of the finite-width plate can be calculated with the specification of the finite-width correction factor and the stress concentration of the infinite plate with a similar opening and material. Since only trends and relative contributions from the inherent parameters to stress concentration were in focus in this work, the finite-width corrections were not used here.

Results and Discussions
The stress distributions around the screw hole depend on the elastic properties of cortical and cancellous bone, the friction coefficient, and the loading direction. To illustrate the influence of different parameters, we determined the numerical values of the stresses for the different regions in the mandibular bones based on the analytical model. The friction coefficient between the bones and the porous titanium alloy was around 0.6 [36]. Thus, the friction coefficient f in the present work varied from 0, 0.2, 0.4, 0.6 to 1, where the frictionless case would correspond to an unbonded screw. With the increase of f , the screw is considered to become increasingly bonded to the bone material. For each friction coefficient value, the profiles of the radial normal stress σ rr and the circumferential shear stress τ rθ were calculated for different cross-angles ϕ, i.e., 0 • , 15 • , 30 • , 45 • , 60 • , 75 • , and 90 • , between the loading direction and the direction of the main material axis.
The stress distributions around the screw holes in different regions on the mandible were obtained for different combinations of the friction coefficient and the cross-angle. The analytic model predicts the average stresses for the homogenized bone laminate. The dimensionless relative thickness t r of the cortical bone layer is defined as where t cort is the thickness of the cortical bone layer and t can is the thickness of the cancellous bone layer. The relative thickness varies for different regions of the mandible bone. For the sake of simplification, the relative thickness for each of the superficial cortical bone layers was assumed as 0.2, 0.4, 0.5, 0.6, 0.8, and 1.0 for the symphysis, body, angle, ramus, condyle, and coronoid bones, respectively. For a relative thickness of 1.0, the bone is comprised only of the stiff cortical bone, which can be assumed to be the case for the coronoid part of the mandible. The effective elastic property of the mandible bone is orthotropic due to the symmetry and the stiffness values calculated from classical laminate theory, as listed in Table 2. Table 2. Effective orthotropic elastic properties of different regions of the mandible. Mandible

Average Stresses over the Cross-Section of the Laminated Bone
It was found that the differences between the stress distributions in different regions were not significant for the various values of f and ϕ. In order to limit the number of combinations of f and ϕ, only the distributions of the normalized stresses σ rr /(P x /2r) for f = 0 and ϕ = 0 • and τ rθ /P x /2r for f = 1 and ϕ = 90 • , which are two extreme cases, are illustrated in Figure 4. Shear stress τ rθ does not exist when the friction coefficient f = 0 and, thus, was not demonstrated. As depicted, different line colors imply the stresses in different regions, which are indicated by red roman numerals. It is clear from Figure 4 that, for the two extreme cases, the distributions of both the normalized radial and shear stresses were not distinguishable for different regions of the mandible, except for the coronoid region. Notably, the computational method was deemed to be only valid under the assumption that the screw was located far from the bone edge, which means that the stress distributions in the condyle and coronoid regions were less reliable. When there was no friction, the maximum radial stress occurred at the location along the loading direction (Figure 4a), as expected. With the presence of friction, the shear stress increased (Figure 4c), and the distribution of the radial stress varied with its maximum location moving away from the loading direction, as can be seen in Figure 4b.
Since the differences of the average stress distributions over the mandible were fairly small, the angle region in the mandible, where the largest percentage of mandibular fractures occur [37], was chosen to study the influence of the friction coefficient f and the cross-angle ϕ in more detail. In Figure 5, the stress distributions are demonstrated for the angle region (denoted as III), where different colors are representations for different values of the friction coefficient f and different brightness intensities stand for different values of the cross-angle ϕ. As can be observed, the stress distribution curves could be divided approximately into five groups with different levels of friction, for both the normalized radial and shear stresses. With the increase of the friction coefficient, the radial stress level around the screw hole gradually became lower [18,22]. Simultaneously, the distribution for the radial stress in radial direction σ rr changed with the high-stress region spreading out bilaterally from θ = 0 • . Comparing the curves in each group, the effect of the loading direction, ϕ, can be observed only under closer scrutiny, and the stress distributions remained relatively constant. The variation of the friction coefficient, f , obviously had a larger impact on the stress distributions [38]. However, readers should bear in mind that the stress distributions derived in this work are only valid for screw-bone contact with large values of the friction coefficient, i.e., the bond between the screw and the bone should be tight.  Since the differences of the average stress distributions over the mandible were fairly small, the angle region in the mandible, where the largest percentage of mandibular fractures occur [37], was chosen to study the influence of the friction coefficient f and the crossangle φ in more detail. In Figure 5, the stress distributions are demonstrated for the angle region (denoted as III), where different colors are representations for different values of the friction coefficient f and different brightness intensities stand for different values of the cross-angle . As can be observed, the stress distribution curves could be divided approximately into five groups with different levels of friction, for both the normalized radial and shear stresses. With the increase of the friction coefficient, the radial stress level around the screw hole gradually became lower [18,22]. Simultaneously, the distribution for the radial stress in radial direction changed with the high-stress region spreading out bilaterally from = 0°. Comparing the curves in each group, the effect of the loading direction, , can be observed only under closer scrutiny, and the stress distributions remained relatively constant. The variation of the friction coefficient, , obviously had a larger impact on the stress distributions [38]. However, readers should bear in mind that The critical point for possible damage formation, e.g., cracking, is at the location of the maximum stress levels. To demonstrate the influence of the friction and loading direction, the curves of the maximum stresses for the angle region with respect to f and ϕ are plotted in Figure 6. With the increase of the friction coefficient, the maximum radial stress decreased while the maximum shear stress increased steadily. On the contrary, the maximum stresses varied only within a small range when the loading direction varied. Not only for the stress distributions in Figure 5, but also for the maximum stresses, the friction played a more important role than the direction of the load.

Potential Applications of the Analytical Model
It is important to keep in mind that the current analytical method could be used together with other analysis methods to implement a complete analysis of failure prediction in the case of in-plane loading. A schematic procedure is proposed and shown in Figure 7. The present study was limited to the highlighted orange box in Figure 7, i.e., the two-dimensional stress analysis under in-plane loading. The schematic shows how this analysis could fit into a larger framework. As can first be seen, the kinetics information of the mandible and the implant could be obtained through a macro-level analysis using commercial software, e.g., in [39]. The in-plane loads applied to the local screws were then used as the input for the 2D analytical stress analysis. The overall stress distributions around the screw hole were determined by superposing the contributions from the in-plane loads and other loads, e.g., moments and transverse loads from the muscular forces, which most likely should be treated numerically. Eventually, the stress distributions will be further evaluated, using currently available bone fracture criteria such as the Rankine criterion, the Mohr-Coulomb criterion, and the Tsai-Wu criterion [40], to predict bone failure.
Bioengineering 2023, 10, x FOR PEER REVIEW 9 of 15 the stress distributions derived in this work are only valid for screw-bone contact with large values of the friction coefficient, i.e., the bond between the screw and the bone should be tight. The critical point for possible damage formation, e.g., cracking, is at the location of the maximum stress levels. To demonstrate the influence of the friction and loading direction, the curves of the maximum stresses for the angle region with respect to f and are plotted in Figure 6. With the increase of the friction coefficient, the maximum radial stress decreased while the maximum shear stress increased steadily. On the contrary, the maximum stresses varied only within a small range when the loading direction varied. Not only for the stress distributions in Figure 5, but also for the maximum stresses, the friction played a more important role than the direction of the load.

Potential Applications of the Analytical Model
It is important to keep in mind that the current analytical method could be used together with other analysis methods to implement a complete analysis of failure prediction in the case of in-plane loading. A schematic procedure is proposed and shown in Figure  7. The present study was limited to the highlighted orange box in Figure 7, i.e., the twodimensional stress analysis under in-plane loading. The schematic shows how this analysis could fit into a larger framework. As can first be seen, the kinetics information of the mandible and the implant could be obtained through a macro-level analysis using commercial software, e.g., in [39]. The in-plane loads applied to the local screws were then used as the input for the 2D analytical stress analysis. The overall stress distributions around the screw hole were determined by superposing the contributions from the inplane loads and other loads, e.g., moments and transverse loads from the muscular forces, According to different fracture criteria, it is preferable that the stress states in each bone layer are known, since the failure conditions will be different for different bone structures. The average stresses over the thickness of the laminated bone were calculated in Section 3.1. The local stresses in each layer, e.g., the cortical bone layer, are determined and analyzed hereafter. With only in-plane deformation considered, the average strain vector ε lam of the bone laminate analogy (cf. Figure 2) is expressed as ε lam = HA −1 σ lam (11) where H is the thickness of the bone, A −1 is the inverse of stiffness matrix A, and σ lam is the average stress vector of the laminated bone in a Cartesian coordinate system. Then, the constitutive equation of the cortical bone was used to calculate the local stresses in the cortical bone layer from the average strains: (12) in which D cort is the stiffness matrix of the cortical bone material under the plane stress state in the x-y coordinate system.
Bioengineering 2023, 10, x FOR PEER REVIEW 11 of 15 which most likely should be treated numerically. Eventually, the stress distributions will be further evaluated, using currently available bone fracture criteria such as the Rankine criterion, the Mohr-Coulomb criterion, and the Tsai-Wu criterion [40], to predict bone failure. According to different fracture criteria, it is preferable that the stress states in each bone layer are known, since the failure conditions will be different for different bone structures. The average stresses over the thickness of the laminated bone were calculated in Section 3.1. The local stresses in each layer, e.g., the cortical bone layer, are determined and analyzed hereafter. With only in-plane deformation considered, the average strain vector of the bone laminate analogy (cf. Figure 2) is expressed as where is the thickness of the bone, is the inverse of stiffness matrix , and is the average stress vector of the laminated bone in a Cartesian coordinate system. Then, the constitutive equation of the cortical bone was used to calculate the local stresses in the cortical bone layer from the average strains: The stress profiles around the screw hole in the cortical bone layer are demonstrated in Figure 8. As can be seen in the diagram, the stress distributions in the cortical bone layer were similar to those for the average stresses over the entire cross-section. Nevertheless, the differences between the stresses in different regions of the mandible were more distinguishable. Applying the same in-plane force to the screw, the stress levels in the cortical bone of various regions were different depending on the relative composition of the cortical and cancellous bone. The largest maximum stress level was observed in the symphysis region and decreased with the screw location moving towards the coronoid region. If the same stress-based failure criterion was applied for the cortical bone, failure was more likely towards the symphysis region since the stresses were higher there. the differences between the stresses in different regions of the mandible were more distinguishable. Applying the same in-plane force to the screw, the stress levels in the cortical bone of various regions were different depending on the relative composition of the cortical and cancellous bone. The largest maximum stress level was observed in the symphysis region and decreased with the screw location moving towards the coronoid region. If the same stress-based failure criterion was applied for the cortical bone, failure was more likely towards the symphysis region since the stresses were higher there.

Conclusions
The effects of the friction coefficient and the loading direction, on the stress distributions around the screw hole in the mandible, were investigated analytically in two dimensions. The loading direction with respect to the material axis of the anisotropic bone turned out to have a negligible influence on the stress distributions. The parametric studies of friction suggested that the friction coefficient had a significant influence on both the radial and shear stress distributions and their maximal values. In crushing failure, the radial normal stresses were the dominating ones, and a high degree of bonding, manifested here by a high friction coefficient, is desirable to mitigate the stresses and increase the strength. The stress distributions in each bone layer were also obtained from the average stresses through classical laminate theory. It was concluded that the analytical method was much easier and more general compared with the finite element analysis, which consists of tedious work such as creating the geometry, meshing, simulating, etc., and it could allow for the parametric study of the design parameters given that a two-dimensional approach is viable. In this case, the present model may be implemented in a larger framework to predict failure based on a suitable criterion. The long-term goal is the quantitative design of improved implants from a mechanical perspective. However, one should keep in mind that the proposed analytical method is only for in-plane loading of flat parts on the bone, where the effects of the finite dimensions, bending moments, and out-of-plane loads were not accounted for.

Acknowledgments:
The authors would like to thank César Antonio Pérez Aranda for technical support and diligent editorial work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The expressions of the coefficients in the expressions of σ rr and τ rθ are provided as follows: where the coefficients A 1 , B 1 , C 1 , D 1 , A c , B c , A ϕ , and tan θ 1 and different terms in the expressions of σ θθ are expressed in detail in [18].