Interpretation of parameters in strain energy density bone adaptation equation when applied to topology optimization of inert structures

The internal structure of bones is naturally design to fulfil their mechanical function with a minimum mass. The structure is continuously adapted and maintained during life through a local control process known as bone remodelling (BR). Created to model this process, there are now several mathematical theories, e.g. [1, 2], which, in conjunction with the finite element method (FEM), proved to be able of generating structures resembling the trabecular bone morphology. In essence, in the case of isotropic material assumption, BR simulations aim to redistribute the material within a given analysis domain, by varying the local density, such that higher density is associated with higher loads. This goal is achieved using control parameters which are defined based on assumptions and experimental evidence regarding the mechanobiology of bones. In structural design, the same goal of determining an optimum material distribution within a certain domain, given the boundary conditions and loads, is achieved using topology optimization (TO) techniques [3-5]. The optimum distribution is defined relative to an objective, such is the minimization of compliance. Although developed independently, the two techniques mentioned above share significant similarities. Particularly, the similarities between the strain energy density (SED) approach (discussed in paragraph 2) for BR and the solid isotropic material with penalization (SIMP) (synthetized in paragraph 3) for TO, have been investigated. A detailed comparison of the two methods can be found in the work of Jang et al. [6]. Also, Jang et al., demonstrated that, for certain conditions, there is a mathematical analogy between SED based BR and SIMP based TO [6]. Nowak [7] shown the usefulness of applying BR theory to structural design optimization. Two main advantages were emphasized: the independence of the optimization process of the design domain and the ability of allowing multiple load case optimization, a common issue in TO. According to the above arguments, the BR SED based theory is a potential candidate for developing algorithms that aim to optimize inert structures. However, when applied to structural optimization, some parameters that describe the BR equations lose their original significance, related to bone mechanobiology. This paper aims to give an interpretation of such parameters, in the case of the SED based BR equation with a spatial influence function (SIF). The purpose is to facilitate the understanding of the BR mathematical model when applied to the optimization of continuum inert structures and to evaluate whether values of the parameters which are not necessarily related to some biological aspect, can lead to optimum solutions important from the mechanical design perspective. A 2D numerical example is given, in order to test the influence of the interpreted parameters upon the topological final solutions.


Introduction
The internal structure of bones is naturally design to fulfil their mechanical function with a minimum mass.The structure is continuously adapted and maintained during life through a local control process known as bone remodelling (BR).Created to model this process, there are now several mathematical theories, e.g.[1,2], which, in conjunction with the finite element method (FEM), proved to be able of generating structures resembling the trabecular bone morphology.In essence, in the case of isotropic material assumption, BR simulations aim to redistribute the material within a given analysis domain, by varying the local density, such that higher density is associated with higher loads.This goal is achieved using control parameters which are defined based on assumptions and experimental evidence regarding the mechanobiology of bones.
In structural design, the same goal of determining an optimum material distribution within a certain domain, given the boundary conditions and loads, is achieved using topology optimization (TO) techniques [3][4][5].The optimum distribution is defined relative to an objective, such is the minimization of compliance.
Although developed independently, the two techniques mentioned above share significant similarities.Particularly, the similarities between the strain energy density (SED) approach (discussed in paragraph 2) for BR and the solid isotropic material with penalization (SIMP) (synthetized in paragraph 3) for TO, have been investigated.A detailed comparison of the two methods can be found in the work of Jang et al. [6].Also, Jang et al., demonstrated that, for certain conditions, there is a mathematical analogy between SED based BR and SIMP based TO [6].
Nowak [7] shown the usefulness of applying BR theory to structural design optimization.Two main advantages were emphasized: the independence of the optimization process of the design domain and the ability of allowing multiple load case optimization, a common issue in TO.
According to the above arguments, the BR SED based theory is a potential candidate for developing algorithms that aim to optimize inert structures.However, when applied to structural optimization, some parameters that describe the BR equations lose their original significance, related to bone mechanobiology.This paper aims to give an interpretation of such parameters, in the case of the SED based BR equation with a spatial influence function (SIF).The purpose is to facilitate the understanding of the BR mathematical model when applied to the optimization of continuum inert structures and to evaluate whether values of the parameters which are not necessarily related to some biological aspect, can lead to optimum solutions im-portant from the mechanical design perspective.A 2D numerical example is given, in order to test the influence of the interpreted parameters upon the topological final solutions.

