The Biomechanical Analysis of Tibial Implants Using Meshless Methods: Stress and Bone Tissue Remodeling Analysis

: Total knee arthroplasty (TKA) stands out as one of the most widely employed surgical procedures, establishing itself as the preferred method for addressing advanced osteoarthritis of the knee. However, current knee prostheses require refined design solutions. This research work focuses on a computational analysis of both the mechanical behavior of a knee joint implant and the bone remodeling process in the tibia following implantation. This research study delves into how specific design parameters, particularly the stem geometry, impact the prosthesis’s performance. Utilizing a computed tomography scan of a tibia, various TKA configurations were simulated to conduct analyses employing advanced discretization techniques, such as the finite element method (FEM) and the radial point interpolation method (RPIM). The findings reveal that the introduction of the implant leads to a marginal increase in the stress values within the tibia, accompanied by a reduction in the displacement field values. The insertion of the longest tested implant increased the maximum stress from 5.0705 MPa to 6.1584 MPa, leading to a displacement reduction from 0.016 mm to 0.0142 mm. Finally, by combining the FEM with a bone remodeling algorithm, the bone remodeling process of the tibia due to an implant insertion was simulated.


Introduction
Total knee arthroplasty (TKA) stands as one of the frequently conducted surgical interventions in the field of Orthopedics, offering the promise of enhancing functionality, alleviating pain, and reinstating the quality of life for patients [1].It is estimated that in the USA, the growth in surgical requests from 2005 to 2030 will be around 673% (3.48 million) [1].The main objective of a TKA is to reduce knee pain, based on the replacement of most of the damaged components of the joint with an implant, and the success of the surgery is estimated by the absence of pain and functional recovery of the knee in a short period [2].TKA can encompass the surface of up to three bones: the femur, the patella, and the tibia; however, the latter is the area where most injuries leading to this surgical procedure occur [2].In certain cases, for enhanced implant fixation and stability, a cylindrical extension known as the tibial stem is incorporated at the distal end of the tibial metal component whose primary function is to increase implant stability, particularly in challenging cases, by minimizing micromovements at the bone-implant interface, thereby mitigating the risk of aseptic loosening [3].These tibial extensions come in various lengths and widths and can be secured through different methods, including total cementation, proximal cementation only, or a press-fit approach without the use of cement [3].The existence and design of the tibial extension, including its size and method of fixation, can induce a phenomenon referred to as stress shielding, resulting in bone loss in the regions most affected, which may weaken the implant fixation, requiring the placement of a new prosthesis [4][5][6][7].Stems are fundamental in most TKAs in aiding the transfer of loads from the damaged articular and metaphyseal bone to the tibial cortical limit and in distributing the increased stress of a specific joint [8].Stems, at the cost of stress shielding, improve mechanical stability through shearing resistance, reduction in lift-off, and a decrease in micromovement [8].In cases where the available bone reserve is inadequate to sustain the prosthesis, the incorporation of stems is generally advised [9].When addressing substantial volume defects through bone grafting, the use of a stem becomes essential to shield the graft from excessive loads, since the knee joint bears loads several times the body weight and a failure in load transfer by the stem can subject the remaining trabecular bone to a load surpassing its maximum strength, resulting in a compromise in the fixation of the component [9].While there is today a large consensus about the necessity of a stem for enhancing both the initial mechanical stability and the ultimate longevity of the component, the recommended lengths and diameters and fixation methods continue to be subjects of controversy [8].It is noteworthy that it is highly unlikely that an ideal length, for example, will ever be determined, due to a huge heterogeneity in the anatomical characteristics of patients [8,10].
In cemented arthroplasty, the oldest fixation method, the cement ensures the union between the implant and the bone tissue, providing a uniform load distribution throughout the extent of the bone-implant interface [11].However, the cement, usually composed of polymethyl methacrylate (PMMA), does not provide adhesive properties; it only fills the gaps located between the bone and the prosthesis [11].Many studies have indicated that cemented fixation shows better short-and medium-term results, as it reduces patient pain and increases mobility and, consequently, their quality of life [11].This method is mainly employed in cases of extensive bone damage and a significantly compromised internal cortex, particularly observed in individuals diagnosed with advanced osteoporosis [8].However, some incidents with the use of cement have also been reported, as it degrades over time, causing the release of debris that leads to inflammation of the surrounding bone [12].Another point to consider is the surgical intervention of revision, recurrent in younger patients, since the removal of the cement-fixed prosthesis is quite complicated, and therefore, this technique is then recommended only for patients over 65 years of age or patients with poor bone condition, in order to ensure a strong primary fixation [12,13].The adversities and complexity related to cemented arthroplasty have encouraged the search for new fixation methodologies.
Unlike cemented arthroplasty, which results from the mechanical fixation of the cement, uncemented arthroplasty is based not only on mechanical fixation but also on the biological junction of the implant to the bone tissue [14].This new technique aims to improve, in the long term, the success of implants in young patients, who demonstrate good bone quality.The main drawback of this approach is the pain in the lower limbs, which is a consequence of the poor fixation of the implant to the bone tissue [12].
Computational tools such as the FEM have been very useful for researchers in the field of arthroplasty.Studies using these tools focus on the implant geometry and materials as well as surgical techniques.With these studies, issues such as the minimization of stress shielding or micromotions, which can heavily hinder the success of the treatment, are dealt with.Completo et al. [15] conducted a study using the finite element method, which evaluated the distribution of loads and stability at the cement-bone interface for two fixation techniques, cemented and press-fit.The findings revealed that the load transmitted from the cemented stem to the bone was fourfold higher [15].Quevedo González et al. [16] used a computational approach to study the effect on tibial stress of using a smaller tibial baseplate in order to prevent excessive rotation, having concluded that the consequences of this undersizing are minimal.Quevedo Gonzalez et al. [17] used a computational approach to determine that adding an anterior spike to the tibial baseplate would help to decrease the micromotions between the bone and the implant.Other approaches addressing implant geometry include the work of Liu et al. [18], who used a topology-optimization approach to improve the design of the metal plate, through a porous design, used to treat defects in the proximal tibia with TKA.Liu et al. [19] compared two techniques used in TKA, namely, the cement screw and the metal block, and concluded that the use of a metal block is more suitable for larger defects.Regarding the study of the most suitable materials, studies such as those of Bhandarkar and Dhatrak [20] and Apostolopoulos et al. [21] deal with the cushion material and tibial insert, respectively.In the first study, it is concluded that using a UHMWPE cushion leads to lower stresses in the inserts, and in the second work, an all-polyethylene TKA is compared to a metal-backed TKA and it is stated that the all-polyethylene one is cheaper and yields similar results to the metal one, with the disadvantage of requiring cement and having less flexibility.Finally, some studies concerning the arthroplasty of the ankle can also be mentioned, such as that of Jyoti and Ghosh [22] concerning the optimal resection length, which concluded that stress shielding and micromotion increase with resection depth, and thus, that this should be kept to a minimum; and the work of Jyoti et al. [23] concerning the interface of the implant, in which it is proposed that both the friction coefficient between the implant and the bone, and the implant design, along with the bone quality, heavily influence micromotion, while stress distribution is less influenced by the friction coefficient and more by the design of the implant and the quality of the bone.
The present computational approach aims at complementing experimental studies and international standards such as ISO DIS 22926 [24] and ISO/CD 5092:2023 [25].Furthermore, this research employs meshless methods to achieve higher accuracy in obtaining solutions, thereby enhancing the quality of computational studies on implant biomechanics.Moreover, it incorporates a bone remodeling approach to assess if the implant could cause considerable bone loss as a result of stress shielding.

