key: cord-0070433-ly1v5pky
authors: Pais, A. I.; Alves, J. L.; Belinha, J.
title: Using a radial point interpolation meshless method and the finite element method for application of a bio-inspired remodelling algorithm in the design of optimized bone scaffold
date: 2021-11-24
journal: J Braz
DOI: 10.1007/s40430-021-03280-2
sha: 81fe444b1afa6cca999d92276c431a3d7590bd97
doc_id: 70433
cord_uid: ly1v5pky

The design of bone scaffold involves the analysis of stress shielding, which can occur when the Young’s modulus of the implant is higher than the Young’s modulus of the bone it is replacing, leading to bone decay in the surrounding tissue. It is therefore very important that the material is adequately designed to match the properties of the surrounding tissue, allowing an appropriate load transfer. While some approaches exist in the literature exploring functional gradients of material density, there are much less solutions based on biological laws. A homogenized model of gyroid infill obtained with PLA ([Formula: see text] MPa) was obtained through mechanical tests of 3D printed specimens, namely tensile and compression, and the obtained model was implemented in a bone remodelling algorithm. The homogenized law was compared to the results obtained with a bone tissue law to assess the equivalence of density distribution and mechanical properties. Through a radial point interpolation method, it was found that similar density fields were obtained for the gyroid infill and for bone tissue when subject to the same boundary conditions. The finite element method was also used for comparison and validation. With the density field results, the gyroid mechanical behaviour was extrapolated to other materials, and similar stiffness values were obtained for bone tissue and titanium alloy ([Formula: see text] GPa) scaffold, which justify this proposal of gyroid scaffolds for mimicking bone properties.

The gyroid shape, first discovered in 1970 [1] , is based on a mathematical equation and is classified as a triply periodic minimal surface (TPMS). These surfaces present the characteristic of being defined by an implicit function and being periodic in its three main directions. Therefore, it is possible, through the equation parameters to edit the unit cell in a way that fits the necessary constraints. Based on these shapes, it is possible to obtain new lighter materials suitable for structural applications [2, 3] .

Bio-scaffold for tissue engineering is another possible application for these cellular materials. Conventional methods for manufacturing bone scaffold implants include poreforming agent leaching, gas foaming or phase separation methods. It has been shown that scaffolds produced by AM are more advantageous since the mechanical properties are more optimized [4] .

TMPS possesses a mean curvature of zero, and approximately, trabecular bone also presents a null mean curvature [5, 6] . Moreover, gyroid scaffolds are, according to Melchels et al. [7] , more easily wetted and intruded, which results in a more homogeneous cell distribution and deeper cell colonization.

Mainly, TPMS has broad application in the development of bone scaffold [8, 9] due to being capable of developing structures with the adequate mechanical properties to fix large bone defects [10] while avoiding stress shielding. In [6] , it was showed that TPMS scaffolds present comparable properties to those of bone depending on the building material and manufacturing process. Thus, combined with adequate density gradients [2, 3] from topology optimization algorithms, it is feasible to match the mechanical properties of bone patch. Similar conclusions have been obtained in [11] .

Homogenization-based optimization approaches consist of first obtaining the constitutive law which determines the elastic constants of the representative volume element (RVE) as a function of relative density and then running the optimization itself [2, 12, 13] . With this approach, it is not necessary to use high-fidelity discretization of the cellular geometry, which decreases the computational cost of the analysis. These approaches take significant interest in the field of biomechanics where unit cells can be designed in a way that its homogenized properties will fit the properties of the tissue it replaces. For instance, Lin et al. [14] developed an homogenization-based optimization procedure which aimed at obtaining the unit cell matching the necessary mechanical properties. Therefore, it is possible to obtain a homogeneous material law which simplifies the RVE in the form of C C * = f ( r ) where C is the homogeneous material constant, C * is the respective constant of the building material, and r is the relative density of the RVE. In this work, the material law being used to characterize the gyroid infill was obtained experimentally and is the form of E = f ( app ) , that is, the law is defined for the building material.