The BR equation with SIF
In the simplest form, the BR models are expressed by an ordinary differential equation of the form: where dρ / dt designates the derivative of a mass function (commonly, the apparent density) with respect to time, B is a constant meant to regulate the speed of the process, S(x,t) is the control variable at the location x, thought as the tissue level mechanical stimulus, which is somehow sensed by specialized cells and So is the stimulus value corresponding to a homeostatic equilibrium state, at which level no remodeling occurs.The density is only allowed to vary between zero value and a maximum value which is either the cortical bone density or unity, if the relative density is used [1].
The bone material behaviour is considered linearelastic, in most BR studies.The Young modulus is iteratively updated based on its correlation with the local density, via an expression of the form: with a and m empirically determined constants.A review of different expressions can be found in [8].Eq. ( 2) indicate that there is a dependence of the structure stiffness matrix on the material elastic properties, via the local bone density.
When coupled with FEM, Eq. ( 1) is usually solved per each finite element (FE).It is assumed that there is one sensor point (where the mechanical stimulus is evaluated) per each FE.
In order to estimate the local load, S(x, t), Huiskes et al. [9] used the SED.Later, Weinnans et al. [10] investigated models defined by Eq. (1) using as mechanical stimulus the strain energy per unit of mass (SEM), determined from: where U(x, t) and ρ(x, t) represent the SED and the apparent density at the location x.They proved that BR equa-tions of the form Eq. ( 1) are mesh dependent and checkerboard pattern generators.To overcome these issues, Mulender et al. [11] proposed a model based on the previous version described, generally, by Eq. ( 1), but decoupling the sensor points from the FE mesh.In order to model the range of bone sensor cells (osteocytes) influence they assumed an exponential decay of the signal with respect to the distance to the remodeling location, according to the function: where di(x) is the distance from the sensor i to the location x and D represents the distance at which the signal is attenuated by 36.8 % (e -1 ).Using the SIF (4) the remodeling equation is written: with N being the total number of the sensors within the vicinity of the location x and Si(t) the SEM at the location of the sensor i.The authors shown that, within certain limits of the mesh density which depend on problem dimensionality, the model described by Eq. ( 5) is mesh independent.It also eliminates the checkerboard effect.

The SIMP optimization model
As in BR, in the case of TO, finding the optimum layout of the material is also based on the iteratively change in local rigidity, evaluated at each FE level via the element stiffness matrix.The optimum solution is achieved by correlating the stiffness matrix of the structure with a vector of design variables, which are used to search the solution.In the SIMP approach, the stiffness matrix, K, depends on a fictitious density vector, ρ, according to [12]: where NE is the number of the finite elements in the design domain and Ke is the stiffness matrix of the finite element e, corresponding to a unit value of density, ρe and q represents a so called penalty exponent, which has the role of diminishing towards zero (penalizing) the subunit densities.The basic idea beyond the stiffness matrixdensity relation, is the correlation of the Young modulus with the density, following [13]: with Ee and Eo are the Young moduli of the element e and of the base material (with unit density), respectively.
According to [13], the term "density" is used because the total volume of the structure is calculated as: with Ve being the volume of the finite element e.The volume, V, is thought as a cost function, on which a constrained is imposed, as an upper bound, Vmax.
Based on the parameterized stiffness matrix from Eq. ( 6), for a linear elastic discrete structure, the equilibrium equation is written: where F and u are the vectors of nodal forces and nodal displacements, respectively.Unlike BR problem, which is formulated using ordinary differential equations, the TO problem lies on the minimization of an objective functional, g, as follows [12]: A commonly used objective has the form [14]: where the first term represents the structural compliance while the second is the volume of the structure.The relative importance of the two terms is weighted by the positive coefficient, µ.