Meshless Methods Formulation
Generally, in meshless methods the process initiates with the discretizaton of the problem's physical domain with a nodal distribution.Although, as in other discretization techniques, regular nodal distributions allow more accurate and stable results to be achieved, in meshless methods nodes can be irregularly distributed over the problem domain.Geometry features such as cracks and holes can have a higher nodal density to anticipate the stress concentration that may occur.Usually, a very dense nodal distribution will lead to more accurate results, with the consequence of presenting a much higher computational cost [26].
Next, to integrate the integro-differential equations governing the physical phenomenon under study, it is necessary to create a background integration mesh.Commonly, background integration cells are created and then filled with integration points representing volume portions.Most meshless methods requiring integration apply the Gauss-Legendre integration scheme to the integration cells [26].
The previous stage is followed by the imposition of nodal connectivity through influence domains.Thus, each integration point searches for a certain number of nodes within its vicinity (generally, following a radial search).The nodes within its vicinity form the influence-domain of the integration point.The literature shows that the size of the influence domain significantly affects the performance of the method.
Then, using the nodes forming the corresponding influence domain, the shape functions (and their partial derivatives) of a given integration point are built, and matrix B, from Equation (18), can be established.
Taking the example of the displacement field u, the displacement at an integration point x I can be interpolated using the displacement values of the nodes within the influence domain of x I , where the number of nodes inside the influence domain of x I is represented as n d , the displacement values of each node i belonging to the influence domain of x I is indicated as u(x i ), and φ i (x I ) is the value on node i of the shape functions of x I .
For each integration point x I , a local stiffness matrix is established, K I , Equation (22).Then, all local stiffness matrices are assembled into a global stiffness matrix, allowing the final discrete system of equations Ku = F to be obtained, which can be solved to obtain the displacement field u = K −1 F.