With regard to the use of bone remodelling and topology algorithms to develop and study bone implants, there are several examples in the literature. In 2007, Harrysson et al. [15] mentioned the possibility of using optimization algorithms to develop non-stochastic mesh implants to improve the properties of the hip stem prosthesis. Ghaziani et al. [16] studied a bone growth through a bone remodelling algorithm to test the design of several functionally graded material (FGM) prostheses. The algorithm uses a strain energy density (SED)-based approach, where the density is updated according to the stimulus and Wu et al. [17] evaluated bone ingrowth in the scaffold as a function of the mechanical stimulus. These approaches are concerned with the osteointegration of bone in the implant, by considering previously defined scaffold designs. Wu et al. [18] used topology optimization to develop bone plates with optimal design to maximize bone density where the bone remodelling step is followed by an optimization step. On a different approach, Chuah et al. [19] applied topology optimization to design a spinal interbody cage with the aim of reducing stress shielding. This work presents the main difference that it uses the bone remodelling inspired algorithm to design the scaffold itself. The principle behind it is that as in bone remodelling, the material will develop as a function of the mechanical stimulus, and in accordance with its properties. It is proposed that the gyroid foam will follow a gradient similar to the gradients verified in bone tissue.

The use of meshless methods can improve the quality of the solutions, because meshless methods obtain smoother variable fields than the FEM [20] .

For solid mechanics applications, the diffuse element method (DEM) was one of the first meshless methods being developed, using the moving least square approximations [21, 22] . Then, DEM evolved and the element-free Galerkin method [23] was developed. The smooth particle hydrodynamics method [24] (a pioneer particle/meshless method for particle interactions) was modified and originated the reproducing Kernel particle method [25] . The meshless local Petrov-Galerkin [26] method evolved into the method of the finite spheres [27] , the finite-point method [28, 29] , the radial basis function method [30, 31] . The shape functions constructed through these methods lack the Kronecker delta property [20] , and thus, the essential boundary conditions cannot be directly imposed in the system of equations, increasing the computational cost of the analysis. Thus, other methods have been more recently developed in order to overcome this obstacle such as: the point interpolation method [32] , the point assembly method [33] , the radial point interpolation method (RPIM) [34, 35] , the meshless finite element method [36] , the natural neighbour finite element method [37, 38] and the natural element method (NEM) [39] [40] [41] , the natural radial element method (NREM) [42] [43] [44] and the Natural neighbour radial point interpolation method (NNRPIM) [45] , which is a combination of the NEM and the RPIM.

Recently, partial differential equations have been solved using a non-local operator based on the variational principle of solving partial differential equations. The main advantage of this method is that no discretization is required [46] . Furthermore, with another meshless approach, it is possible to use deep neural networks (NN) to solve the PDE by using the energy as the loss function of the NN. Using this approach it is also not necessary to define the discretization, only the boundary value problem [47] .

The general procedure to solve any problem using a meshless method can be summarized as following [20] :

1. Discretize the problem domains into nodes. The way this process is conducted will heavily affect the performance and accuracy. 2. Build the background integration mesh to integrate the weak form equations, using a Gaussian integration scheme or other integration numerical techniques. 3. Impose nodal connectivity through influence domains.

Each integration point will possess its own influence domain, which is formed the set of nodes near that point. The size and shape of the influence domain will also heavily affect the performance of the method. 4. Constructing the shape functions from the influence domain of that integration point. 5. Obtaining the field variables using approximation functions. 6. Establishing the nodal-based system of equations by assembling the local nodal matrix into the global equation system matrix. 7. Solving the discrete system of equations Ku = f using an appropriate solver.