Methods and interpretation of BR parameters
On the basics of the Eq. ( 5), it was developed a code that couples MATLAB R7 (MathWorks, Inc.) and ANSYS Mechanical APDL R15 (ANSYS, Inc.) simulation environments.The remodeling equation is integrated using forward Euler scheme.Although known as imprecise, it was proved [15] that, in BR simulations, if a proper integration step is used, it generates enough accurate density distributions.Thus, at iteration n+1, the density update is given by: where h is the integration step.The remodelling constant B has no significance with respect to the inert structures optimization, because there is no relevance to correlate the speed of the process with some biological aspect.Therefore, it is convenient to couple it with the integration step h, by choosing a unit value.Thus, Bh = h.Eq. ( 5) needs an initial condition, i.e. an initial value, ρo, of the density.In BR simulations, it is usually adopted half of the maximum value, determined by the cortical bone density.In this paper, the initial condition is defined by ρo = 1, i.e. the maximum relative density, considering that the purpose of the TO algorithm is to reduce the material within a given domain.Hence, the starting value should correspond to the material from which to reduce.It is however important to verify whether the initial condition affect the results in terms of final density distribution, including the shape and the total mass.According to Weinnans et al. [10], by changing the initial value of the density similar solutions are obtained but not identical.The same conclusion was formulated by Mulender et al. [1].They also noticed that the mass remains constant, while the shape changes.
Due to the lack of empirical Young-density relations for inert materials, the Young modulus can only be expressed as a function of relative density.As described in paragraph 3, such an approach is applied within the SIMP method and is adopted herein.Therefore, the elasticity modulus is iteratively updated based on the relation: where Eo is the Young modulus corresponding to the material with the unit relative density, thus the material defining the final topology.According to Bendsoe and Sigmund [16], in the case of SIMP method, the minimum admissible value of m should be 3, such that the intermediate Young moduli, E, to satisfy the Hashin-Shtrikman bounds for a composite material with one phase being void.
Another important parameter of the remodeling equation is the threshold So.For bones, it corresponds to a level of local loading that keeps the remodeling process inactive [17].Its significance should be adapted to the purpose of achieving a structure with minimum material consumption.Hence, assuming that the optimum material distribution include only areas of unit densities, the resulting structure is defined by FE of constant rigidity, having the Young modulus equal to Eo.Under the linearelasticity assumption, So can be determined from the stress-strain curve of the material with the Young modulus equal to Eo.For an allowable limit of stress, σa, using the von Misses criteria, the target value of the mechanical stimulus is given by: with υ being the Poisson's ratio of the material.Because Eq. ( 5) is thought to drive the model towards a relative uniformly distributed SED in FE with Eo, the value of So also determines the value of maximum stress, playing thus the role of a global stress constrained.generally, one cannot expect a final topology characterized by uniformity of SED distribution and thus uniformity of von Mises stress distribution.It is already known from BR [10] that the threshold value cannot be achieved in all the elements.But a von Mises stress distribution with values approaching σa is expected.
In the case of application to inert structures optimization, the osteocytes influence function also loses its original significance.The choice of the D parameter is no more correlated with the thickness of the trabeculae [11].Its importance, however, remains since it contributes to the elimination of the checker boarding and to generation of a mesh independent solution.As pointed out in the original paper [11], the parameter D should be taken such that at least one FE dimension to be covered.In the case discussed here, its variation may lead to optimal topologies that are structurally relevant.
The density at the location x is regulated by the SED determined in all the elements within a certain vicinity.The elements for which the distance to the location significantly exceeds the value of D are not important, because their influence decreases towards zero.Therefore, it is useful not to take into account the contribution of these elements in order to reduce the computational cost.In this respect, a new parameter is introduced, named radius of influence (RIF), that defines the vicinity of location x containing active elements, i.e. elements which contribute to the calculation of SED at location x.Denoting the RIF with R, than the minimum value of the SIF is: which leads to a correlation between R and D given by: In order to control the number of active FE, one can either choose to modify R or fmin.The first parameter explicitly controls the number FE vicinities contributing to the SED calculation, while the latter gives the percent reduction of the contribution.
In BR studies [1], during the simulation process, the elements for which the density decreased to the minimum value were excluded from SED contribution.The reason is based on the idea that the mechanical signal comes from osteocytes, which can be found within the bony tissue.In the case discussed herein, the same approach is adopted.Although no osteocytes are present, excluding the low density elements reduces the computational time.More, in these elements the stresses tend to zero due to their low rigidity and very small deformations occur.Hence, the SED also tends to zero, becoming irrelevant.
For a clear presentation of the connection between the BR mathematical model and the FEM, in Fig. 1 is presented a general block diagram of the implemented algorithm.In essence, the diagram shows how BR algorithm is applied to each finite element of the discretized structure.

Fig. 1 General block diagram of the simulation algorithm
Under the linear elasticity state equation and following the algorithm block diagram presented in Fig. 1, one can deduce that stresses and strains, at the finite element level, play the role of problem state variables.