Influence Domains
The RPIM uses the concept of influence domains, which is shown schematically in Figure 1.The influence domain consists of the set of nodes in the vicinity of the interest point.The nodes within the spatial vicinity of integration point x I constitute the influence domain of x I .Influence domains can all have the same number of nodes or instead present different sizes and different numbers of nodes.Moreover, the shape of the influence domain can be any.In Figure 1a, the influence domains present the same size (r 1 = r 2 ), which leads to a different number of nodes in each influence domain, while in Figure 1b the number of nodes inside each influence domain is the same due to their different sizes: (r 1 ̸ = r 2 ).

Shape Functions
Consider the integration point x I ∈ R 3 and its influence domain: with n being the total number of nodes within the influence domain of x I .The shape function of x I is calculated using only the n nodes belonging to its influence domain.Thus, the variable field u(x I ) can be interpolated at x I using where a(x I ) and b(x I ) are non-constant coefficients of the radial basis function (RBF), r(x I ), and of the polynomial basis function (PBF), p(x I ), respectively.Although the literature describes several distinct RBFs [27], generally, radial point interpolation techniques apply the multi-quadrics RBF (MQ-RBF) [28]: The parameter d c is the size coefficient, directly correlated with the numerical weight w I of the integration point x I , and the parameter d iI is the Euclidean distance between the integration point x I and the node x i inside the influence domain of x I .The parameters p and γ are shape parameters and, according to [27], these should present values so that γd a is almost null and p is almost unity, which maximizes the method's performance.
The PBF possesses m terms (the monomials), which can be established using Pascal's triangle.