As in FEM, the RPIM requires a background integration mesh, which is obtained by first dividing the domain in regular cells. Each cell is filled with integration points according to a Gauss-Legendre quadrature scheme. It is possible to construct the integration mesh with a grid which fits the problem domain and then transforming each cell in an isoparametric cell. Using the isoparametric functions the quadrature points are obtained in Cartesian coordinates (following the same procedure used by FEM to establish the integration points inside each element). The isoparametric weight of the integration point is obtained through multiplication of the inverse of the Jacobian matrix with the determinant of the respective grid cell. More detailed information of the numerical integration procedure can be found in [20] .

Nodal connectivity in the RPIM is established through influence domains. Influence domains are formed by the nodes in the vicinity of a said interest point (which is commonly a integration point). Generally, the influence domain is found by a radial search inside a fixed area or volume. Moreover, the influence domains must overlap to establish nodal connectivity. Due to its simplicity, several meshless methods apply this concept [23, 25, 26, 34, 48] . The performance of the method is heavily dependent on the size and shape of the influence domains. Preferably, each domain should have approximately the same number of nodes. For a 2D problem, the number of n d nodes inside the influence domain should be n d = [9, 16] [20]. The 3D problems analysed in this work, considered an influence domain containing n d = 27.

Considering an integration point x I ∈ ℝ 3 , the influence domain of the integration point is the set of n nodes X I = {x 1 , x 2 , ...x n } ∈ ℝ 3 which are used to calculate the shape function for that said integration point. At any given integration point x I , the variable field u(x I ) , which can be any variable field, is interpolated by:

where a(x I ) and b(x I ) are non-constant coefficients of the radial basis function (RBF), r(x I ) and the polynomial vector p(x I ) , respectively. There are several radial basis functions, in this case, the multi-quadrics RBF (MQ-RBF) is used, which is a function of the Euclidean norm d iI of the distance between the interest point x I and the node x i in the support domain of that same interest point.

The parameter d a is the size coefficient, related to the numerical weight ⌢ W I of the integration point x I [20] . In the RPIM, this parameter is calculated as the average nodal spacing of the nodes in the support domain of x I . The parameters and p are shape parameters which have been optimized for best performance, namely = 1.42 and p = 1.03 [20] .

The polynomial vector will present m terms in its base, the monomials obtained from the Pascal's triangle.

If, for example, a constant polynomial base is used, p(x) = {1} and if a linear polynomial base is used p(x) = {1 x y z} T . This work used a linear polynomial base.

The number of unknowns in t he system Ra(x I ) + Pb(x I ) = u s is n + m so a new set of equations must be added. Hence, it is enforced that P T a(x I ) = 0 so that m equations are included, and an unique solution is obtained. Therefore, the system of equations to be solved to obtain the coefficients is:

where R [n × n] is the radial moment matrix determined with (5), P [n × m] is the polynomial moment matrix determined with (6), Z is a null matrix with [m × m] , u s is a vector with the field function nodal parameters [n × 1] , and z is a null vector [m × 1].

The radial basis function vector and the polynomial vector can be combined in a single line vector. Thus, the system of equations becomes:

being the RPI shape functions the vector (x I ) T , and (x I ) T a by-product with m terms and without physical meaning.

Let us consider a body described by the domain

being that Γ u is the essential boundary, and Γ t is the natural boundary. The equilibrium equations which govern the linear elasticity problem are defined by Eq. (9), where ▽ is the gradient operator, is the Cauchy stress tensor for a kinematically admissible displacement field u , b is the body force per unit volume.

The boundary conditions are given by n = t on Γ t and u = u on Γ u , where u is the prescribed displacement on the essential boundary Γ u , t is the traction on the natural boundary Γ t , and n is the unit outward normal to the boundary domain Ω . The constitutive relation is given by = C ∶ , with C being the material compliance tensor defined for a homogeneous isotropic material and is the small strain tensor, = ▽ s u , where ▽ s is the symmetric gradient operator. The governing weak form equation and associated boundary conditions can be presented as

where v(x) is the test function. Through the Galerkin procedure, the weak form for the discrete linear-elastic problem can be written as:

In the Galerkin weak form, the trial and test functions, respectively, u h and u h are represented as linear combinations of the interpolation functions and thus, for an interest point

From Eq. (11), the discrete system of equations Ku = f is obtained, where the stiffness matrix is defined as

For an isotropic material, considering Ω ⊂ ℝ 3 and using the Voigt notation, 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.

Finally, since the interpolation functions used in the meshless method present the Kronecker delta property, the boundary conditions can be imposed directly.

Before explaining the remodelling process, it is important to explain how the mechanical properties of bone are dependent on the apparent density. The apparent density of bone is given by:

where app is the apparent density, w sample is the wet mineralized mass of a said sample, and V sample is the volume said

sample presents. The apparent density can also be presented in terms of the porosity p as following:

where 0 is the density of compact bone, about 2.1 g/cm 3 and the porosity p is given by V holes ∕V sample , being V holes the volume of holes.

The mechanical properties, namely the elasticity modulus and the ultimate stress of bone tissue, can be correlated with the value of apparent density. In addition, bone tissue presents anisotropy so the correlations will differ for the axial and transversal directions.

Lotz et al. [49] described the properties of bone using power laws in cortical and trabecular bone with distinct expressions.

Belinha et al. [20] considering a transition density of 1.3 g/cm 3 to differentiate between trabecular and cortical bone, correlated the mechanical properties of bone tissue with the apparent density using polynomial laws.

Since the proposed algorithm adapts a bone tissue remodelling algorithm, first the fundamentals of Carter's model [50, 51] for bone tissue are presented. Carter's model states that the functional adaptation of bone to the applied stress works an optimization tool by maximizing structural integrity for a configuration with minimal mass. Hence, it is assumed that bone remodelling occurs as a response to a mechanical stimulus S n . This stimulus must be quasi-stationary, constant in the entire bone, and proportional to the stress. which considers l load cases and for each m i load cycles, and for each load case, the exponent measures the influence of the magnitude. Thus, i depends on the apparent density and stress state.

Alternatively, it is possible to use a strain energy density (SED) approach instead of a stress approach, obtaining a correlation between the apparent density and the local strain energy U.

It is assumed in this model that the trabecular arrangement occurs in the direction of principal stress which corresponds to an optimal configuration for stiffness maximization. This model was adapted by Belinha et al. [52] and combined with meshless methods. The bone remodelling phenomenon is

nonlinear. It is presented as a differential equation where a temporal-spatial-based functional app (x, t) is minimized with respect to time:

being that app (x, t) ∶ ℝ d+1 ↦ ℝ is defined for the spatial d-dimensions and the temporal one-dimension. The spacial domain is discretized into N nodes X = {x 1 , x 2 , ...x N } ∈ Ω , and thus, originating Q integration

For a fictitious time step t j , the medium apparent density of the model is given by:

where app I is the apparent density of the interest point x I , defined through I = g I , being g I ∶ ℝ 3 ↦ ℝ .

The stresses 1 , 2 and 3 are the principal stresses, and which define the set of points being subject to the remodelling process, where U m is the minimum U and U M is the maximum U verified. Growth and decay rate of the apparent density are represented by and , respectively, and must be defined for each problem, as well as the control density control app . When equilibrium is achieved, the remodelling stops.

Initially, on iteration i = 0 , linear-elastic analysis is performed considering the constitutive matrix from the building material with the objective of aligning the constitutive matrix with the principal stresses. In order to aid in the visualization of the algorithm, a flowchart is shown in Fig. 1 .

For each iteration i >= 1