Numerical example
The theory presented above is tested on a common structure (Fig. 2) in TO literature [18,19], in order to have a suited reference for verifying the validity of results.
A hypothetic material is considered with a Young modulus, Eo, of 10 GPa corresponding to the unit relative density.The intermediate Young moduli are estimated according to Eq. ( 14), with m = 3.The Poisson's ratio is taken equal with 0.3.Assuming an allowable stress, σa, of 100 MPa, from Eq. ( 15) it follows that So = 0.43 MPa.The kinematic properties of the used finite elements are based on small displacements assumption.
The plate dimensions used within the simulations are L = 50 mm and H = 20 mm.A mapped mesh is generated using four nodded quadrilaters.All the elements are congruent squares having the edge length of 0.5 mm.The center of each element also represents a sensor point, but the sensors are decoupled from the FE mesh, as described in paper [11].A plane stress state is imposed.The loading force is plp = 400 N, selected such that the stresses far from the nodes were supports and loads are applied, not to exceed σa.In order to reduce the stress concentration effect, the force is distributed over several nodes selected on the length lp = 4 mm.First, a value of fmin is selected.In this respect, the contribution to the mechanical stimulus calculation of elements for which the value of SIF is below 0.05 is not considered.Therefore, fmin = 0.05 (corresponding to the distance to which the influence has reduced to 5%) which, according to Eq. ( 17), leads to: Note that relation (18) gives a continuum correlation between the two parameters, covering any combination of sensors and finite elements.Because the algorithm localizes the finite elements and the sensors by their center, the values of R should be multiples of the distance between two consecutive centers.Hence, for the Mitchel plate mesh adopted here, R should be either a multiple of 0.5 mm or a multiple of 0.72 mm.The second value accounts for the diagonal neighbors.Following these observations, for each selected D different values of R are imposed as multiples of 0.5 and 0.72.Then, one solution is selected on the basics of density distribution, uniformity of von Mises stresses and mass reduction, as discussed in the next section.For the selected solution, the effect of the other parameters mentioned above is investigated.

Results and discussions
In Table 1, there are presented several variants of topologies, obtained for different combinations of R and D. The results are given in terms of final density distribution, where the white finite elements correspond to minimum allowed density, ρmin = 0.01, while the black ones correspond to maximum density, ρmax = 1.All the solutions were determined with h = 0.01.
One can notice the similarity between the solutions from Table 1 and the ones presented in the literature using other methods for TO [18,19].The first topology in Table 1 is achieved for values of D and R that reduce the contribution to SED calculation to be equivalent with the model without a SIF, described generally by Eq. ( 1).The checkerboard effect appears and the final topology tends to a structure with many branches.As known from Mulender et al. [11], increasing the parameter D results in reducing the number of branches and in the elimination of the checkerboard patterns.Note that, for a given value of D, the number of branches is also reduced by increasing R. It is interesting that the final mass reduction is around 70% for all the configurations.The differences consist in shape rather than in mass quantity.
Increasing the parameter R over the value given by 3D does not significantly change the final density distribution (note solutions 3 and 4 or 7 and 8), but implies more computational time due to higher number of neighbors that has to be identified and accounted for within the SED calculation process.
An important result is that similar solutions can be achieved for different combinations of R and D. Compare, for instance, solutions 3 and 5 or 8 and 9.By increasing D and decreasing R, the number of contributing elements is reduced, but their influence is increased.Thus, the mechanical signal, which is the sum of the contributions, can be the same for different values of R and D. The fundamental advantage of this result consist in the computational cost reduction.For instance, the time for solving the solution 9 was reduced with 14% compared to solution 7 and with 26% with respect to solution 8.
The von Mises stresses for the initial homogenous plate are represented in Fig. 3.The scale is limited to maximum 100 MPa, in order to emphasize the stresses far from the boundary conditions areas.The maximum value above this limit is of 230 MP and is located in the nodes were load is applied.The von Mises stresses in the solutions given in Table 1 vary between 70 MPa and 400 MPa.But the relevant values should be considered on the median paths on each strut defining the final topology.Due to unrealistic local modeling of boundary conditions, in the corresponding areas the stresses are excessively high and should be neglected.Another local stress concentrations appear on the edges of the resulting struts, in some notches between consecutive rigid elements.These stress raisers are also irrelevant being a result of FE discretization.If the local high stresses are neglected, than within the final density configurations presented in Table 1 (except the configuration 1), the von Mises stresses vary in the range 70-160 MPa, as shown in Fig. 4 for the solutions 5 (a) and 9 (b).
In order to visualize only the areas with stresses around 100 MPa, the color map scale is limited to the interval 80 MPa -120 MPa, for both the solutions 5 and 9.The resulting von Mises plots are shown in Fig. 5.
From Fig. 5 it can be noticed that the solution 9 distributes the loads within the material more uniformly.For instance, within the superior strut the stresses are clos-er to 100 MPa in the solution 9 than they are in the 5 one.It is concluded that, the smaller the number of density struts, the smaller the stress variation is recorded.Diminishing the number of struts by keeping the mass constant, the algorithm generates solutions with thicker struts, in which the stresses are reduced towards the targeted value.The results in terms of von Mises stresses, show that, by choosing the value for the threshold according to expression (15) one can determine solutions with a stress distribution around an admissible value.However, it is not ensured that the admissible stress value is constantly distributed over the entire final topology.
Because the mass reduction is the same for all solutions, the criteria for selecting one for further simulations is based on the configuration, stress uniformity and computational time.As discussed above, the solutions with fewer and thicker struts ensure better uniformity and smaller values in stress distribution.Between similar topologies, as the solutions 7, 8, 9 are, the best is considered the one ob-tained with the smallest computational time, thus, the solution 9 (the smallest value of parameter R).Using this solution, the influence of the other parameters is further investigated.It is to be mentioned that a rigorous measure of time needed for each computation to be finalized is not performed, because of the dependence on the computer configuration and resources allocation.Thus, this study is intended to demonstrate the efficiency of the RIF in terms of relative time reduction, not to provide quantitative time information.
Several simulations are performed based on the same conditions of solution 9, but changing the initial relative density, ρo.In Table 2 there are presented corresponding final density distributions.Except for the initial density values within the range 0.4-0.6,all the others determined similar solutions.The final mass reduction is the same for all variants.The result shows, however, that some sensitivity of the mathematical model to initial conditions exists.From the structural point of view it is important, as different configurations can be obtained, starting from different initial conditions.On the other hand, if the initial density is, for instance, 0.1, relative to this value, the final mass increases.Thus an optimum solution from the mass reduction perspective is no more achieved.Such an interpretation is however only formal because an initial condition, from the practical point of view should correspond to the material density from which to reduce.More, low values of the initial density are associated with low initial Young modulus.In this context, large deformations occur and the linearity assumption is violated.

Solution
In Table 3 there are presented the solutions determined with different values of So, in order to verify its influence upon the final configuration.It can be noticed that, the So influence only consists of modifying the struts thickness.No change in shape of the final topology is determined.Therefore, So has the role of controlling the total final mass and the maximum stress, as a higher threshold leads to lower stresses over the entire model.
It is to be noted that, the calculation of the threshold So based on the relation (15) assumes that the applied load determines stresses in the initial structure in the range of σa.If the load is either to low or to high relative to σa, than the value So is rather indicative.Depending on the topology obtained, So can be adapted to correspond to the load.

Conclusions
The paper presents an interpretation of the control parameters that define the SED based BR equation, in order to adapt their significance to topology optimization of inert structures.The parameters discussed are the homeostatic equilibrium threshold, the constants involved in densityelasticity modulus expressions, the initial density and the distance attenuation parameter.A new parameter is also introduced in order to control the number of finite elements that participate to the calculation of the mechanical signal for each spatial location.The parameter is useful for reducing the computational time.
It is shown that, neglecting the significance of the BR parameters, originally related to bone mechanobiology, other values become valid.In this respect, the mathematical model of BR leads to optimal structural topologies useful in mechanical design.
For combinations of D and R within the range of FE dimension, the final mass tends to be similar.However, different values of the parameters D and R can produce solutions with different configuration.Therefore, one cannot expect on the existence of a unique solution from the final shape perspective.
The value of the mechanical stimulus threshold is not enough to control the stress distribution.But controlling the parameters R and D, the model can be oriented towards a solution with a more uniformly distributed stress over the entire topology.

Fig. 2
Fig. 2 Test structure used for TOThe influence of different parameters is investigated, including D and R, the initial condition ρo and the threshold So.Each simulation is performed until no significant total mass changes are registered.The integration step is adapted to each group of parameters after several numerical experiments.First, a value of fmin is selected.In this respect, the contribution to the mechanical stimulus calculation of elements for which the value of SIF is below 0.05 is not considered.Therefore, fmin = 0.05 (corresponding to the distance to which the influence has reduced to 5%) which, according to Eq. (17), leads to:

Fig. 3
Fig. 3 The initial von Mises stress distribution with the scale limited to 100 MPa (the value of σa )

Table 2
Structural topologies obtained for different values of the initial density, ρo.The parameters D and R correspond to solution 9 from Table1.Above each solution, the initial relative density is indicated

Table 3
Structural topologies obtained for different values of the threshold So