p(x
In radial point interpolation formulations, the constant PBF p(x) = {1} or the linear PBF p(x) = {1 x y z} T is generally used.
The number of unknowns in the system Ra(x I ) + Pb(x I ) = u s is n + m.Therefore, to build a system of equations it is necessary to include a new set of equations.The literature shows that it is necessary to impose P T a(x I ) = 0 to achieve a unique solution [27].This leads to the following system of equations: where R [n × n] is the radial moment matrix (Equation ( 6)), P [n × m] is the polynomial moment matrix, (Equation ( 7)), and u s is a vector with the field function nodal parameters By back-substitution in Equation ( 2), it is possible to obtain In which ϕ(x I ) T represents the RPI shape function vector, and ψ(x I ) T is a neglectable vector, a by-product of the additional set of equations: P T a(x I ) = 0.
Finally, the RPI shape functions verify the Kronecker delta property, which allows for the direct imposition of boundary conditions, as they pass through all the nodes in the influence domain: as well as satisfying the partition of unity:

Weak Form and Discrete System of Equations
Assuming well-known essential boundaries and the initial and final time conditions (i.e., the compatibility conditions), the energy principle dictates that out of all possible displacement configurations satisfying such compatibility conditions, the unique final solution is the one minimizing the Lagrangian functional.The minimization of the Lagrangian functional for a solid with domain Ω and boundary Γ can be represented as where ρ is the solid-mass density, u is the velocity, ε is the strain tensor, σ is the stress tensor, u is the displacement vector, b are the body forces, and t are the traction forces applied to the boundary Γ ∈ Ω.The first term of ( 12) refers to the kinetic energy and is, therefore, discarded in static problems.The second term is referent to the strain energy and the last two terms are referent to the work produced by the external forces.The variational operand can be moved inside the integral and the second integral in ( 12) is rearranged.Thus, for Equation ( 12) to be valid for the compatibility conditions, its integrand must be equal to zero.Such an equality results in the Galerkin weak form.
where δε is the virtual strain tensor and δu is the virtual displacement.
The principle of virtual work states that if a solid body is in equilibrium, the virtual work produced by the inner stresses and the body applied external forces are null when the body experiments a virtual displacement.
The energy conservation principle states that if the work produced by the inner stresses is equal to the work produced by the applied external forces, then the body is in equilibrium.If the work produced by the applied external forces is calculated with a virtual displacement field and the work produced by the inner stresses is calculated using a virtual strain field resultant from the same virtual displacement field, then the principle of virtual work is obtained.Thus, knowing that stresses and strains are related through the generalized Hooke's law equation, Equation ( 14), σ = Dε (14) and strains are linearly related to displacements, Equation (15), the weak form Equation ( 13) can be re-written in terms of displacement as follows: Inverting the material compliance matrix C, it is possible to obtain the constitutive material matrix D = C −1 .Assuming the Voigt notation, the material compliance matrix C can be presented as where c 1 = 1/E and c 3 = −ν/E, being that E and ν are the Young's modulus and the Poisson ratio of the material, respectively.In Equation ( 16), B is the deformation matrix: so that the strain and displacement relation is for any integration point x I .H is a diagonal matrix containing the shape functions φ(x I ).
The virtual displacement can be removed from (16), obtaining the discrete system of equations.Equation ( 20) can be written in the matrix form as where the elemental stiffness matrix K e is and the force vector F e is given by Finally, the global system of equations is assembled from the element matrices so that the displacement field can be solved.

Bone Tissue Remodeling Algorithm
Bone remodeling is a vital process for bone health and maintenance, allowing bones to develop, repair, and renew throughout life.This process is complex and highly coordinated, involving the removal and replacement of deteriorated bone tissue through the activity of various cellular components.These components ensure mineral homeostasis and structural integrity [29].The cellular groups that sequentially carry out the process of bone resorption and formation compose the bone remodeling unit, which has been termed the basic multicellular unit (BMU) [30].Initially, a BMU is formed by osteoclasts responsible for bone resorption.Subsequently, the bone surface is covered by reversal cells and prepared for bone replacement.In the next step, osteoblasts release and deposit osteoid, the unmineralized bone matrix [29].
Several studies [31][32][33][34] affirm the correlation between the mechanical properties of bones and their apparent density.Additionally, various bone tissue models [35][36][37][38] in the literature delineate the evolution of trabecular architecture based on mechanical stimuli, employing diverse formulations and assumptions.
The apparent density of bone is given by where ρ app corresponds to the apparent density, w sample represents the wet mineralized mass of a given sample, and V sample represents the volume of the same sample.Alternatively, it is also possible to represent this property as a function of its porosity p: where ρ 0 corresponds to the density of compact bone, approximately 2.1 g/cm 3 , and porosity p is calculated from V holes /V sample , where V holes is the total volume of holes.The apparent density holds significance in its connection to the mechanical attributes of bone tissue, specifically the elasticity modulus and ultimate stress.Belinha [26] employed polynomial laws to establish a relationship between the mechanical properties of bone tissue and apparent density, incorporating a transition density of 1.3 g/cm 3 to differentiate between trabecular and cortical bone.
Thus, in this work the following phenomenological law was considered [26]: In Equations ( 26) to ( 28), E bone cortical corresponds to the elastic modulus of cortical bone, E bone trabecular corresponds to the elastic modulus of trabecular bone, and σ bone c corresponds to the ultimate compressive stress.
In the bone remodeling process, there is strong evidence of the correlation between the loads (stress or strain) that the area is subject to, and the shape that the bone takes.This has led to the requirement to develop some semi-empirical laws in order to model how the bone will functionally adapt to the load case it is subject to.These laws are the basis for computational tools which aim at predicting bone adaptation under stress or strain states or changes in stiffness [26].The model employed in this study [39][40][41][42][43] assumes that, regardless of the material law applied, the induced stress serves as an optimization tool.It aims to maximize structural integrity while minimizing mass.This perspective is analogous to viewing induced stress as an optimization tool, where the objective function is minimized, as discussed by Belinha [26].
Figure 2 summarizes the bone remodeling iterative algorithm.
for each load case i

Weighted critical variable field
Point subject to remodeling? yes Update density of point x i using the material law for all interest points The bone remodeling nonlinear equation is presented as a differential Equation ( 29) in which a temporal-spatial based functional, ρ app (x, t), is minimized with respect to time, and where ρ app (x, t) : R (d+1) → R is defined for the one temporal dimension and the d spatial dimensions.
It is assumed that the d-dimensional domain is discretized in N nodes, X = {x 1 , x 2 , . . .x N } ∈ Ω, leading to Q integration points, where The temporal domain is discretized in iterative fictitious time steps t j ∈ R, where j ∈ N. Within the same iterative time step t j , the average apparent density of the iteration point x I is defined by ρ I = g(σ I ): Then, g(σ I ) : R 3 → R is defined as in (31).
Here, σ j are the principal stresses obtained for integration point x I and σ −1 j (ρ I ) are the inverse functions of σ j (ρ I ), defined with the material law.
The inverse equation of stress as a function of apparent density ( 28) is then applied to the integration points with von Mises stress values inside the following interval, where The parameters α and β define the growth rate and the decay rate of the apparent density.The remodeling equilibrium is achieved when Equation ( 33) is verified, where α, β and ρ control app depend on the analyzed problem.
In summary, at each iteration the linear elastic analysis is run on the model.The elastic properties of the bone material are obtained from the density at each point using (26) and (27).With the obtained variable fields, the critical values for the variables are determined.The points with higher values of stress will have their density increased and the points with lower stress will have their density decreased, being that the new density value is determined through (28).

Tibia and Implant Numerical Models
In this section, the construction of 3D models for the numerical applications analyzed in this study is presented.Note that all the calculations were performed using an original code developed by the authors using the Matlab© (Natick, MA, USA: The MathWorks Inc.) environment.All the routines and images produced were programmed by the authors without using any toolbox of Matlab©.Nevertheless, the 3D models were created in Autodesk Fusion© (San Francisco, CA, USA: Autodesk) and the 3D discretizatons were produced with ABAQUS© (Johnston, RI, USA: Dassault Systemes Simulia Corp).
The stress and displacement analysis of the proximal tibia was run in three-dimensional models discretized into 4-node tetrahedral elements.The analysis of osteointegration to the implant was run on a two-dimensional model in order to allow for a finer mesh and, consequently, a finer trabecular arrangement.The two-dimensional model was discretized into 3-node plane-stress elements.Table 1 summarizes the number of nodes and elements used for each model.The knee joint is subjected to loads that vary significantly with the activity the individual performs.The loads applied in the model took into account another study, which analyzed the gait of an individual, where it was concluded that the average peak force was 1645 N [18].According to the literature, the axial force is considered the primary point of evaluation of the model's response, since it was previously determined as an indicator of injury risk.Thus, it was found that the axial load was 987 N and 658 N for the medial and lateral platforms of the tibia, respectively, since the load ratio between the medial and lateral platforms is 60%:40% [18,44].It is important to note that the applied loads do not take into account the forces of the ligaments, muscles, and tendons.The tibia is kept in balance by the surrounding tissues and other bones of the lower limbs.The difficulty in representing this system in a finite element model leads to simplifications in terms of boundary conditions.Therefore, the nodes located at the most distal part of the model were fixed, effectively constraining any rigid-body movement.
Figure 3a shows the boundary conditions applied to the three-dimensional model without the implant.Equivalent boundary conditions were applied to the implant model as well.The loads labeled as F1 and F2 were distributed among the nodes in the vicinity of the positions indicated in Figure 3a, with moduli of 658 N and 987 N, respectively.For the remodeling analysis, the same two forces were considered.However, in order to more accurately emulate the shear loading which occurs during locomotion, these two loads were applied considering a 10º angle, as shown in Figure 3b.Concerning the models featuring an implant, these were developed by converting the mesh model into a spline model, enabling the precise modeling of the implant's shape.When the implant comprised three components, as depicted in Figure 4a, the model was transformed into a spline model as illustrated in Figure 4b.This simplification resulted in a unified structure with distinct material properties assigned to the bone and implant areas, as visible in Figure 4c.Healthy cortical bone has an elastic modulus of 17 GPa; however, in this study, this property was varied (Table 2) in order to check the impact of this change on the distribution of the effective von Mises stress, as well as on the displacement field.The implant was considered to be Ti-6Al-4V.
The initial implant model shown in Figure 4a was modeled based on the the commercial implant by DePuy Synthes®.

Structural Analysis of the Proximal Tibia
The Von Mises stress does not depend on the Young's modulus of the material but it will differ between numerical methods.Thus, Figure 5 shows a set of points to evaluate the von Mises stress and the respective results.
The displacements, however, will depend on the mechanical properties.Figure 6 shows the displacement field results for the three conditions and both methods.Because the results for both methods are very similar, only one of them is shown.
Additionally, the displacement was analyzed at two points in the proximal part of the tibia, labeled in Figure 5a as point A and B. For the three materials and both methods, the displacement values are shown in Table 3.

Structural Analysis of the Implant and Influence of Implant Length
The Von Mises stress fields for all the implant lengths considering healthy bone only are shown in Figure 7.
Concerning the stress analysis of the model without the implant, the FEM analysis concluded that the highest stress value occurs for point 2 while the RPIM analysis concluded that the highest stress value occurs for point 3.This phenomenon is due to the fact that the stress field is a field that has been extrapolated to the nodes, which means that the stress value is being smoothed at the nodes.By increasing the number of nodes in the mesh, this effect starts to diminish, and the maximum stress begins to approach the base of the fixture.
The maximum stress values verified after the insertion of the implant, between 5.0568 and 5.9497 MPa for the FEM and between 4.4690 and 6.1584 for the RPIM, are similar to the values verified at the unimplanted model of 5.0705 MPa for the FEM and 4.7880 MPa for the RPIM model, which indicate that stress shielding would not occur.The maximum stress values occur at similar locations for the implanted and unimplanted models.Stress is an important measure to evaluate the occurrence of stress shielding (SS) as it indicates whether adequate loading is being applied to the bone.It can be quantified as the stress change before and after implant insertion [45] (34): where σ bone and σ implantedbone denote the Von Mises stress of the bone without and with the implant, respectively.Considering the maximum Von Mises stress verified at all stem lengths, the occurrence of SS can be discarded since the stress value increases after insertion.The displacement fields for all implant lengths are shown in Figure 8.The maximum displacement, which occurs at the medial side, is similar with and without the implant.The unimplanted model, taking into consideration the healthy bone, is 0.016 mm, while for the FEM this value ranges between 0.0142 and 0.0151 mm and 0.0156 and 0.0171 mm for the RPIM.It should be noted that the mesh model for the tibia varies with each implant length, which could account for some of the discrepancies observed between the two methods as the essential and natural boundary conditions were not being applied to the exact same nodal coordinates.Thus, excluding the RPIM analysis of the 30 mm length implant, the stress verified at the tibia would increase with implant length.
The greater the stiffness of the implant, the less load will be transferred to the bone from the deformation suffered by the stem [45].The evaluation of displacements suffered by the implant is, therefore, an important measure of stress shielding.Since it was verified that the implant displacements are in the range of the displacements of the intact bone, similar stiffness values were achieved.
The measured displacements with the RPIM are higher than the FEM while the stress values are lower.It is expected that the results obtained through RPIM will be closer to the actual results, after conducting a convergence study (which helps to ensure that the results are more accurate), or if a denser mesh is generated, as was the case in the study developed by Marques et al. [46].
In conclusion, the tested stem sizes for an implant in the selected material are not likely to induce stress shielding.Nevertheless, this study does not consider bone heterogeneities, which alter the mechanical properties and consequently the structural response.

Bone Remodeling Analysis
Regarding the osteointegration evaluation, two analyses were considered.First, the medial and lateral forces were considered as two separate forces acting at two different moments, and for the second case, those forces acted simultaneously.Figure 9a shows the apparent density maps for the two tested load combinations.Only a finite element analysis was considered for this stage.
Concerning the bone remodeling analysis, there is no significant difference between the imposition of the medial and lateral loads simultaneously or as different loads.Regarding the trabecular arrangement, there are three major bone resorption areas, as indicated by the numbers in Figure 9b.Additionally the main trabecula formed by the algorithm is aligned with the load angle as expected, since the algorithm states that higher stress areas are subject to having their density increased and lower stress areas are subject to having their density decreased, leading to alignment of the density map with the stress map.

Conclusions
Utilizing computational approaches enhances the study of implant design and surgical techniques alongside clinical and experimental research.In summary, longer stem lengths, such as 40 mm, result in higher maximum stress in the proximal tibia but lead to lower displacements compared to shorter lengths.However, the 30 mm model stood out as an outlier in the maximum displacement due to the use of different models for testing each length.Additionally, bone remodeling analysis confirmed the expected anatomical structures, including areas of bone resorption and the diaphysis, through the observed density distribution.

Figure 2 .
Figure 2. Flowchart for the bone remodeling algorithm.

Figure 3 .
Figure 3. Schematic representation of the essential and natural boundary conditions applied to (a) the three-dimensional model and (b) the two-dimensional implant model.

Figure 5 .
Figure 5. Von Mises stress results: (a) selected points in the model (A and B are two representative points of load application and, therefore, are used to evaluate the displacement and points 1 to 9 are the points selected to evaluate stress since each point is further from the fixed ending); (b) von Mises stress at each point calculated with the FEM and the RPIM.

Figure 7 .
Figure 7. Von Mises stress field at the nodes for the tested implants.

Figure 8 .
Figure 8. Displacement field at the nodes for the tested implants.

Figure 9 .
Figure 9. (a) Evaluation of osteointegration using the FEM along the iterations; (b) indication of the main resorption (indicated by numbers 1 to 3) and growth areas after implant insertion.

Table 1 .
Mesh information for each model.

Table 2 .
Mechanical properties of the implant material and cortical bone.

Table 3 .
Displacement at two points calculated with the FEM and RPIM for all tested materials.