1. Linear-elastic analysis is run on the structure. More information on this process can be found in [20] . The displacement field in the structure is obtained, for each load case j. The field variables (stress 

In 3.1, it was explained how the mechanical properties of the bone tissue are related to the apparent density [50, 50, [53] [54] [55] [56] . It is possible to extend these correlations to other materials, namely cellular materials, since the material is porous and its properties will vary with the apparent density [2, 57, 58] .

An isotropic material law to characterize the behaviour of the PLA gyroid infill was obtained using the results of several experimental tests performed in the work of Silva et al. [58] . The isotropy assumption relies on the triple periodicity of the gyroid shape and makes the analysis more simple.

The variation of infill density corresponds to varying the number of gyroid units per cube and thus, the apparent density of each cube (Fig. 2) . Several [20 × 20 × 20] mm 3 cubic specimens made of PLA with distinct gyroid infill, as well as specimens with [150 × 10 × 4] mm 3 , were produced by fused filament fabrication (FFF) and experimentally tested with uniaxial compression and tension tests, respectively [58] .

With the experimental results of tensile tests, it was possible to obtain a homogenised Young's modulus for each specimen, which was then associated with the apparent density of the specimen, considering that the density of bulk PLA is 0 = 1.25 g/cm 3 . In the end, it was possible to develop a material law for gyroid foam cells correlating the Young's modulus with the apparent density With the results of the experimental compression tests, it was also possible to obtain the ultimate compression stress and correlate such value with the corresponding apparent density In the algorithm, it is with (25) that the inverse functions referred in (20) are obtained, and thus, the apparent density in the points can be actualized. The mechanical properties to obtain the field variables are obtained with (24) . In order to allow a direct comparison of the mechanical behaviour of a structure fabricated with gyroid foam and the same structure made of bone, the following bone tissue phenomenological law was considered [20] :

which was simplified to consider bone tissue as isotropic. As mentioned previously, the density level which marks the threshold between cortical and trabecular bone is 1.3 g/cm 3 . Figure 3 shows the plot of the material laws for the gyroid infill and bone tissue, shown in Eqs. (24) to (28) .

Two different models of the femur head are going to be compared, a voxelized simplified model (model V) discretized into a structured mesh of 3086 nodes and 2284 hexahedral 8-node elements (Fig. 4a) , and a realistic model discretized into an unstructured mesh (model R) of 1969 nodes and 9112 tetrahedral 4-node elements (Fig. 4b) .

To summarize the RPIM model, it was defined that the number of nodes inside each influence domain is n d = 27 nodes, that means, each integration point will have ( associated to it, 27 nodes. The cells in the background integration mesh are defined similarly to the elements in the finite element mesh. Thus, for the voxel model, it was considered an 2 × 2 × 2 integration scheme for each cell (thus having 8 integration points inside each cell), and for the realistic model, the cells were tetrahedral and the integration scheme consisted of one integration point in the barycenter of each cell. Relatively to the construction of the shape functions a linear polynomial base p(x) was considered, and the shape parameters of the RBF are = 1.42 and p = 1.03 , as it was mentioned previously, have already been defined for optimal performance of the method. The applied loads aim to replicate the most prevalent load case at the proximal femur [20] : a compression load at the femur head (from the hip joint) and the traction load at the greater trochanter (from the trochanter inserted tendons). At each iteration, the design variables are brought to their elastic limit, through (25) and (28) , and thus, unitary forces were applied, being that the upwards vertical force corresponds to 1/3 of the downwards vertical force. For simplicity, the forces are applied as being vertical (only applied along Oz). It was considered that the nodes in the base of NR2 are pinned.

Both models, V and R, consist of an area which will undergo the remodelling process and another area which is not allowed to be remodelled. The non-remodelling area (NR1) at the top of the models and the non-remodelling area (NR2) at the base of the models are aimed to assure that the boundary conditions (in which displacement constrains and force are applied) remain solid material until the end of the analysis.

At each iteration, the points eligible to be remodelled, identified in Fig. 4 as R zone are the points with density values lower than the maximum density defined by the user. The mechanical properties of the NR zone are E PLA = 3145 MPa, equivalent to rigid PLA for the gyroid law and E bone = 30GPa , equivalent 2x2x2 gyroids 4x4x4 gyroids 10% density 20% density to very rigid bone for the bone tissue material law. Moreover, in order to define the points whose density is changed at each iteration, using (21) , the and parameters must be defined to have the growth and decay rate for the process. As well as the mesh discretization, these parameters affect the quality of the trabecular definition. The remodelling parameters used in this work are shown in Table 1 . The used values are obtained from the literature [20] and adjusted to each model so that the remodelling process occurs. Finally, it was also defined how many times an integration point could update its density. For both models, it was opted to allow each point to remodel only 5 times (which speeds the global computational analysis without significantly the output). The maximum density defined for the gyroid infill corresponds to an infill density of 85%, and the minimum density corresponds to an infill density of 10%, because the density of the filament which was used to characterize the gyroid infill is 0 = 1.25 g/cm 3 . The maximum, minimum and transition density values defined for the bone model are correspondent to threshold values for the properties of bone [59] .

First, it is analysed how the algorithm creates optimized density fields and the evolution along the iterations. Thus, Fig. 5 shows the evolution of the apparent density app for several values of average apparent density avg app , where N is the number of nodes in the area subject to remodelling, R zone (Fig. 4) .

As the iterative process progress, the density field evolves into two main structures, joined by a thinner bar. This pattern was verified for both the voxel model and the realistic model. When avg app reaches 0.5 g/cm 3 , this thin bar (which joins the two main bars) disappears.

In order to compare the mechanical behaviour of the gyroid infill structures, a stiffness coefficient K is calculated for each analysis,

where n f if the number of nodes in NR1 possessing a reaction force, f z is the component of the force in the Oz direction, n b is the number of nodes in NR1, and d z is the displacement in the Oz direction. The force value is calculated at each step as:

where is the strain tensor in Voigt notation, is the stress tensor in Voigt notation, and u is the displacement vector for a design domain Ω.

where B is the deformability matrix given by where N is the shape function. Thus, it is shown in Fig. 6 how the stiffness coefficient K varies with avg app ∕ 0 . The average displacement is used, since using the maximum displacement can lead to misleading results, because of eventual numerical singularities. In model V, the RPIM calculated higher values of stiffness and for model R the highest values of stiffness coefficient were calculated with the FEM. Due to the properties of the homogenized gyroid model, the stiffness coefficient first decreases slightly and then decreases more abruptly. Thus, when avg app ∕ 0 reaches 0.6, it is observed a high rate of stiffness loss in the material as the apparent density decreases. This threshold is equivalent to avg app = 0.75 g/cm 3 , so, when comparing with Fig. 5 , it will correspond to the iterative step where the structure starts to split into the two main bars, relying on the strength of the thinner bar.

Next, it is evaluated how both material laws will converge into similar structures, for equivalent values of avg app ∕ 0 . This will mean that the porosity configuration will be similar. Figure 7 is shown the apparent density field for structures with 50% of the initial mass, for both materials. Therefore, it can be evaluated if the functional gradients, naturally occurring in bone tissue, can be mimicked by other materials. This is important since bone tissue scaffold design must also include the study of permeability and porosity (which are important for cell proliferation) and thus, viability of the implant [6, 10, 60] . Again, it is visible that the structure evolves into two large bars joined by a thin bar. Since the remodelling algorithm acts as a topology optimization algorithm, the formed structures will present configurations that maximize strength for that load case, hence, the two large trabeculae. The horizontal distance between the load application in NR1 and NR2 generates flexural moment and meaning that strength will depend on the thin bar. Finally, it is studied the possibility of using other materials to manufacture gyroid scaffold with the appropriate mechanical properties. Assuming that the displacements observed in any two materials (defined by a bi-linear elastic-plastic law) can be obtained by a proportionality constant; then, the stiffness coefficient of a titanium scaffold can be approximated from the stiffness coefficient of a PLA scaffold within the elastic regime. For this, the following expression was used to extrapolate the mechanical properties of titanium scaffold: Therefore, using (34) , the mechanical properties of the gyroid infill, adjusted to the Young's modulus of titanium [61] ( E Ti = 110 GPa), were obtained, which are shown in Fig. 8 . This material is often used as a bio-material due to

its bio-compatibility. One of the most important aspects in scaffold design is developing a material with the mechanical properties matching the mechanical properties of the bone tissue it is replacing. In order to avoid stress shielding, the apparent stiffness cannot be higher than the apparent stiffness of bone. However, the material should present enough load-bearing capability to withstand the naturally occurring load cases in the human body [6] . Therefore, in Fig. 8 , the stiffness curves are compared to the stiffness curves of the bone tissue models.

Similarly to the scaffold model, the bone model also showed a higher rate of stiffness loss when avg app ∕ 0 reaches 0.6. In terms of scaffold adequateness, stiffness is lower than the stiffness of bone, for equivalent porosity values. However, the approximation in (34) is not possible to be used to approximate load-bearing capability, because the material strength is independent of the modulus, i.e. a material can present high modulus and low yield stress. Therefore, the next step in the design could be the evaluation of load-bearing capacity of the scaffold. Thus, for values of avg app ∕ 0 = 0.5 , it is evaluated using the FEM and the RPIM, the von Mises stress in both models considering the bone tissue law (Fig. 9) , which presents the stress values that the hypothetical titanium scaffold should be able to withstand. It can be seen in Fig. 9a that with Model V, the stress level measured with the FEM is lower than the stress level measured with the RPIM. During the remodelling process, the applied force is updated at each remodelling stage. Using (21) the set of points being remodelled is chosen, and based on the stress level, the variables are brought to the elastic limit. During the remodelling process, some points can present a high stress level and still be under the elastic limit, while presenting higher stress than the overall structure. This will update the applied force in a way that a much lower force is applied in the following iteration steps. The stiffness coefficient will remain the same because the lower force will generate a proportionally lower displacement. Figure 10 shows that in model V, from iteration 5 forward, the value of force in the FEM is lower. It can be seen that the stress map at iteration 1 is coherent for both methods, as well as the values of F z from iterations 1 to 5.

Due to (25) , which updates the density values in the points, it is not possible to apply a proportionality law (as in (34) ) to obtain the stress field for the titanium scaffold. Thus, the von Mises equivalent stress field obtained with the gyroid PLA is not presented. Thus, the titanium scaffold should be able to withstand at least the 700 MPa which occur in bone tissue. According to the work by Bobbert et al. [6] , the titanium scaffolds show enough yield stress. However, since the scaffold geometry tested by Bobbert is not the same as the scaffold geometry tested by Silva et al. [58] it is not possible to obtain conclusions.

This work proposes a remodelling algorithm to optimize density fields, allowing to achieve improved functional gradients foams for bone tissue scaffolds (locally mimicking the mechanical properties of bone tissue).

Several remodelling analyses using a radial point interpolation meshless method were performed on two different models and considering two distinct material laws. For comparison purposes, the same analyses were performed with the FEM.

The meshless method (RPIM) proved to be capable of achieving solutions of density gradient to aid in the design of bone scaffold, with structure configurations that are similar to the ones obtained with the FEM (validating the numerical tool). It was found that although geometrically less accurate, the model with the voxel (structured) mesh obtained correct values for the stiffness of the scaffold model and bone model. Whether the voxel model or realistic model was used, similar structure configurations were obtained, adapted to the load case.

Due to computational hardware limitations, the mesh density used for this work was not enough to obtain very well-defined trabecular structures. However, it was possible to visualize two main high-density bars, coherent with the literature [20] and medical imaging. Considering that in a near future, this technique is to be combined with additive manufacturing, an extremely intricate trabecular structure is to complex to be transformed into a gradient part. So, for that reason, the level of detail obtained can be enough to design the scaffold, without wasting computational power in an extremely detailed discretization Finally, the hypothetical titanium scaffold obtains scaffolds with lower density than bone which will prevent stress shielding. However, it should be mentioned that the gyroid infill is an approximation of the gyroid geometry and the manufacturing process will influence the mechanical properties, because the titanium scaffold cannot be manufactured the exact same way as the PLA infill.

Infinite periodic minimal surfaces without selfintersections. National Aeronautics and Space Administration

Optimal design and modeling of gyroid-based functionally graded cellular structures for additive manufacturing

Interior structural optimization based on the density-variable shape modeling of 3D printed objects

Applications of 3D printed bone tissue engineering scaffolds in the stem cell field

A three-dimensional model for tissue deposition on complex surfaces

Additively manufactured metallic porous biomaterials based on minimal surfaces: a unique combination of topological, mechanical, and mass transport properties

Mathematically defined tissue engineering scaffold architectures prepared by stereolithography

Biomechanical influence of structural variation strategies on functionally graded scaffolds constructed with triply periodic minimal surface

Laser additive manufacturing of bio-inspired lattice structure: forming quality, microstructure and energy absorption behavior

Mechanical behaviours and mass transport properties of bone-mimicking scaffolds consisted of gyroid structures manufactured using selective laser melting

3D-printed ceramic triply periodic minimal surface structures for design of functionally graded bone implants

Cellular structures from additive processes: design, homogenization and experimental validation

Bending behavior of optimally graded 3D printed cellular beams

A novel method for biomaterial scaffold internal architecture design to match bone elastic properties with desired porosity

Direct metal fabrication of titanium implants with tailored materials and mechanical properties using electron beam melting technology

The effect of functionally graded materials on bone remodeling around osseointegrated trans-femoral prostheses

A time-dependent mechanobiology-based topology optimization to enhance bone growth in tissue scaffolds

Time-dependent topology optimization of bone plates considering bone remodeling

Topology optimisation of spinal interbody cage for reducing stress shielding effect

Meshless methods in biomechanics-bone tissue remodelling analysis

Generalizing the finite element method: diffuse approximation and diffuse elements

Surfaces generated by moving least squares methods

Element-free Galerkin methods

Smoothed particle hydrodynamics: theory and application to non-spherical stars

Reproducing kernel particle methods

A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics

The method of finite spheres

A finite point method for elasticity problems

A finite point method in computational mechanics. Applications to convective transport and fluid flow

Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-I surface approximations and partial derivative estimates

Multiquadrics-a scattered data approximation scheme with applications to computational fluid-dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations

A point interpolation method for stress analysis for two-dimensional solids

A point assembly method for stress analysis for two-dimensional solids

A point interpolation meshless method based on radial basis functions

On the optimal shape parameters of radial basis functions used for 2-D meshless methods

Del Pin F (2003) The meshless finite element method

Natural neighbour finite elements

Natural neighbour Galerkin methods

A numerical method for solving partial differential equations on highly irregular evolving grids

The natural element method in solid mechanics

Imposing essential boundary conditions in the natural element method by means of densityscaled -shapes

The natural radial element method

Analysis of thick plates by the natural radial element method

Composite laminated plate analysis using the natural radial element method

Analysis of 3D solids using the natural neighbour radial point interpolation method

A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem

An energy approach to the solution of partial differential equations in computational mechanics via machine learning: concepts, implementation and applications

Meshless methods: a review and computer implementation aspects

Mechanical properties of metaphyseal bone in the proximal femur

The compressive behavior of bone as a two-phase porous structure

Mechanical properties and composition of cortical bone

A meshless microscale bone tissue trabecular remodelling analysis considering a new anisotropic bone tissue material law

The mechanical behaviour of cancellous bone

The mechanical properties of trabecular bone: dependence on anatomic location and function

On the dependence of the elasticity and strength of cancellous bone on apparent density

Determinants of the mechanical properties of bones

Design of functionally graded gyroid foams using optimization algorithms and the finite element method

Study on 3D printing of gyroid based structures for superior structural behaviour

Some basic relationships between density values in cancellous and cortical bone

Bone tissue engineering: hope vs hype

Properties: Titanium alloys-ti6al4v grade 5

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations