Ivan Giorgio
email: ivan.giorgio@uniroma1.it

Ugo Andreaus

Daria Scerrato

A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials

Keywords: Mechanical-biological coupling, Bone functional adaptation, Continuum mixture model, Poroviscoelasticity, Dissipation Mathematics Subject Classification (2000) 74L15 Biomechanical solid mechanics

In this paper, the phenomena of resorption and growth of bone tissue and resorption of the bio-material inside a bi-component system are studied by means of a numerical method based on finite elements. The material behavior is described by a poro-viscoelastic model with infiltrated voids. The mechanical stimulus that drives these processes is a linear combination of density of strain energy and viscous dissipation. The external excitation is represented by a bending load slowly variable with sinusoidal law characterized by different frequencies. Investigated aspects are the influence of the load frequency, of type of the stimulus and of the effective porosity on the time evolution of the mass densities of considered system.

Introduction

The first scientific papers about the use of vibrations for therapeutic purposes in human bones date back to 1949, when [START_REF] Toscani | Modification of the effects of immobilization upon metabolic and physiologic functions of normal men by the use of an oscillating bed[END_REF] reported the positive effects obtained through the application of vibrations generated by a special bed swinging on bedridden patients in cast immobilization, suffering from metabolic abnormalities. The ability of the bone tissue of changing its conformation depending on the intensity and frequency of the load has been investigated by numerous researchers over the past years [START_REF] Turner | Mechanotransduction in bone: do bone cells act as sensors of fluid flow?[END_REF]. In particular, it is seen that in the cortical bone subjected to bending a fluid passage occurs from the compressed areas to those in tension, thus determining the onset of the fluid flow which induces shear stress on osteocyte, therefore leading to the mechano-transduction. Bone is a mechanically sensitive biological tissue, which adapts its size, shape, mass and density based on its mechanical environment. It has been long recognized that a dynamic stimulus is required for bone adaptation [START_REF] Lanyon | Static vs dynamic loads as an influence on bone remodelling[END_REF][START_REF] Rubin | Regulation of bone mass by mechanical strain magnitude[END_REF][START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF]. Researchers have shown that the adaptation depends on a combination of different mechanical stimuli such as the magnitude [START_REF] Burr | Effects of biomechanical stress on bones in animals[END_REF][START_REF] Lanyon | Mechanically adaptive bone remodelling[END_REF] and frequency of applied load [START_REF] Burr | Effects of biomechanical stress on bones in animals[END_REF][START_REF] Lanyon | Mechanically adaptive bone remodelling[END_REF][START_REF] Hsieh | Effects of loading frequency on mechanically induced bone formation[END_REF][START_REF] Warden | Mechanotransduction in the cortical bone is most efficient at loading frequencies of 5-10 Hz[END_REF], number of cycles [START_REF] Rubin | Regulation of bone formation by applied dynamic loads[END_REF][START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF], and bouts of the applied loading [START_REF] Robling | Improved bone structure and strength after long-term mechanical loading is greatest if loading is separated into short bouts[END_REF][START_REF] Srinivasan | Restinserted loading rapidly amplifies the response of bone to small increases in strain and load cycles[END_REF]. [START_REF] Turner | Mechanotransduction in bone: do bone cells act as sensors of fluid flow?[END_REF] observed significant cortical bone adaptation when the loading frequency exceeded 0.5 Hz. At low frequencies (between 0.5 and 10 Hz), it has been observed that adaptation in rat ulnae follows an approximately linear dose-response relationship with frequency [START_REF] Hsieh | Effects of loading frequency on mechanically induced bone formation[END_REF]. [START_REF] Warden | Mechanotransduction in the cortical bone is most efficient at loading frequencies of 5-10 Hz[END_REF] found no significant increase in the adaptation response when the frequency of loading was increased beyond 10 Hz. [START_REF] Rubin | Anabolism: Low mechanical signals strengthen long bones[END_REF][START_REF] Rubin | Skeletal adaptation to mechanical stimuli in the absence of formation or resorption of bone[END_REF] investigated the effect of very low magnitude high frequency (greater than 30 Hz) loading on sheep. [START_REF] Bentolila | Intracortical remodeling in adult rat long bones after fatigue loading[END_REF] used supraphysiological cyclic axial loading of the ulna of rats in vivo to induce bending with consequent fatigue and microdamage, revealing that rat bone undergoes remodeling in response to high levels of cyclic strain (see e.g. [START_REF] Andreaus | Damage modelling and seismic response of simple degrading systems[END_REF]; [START_REF] Andreaus | Incremental finite element procedure for cumulative damage of two-dimensional continua via two-surface plasticity[END_REF]; [START_REF] Andreaus | Prediction of micromotion initiation of an implanted femur under physiological loads and constraints using the finite element method[END_REF]; [START_REF] Andreaus | Coupling image processing and stress analysis for damage identification in a human premolar tooth[END_REF]; [START_REF] Andreaus | Finite element analysis of the stress state produced by an orthodontic skeletal anchorage system based on miniscrews[END_REF]; [START_REF] Lu | Modeling of an initial stage of bone fracture healing[END_REF]; [START_REF] Placidi | A variational approach for a nonlinear one-dimensional damageelasto-plastic second-gradient continuum model[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF]; [START_REF] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF]; [START_REF] Yang | Higher-order stress-strain theory for damage modeling implemented in an element-free galerkin formulation[END_REF]; [START_REF] Yang | Higher-order continuum theory applied to fracture simulation of nanoscale intergranular glassy film[END_REF] for damage modeling by using second gradient theories particularly suitable for poroelastic materials such as bone). [START_REF] Besdo | Determination of dynamically adapting anisotropic material properties of bone under cyclic loading[END_REF] presented a phenomenological theory aiming to finite element simulation of bone remodelling and to prediction of dynamically changing anisotropic elastic parameters (see also for numerical implementation [START_REF] Federico | An energetic approach to the analysis of anisotropic hyperelastic materials[END_REF]). [START_REF] Tormena | A bone remodelling model based on generalised thermodynamic potentials and optimisation applied to a trabecula with cyclic loading[END_REF] proposed a model based on the thermodynamic framework to describe the process of bone remodeling and modeled a trabecula subjected to cyclic loading and calibrated the model with experimental data.

Cortical tissue is made up of a solid matrix crossed by canals and other porosity hosting vasculature and bone cells. Blood vessels hosted in the Haversian and Volkmann canals, running through the osteon in the longitudinal and transversal directions, respectively, allow the chemical exchanges necessary the transport of tracers and for bone nutrition. Cancellous (or spongy) bone tissue is composed of irregularly arranged trabeculae. In both cortical and spongy tissue, osteocytes are finely dispersed in lacunae of the solid matrix and are connected with each other and with the vascular porosity by cytoplasmic extensions (dendritic processes) running through tiny canals called canaliculi. The lacuno-canalicular network is saturated by an interstitial fluid. In this respect, [START_REF] Sansalone | Interstitial fluid flow within bone canaliculi and electro-chemo-mechanical features of the canalicular milieu[END_REF] numerically assessed the influence of the physical and material parameters affecting the canalicular fluid flow, namely the fluid viscosity, the so-called zeta potential, the ionic concentration, and the peri-cellular matrix permeability. According to [START_REF] Rho | Mechanical properties and the hierarchical structure of bone[END_REF], there are three levels of bone porosity within cortical bone and within the trabeculae of cancellous bone, all containing a fluid. Osteons are cylindrical structures of 100-150 µm radii that lie very roughly along the long axes of a bone and contain at their center an osteonal canal. This canal contains blood vessel(s), a nerve and some space occupied by bone fluid. The canaliculi are passageways that run between lacunae or from the lacunae to the osteonal canal. The porosities include the vascular porosity (PV), the lacunar-canalicular porosity (PLC) and the collagen-apatite porosity (PCA) as well as a porosity that does not contain bone fluid -the porosity of the inter-trabecular space (PIT). There are, therefore, three porosities of bone associated with bone fluid; the fourth level of bone porosity, PIT, contains substances more viscous than bone fluid. The fluid in the space outside the blood vessels and nerves in the Volkmann and osteonal canals is called 'serum' and the fluid in the extracellular space in the lacunae and the canaliculi 'extracellular fluid'. Both of these fluids are called bone fluids here. All the quasi-cylindrical passageways in the bone matrix -the osteonal canals and the Volkmann canals-that contain the vasculature, nerves and bone fluid, are responsible for PV. The characteristic lineal dimension associated with this space is the radius of the Volkmann and osteonal lumens (order 20 µm); by this criterion, this porosity is the largest bone fluid porosity. Bone fluid freely exchanges with the vascular fluids because of the high permeability of the thin endothelium, the absence of a smooth-muscle layer, and the sparse basement membrane layer of the capillary walls. The bone fluid pressure in PV cannot generally exceed the local blood pressure as the vessels would collapse. Thus the pressure of the bone fluids is probably similar to the blood pressure in bone, on the order of 40-60 mm Hg, which is low compared to the bone fluid pressures in the lower characteristic lineal dimension PLC, described below. All the space in the lacunae and the canaliculi determine the PLC. The characteristic lineal dimension associated with this space is the radius of the canaliculus (order 0.1 µm). Since the bone fluid porosity with the largest lineal dimension, PV, is always at a low pressure, the middle characteristic lineal dimension porosity, PLC, appears to be the most important porosity for the consideration of mechanical and mechanosensory effects in bone. Calculations in [START_REF] Cowin | Bone poroelasticity[END_REF] show that the bone fluid in PLC can sustain higher pressures for longer times due to mechanical loading. Thus PV functions as a low-pressure reservoir that interchanges bone fluid with PLC. This interchange is facilitated by the fact that the lineal dimension associated with PV is two orders of magnitude larger than that associated with PLC, and PV is typically at blood pressure, which is low in bone. PLC is the porosity associated with the slower relaxation of the excess pore pressure due to mechanical loading. It is also the porosity associated with osteocytes; the osteocyte is the prime candidate for the mechanosensory cell in bone. PCA, associated with the spaces between the collagen and the crystallites of the mineral apatite (order 10 nm radius), is the lowest characteristic lineal dimension porosity. The movement of the bone fluid in the collagen-apatite porosity is negligible because most of the bone water in that porosity is bound by interaction with the ionic crystal. This portion of the bone water is considered to be part of the collagenapatite structure. The bone porosity with the largest, but quite variable, characteristic lineal dimension (up to 1 mm) is associated with cancellous bone; it is referred to here as PIT of cancellous bone. It is the porosity external to, and surrounding, the trabeculae. This porosity is well connected to the medullary cavity and may contain marrow, fat and blood vessels. These materials have a viscosity one to two orders of magnitude larger than the viscosity of the bone fluid in PV or PLC [START_REF] Bryant | On the mechanical function of marrow in long bones[END_REF][START_REF] Sobotková | Rheological behaviour of bone marrow[END_REF]. The characteristic lineal dimension of PIT varies with anatomical location; it is smaller near the load bearing surfaces and increases to its greatest magnitude as the medullary canal is approached. The permeability associated with this porosity was surveyed by [START_REF] Arramon | Hydraulic stiffening of cancellous bone[END_REF].

It has been shown that osteocytes can more effectively sense the interstitial fluid flow through their processes and respond to fluid-induced shear stress [START_REF] Adachi | Trabecular bone remodelling simulation considering osteocytic response to fluid-induced shear stress[END_REF]. In view of its biological relevance, biochemical and mechanical effects of canalicular fluid flow have been investigated by many authors [START_REF] Tate | Whither flows the fluid in bone?" an osteocyte's perspective[END_REF][START_REF] Fritton | Fluid and solute transport in bone: flow-induced mechanotransduction[END_REF]. Several groups have examined the relationship between mechanical loading on bone and the associated fluid flow inside the various internal porosities. Knothe Tate and Knothe (2000) observed the fluid flow in and out of cortical bone in sheep forearms under applied load. Knothe Tate et al (2000) showed similar load-induced fluid flow in rat tibia subjected to bending, and also indicated the role of the fluid flow in mechano-transduction. Numerous mechanisms have been proposed to explain the effect of fluid flow on the osteocytes e.g., via shear stress on the cell surface [START_REF] Reich | Fluid shear stress as a mediator of osteoblast cyclic adenosine monophosphate production[END_REF], drag force on the transverse fibrils that tether the osteocyte to the canalicular walls [START_REF] Weinbaum | A model for the excitation of osteocytes by mechanical loading-induced bone fluid shear stresses[END_REF], strain generated electric potentials [START_REF] Pollack | An anatomical model for streaming potentials in osteons[END_REF][START_REF] Salzstein | Electromechanical potentials in cortical bone-II. experimental analysis[END_REF], biochemical diffusive gradients [START_REF] Robling | Mechanical stimulation of bone in vivo reduces osteocyte expression of sost/sclerostin[END_REF], or chemical contributions [START_REF] Ganghoffer | Mechanical modeling of growth considering domain variationpart II: Volumetric and surface growth involving eshelby tensors[END_REF]. [START_REF] Fritton | Fluid and solute transport in bone: flow-induced mechanotransduction[END_REF] and [START_REF] Santos | The role of osteocytes in bone mechanotransduction[END_REF] with their extensive reviews provide a more complete description of fluid flow induced mechanotransduction in cortical bone. [START_REF] Kumar | Dissipation energy as a stimulus for cortical bone adaptation[END_REF] presented a finite element study of a poroelastic rectangular beam subjected to oscillatory bending loads. Then they proposed the use of the dissipation energy of the poroelastic flow as a mechanical stimulus for bone adaptation, and showed that it can predict the effect of frequency of the applied load. They showed that the dissipation energy stimulus and the resulting increase in second moment of inertia of the cross section increase linearly with frequency in the low frequency range (less than 10 Hz) and saturate at the higher frequency range (greater than 10 Hz). Similar non-linear adaptation frequency response also has been observed in numerous experiments. [START_REF] Kumar | Modeling of cortical bone adaptation in a rat ulna: Effect of frequency[END_REF] employed the above mentioned model for the adaptation of cortical bone in response to mechanical loading to study the effect of loading frequency on the computed response, and they compared their results to previous experimental measurements on rat ulnae. They represented the cortical bone as a poroelastic material with orthotropic permeability. Bone adaptation in the model is related to a mechanical stimulus derived from the dissipation energy of the poroelastic flow induced by deformation. [START_REF] Baïotto | A study of the viscoelastic effect in a bone remodeling model[END_REF] propose a model of bone remodeling which take viscosity properties of the tissue into account. They employed Zener's law to describe the mechanical behavior of the bone and proposed a specific law of the apparent bone density rate.

The aim of this paper is to develop a numerical model to predict the processes of growth and resporption of bone tissue and resorption of bioresorbable material, modeled as fluidinfiltrated poroelastic materials under repeated loads. In order to achieve this goal, the use of different mechanical stimuli is proposed, namely the strain energy density of the poroelastic material and the dissipation energy of the poroelastic flow, that can predict the effect of frequency of the applied load. Furthermore, a linear combination of the two above mentioned stimuli is considered, which allows for an optimal choice of their relative weights. Osteocytes are interconnected through the processes inside the canalicular space, and it is reasonable to assume that they can communicate with each other through these processes and exchange information on the current state of the stimulus at each location. Thus, adaptation process in the model depends on the weighted average of the mechanical stimulus in a 'zone of influence' near each point, in order to incorporate the non-locality in the mechanotransduction of osteocytes present in the lacunae and to diversify the distribution of porosity experimentally observed in inner and outer parts. It has to be remarked that at least two time scales are of relevance in the considered context: i) the characteristic time of bone growth and resorption and ii) the characteristic time of varying externally applied loads. In the first approximation which we will consider, the second time scale will be assumed to be negligible and mechanical phenomena will be assumed to be quasi-static. In particular, we focus our attention on the remodeling process in a two-dimensional sample constituted by two coupled elements, namely bone living tissue and artificial bio-resorbable material. We consider this sample as a biphasic mixture porous medium which is subjected to a cyclic bending load slowly variable in time and the mass of which may vary in time; a non-linear theory of porous visco-elastic materials is employed, which is based on classical models [START_REF] Cowin | Linear elastic materials with voids[END_REF]Biot, 1962b,a), to describe the poroelastic material, and the Darcy-Brinkman model (Brinkman, 1949a,b;Biot, 1962a) to simulate the fluid flow through the voids of the mixture. The model proposed by Lekszycki and dell'Isola (2012) is utilized to follow the evolution of the remodeling processes. Herein, we employ a model that should result from a suitable homogenization process. The system under study presents high complexity and therefore needs dedicated tools especially if one is interested in getting a homogenized limit through discrete measure functional approximation, as done e.g. in [START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF]; [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF]; [START_REF] Cecchi | Heterogeneous elastic solids: A mixed homogenizationrigidification technique[END_REF].

Modeling

In this work, we consider a solid mixture composed of two phases, living trabecular bone and a graft constituted by artificial bio-resorbable material endowed with a microstructure that takes into account the presence of connected pores (PIT). These last are filled with marrow which very reasonably plays an important role in dissipation phenomena occurred in bone [START_REF] Braidotti | A critical damping approach for assessing the role of marrow fat on the mechanical strength of trabecular bone[END_REF] and here considered as a highly viscous fluid [START_REF] Bryant | Rheology of bovine bone marrow[END_REF][START_REF] Birmingham | Computational modelling of the mechanics of trabecular bone and marrow using fluid structure interaction techniques[END_REF]. Indeed, the graft is used for creating an osteoconductive structural scaffold in bone neoformation after a vascularization process that makes it active. For the particular nature of its use, therefore, a bio-resorbable graft should be structurally similar to bone and such as to avoid the formation of holes that could compromise the structural integrity of the bone during the process of resorption of bio-material and bone and formation of new bone.

Of course new developments in material engineering have to be taken into account when designing newly conceived grafts of the type of those proposed in the present work. Indeed, it is interesting to notice that the possibilities provided by computer-aided manufacturing are widely increasing the range of material properties which are at disposal, and for instance metamaterials [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF] displaying very advantageous properties as strength-to-weight ratio, safe behavior in fracture or smart interactions with piezoelectric micro-actuators (see for instance dell 'Isola et al (2015a,b); [START_REF] Madeo | Wave propagation in pantographic 2D lattices with internal discontinuities[END_REF]; [START_REF] Giorgio | Piezo-electromechanical smart materials with distributed arrays of piezoelectric transducers: current and upcoming applications[END_REF][START_REF] Giorgio | Multimode vibration control using several piezoelectric transducers shunted with a multiterminal network[END_REF]).

The basic kinematic field is the displacement u(X,t) -designated by the components u i -where t is time and X stands for the selected space parametrization, i.e. three material coordinates in a Lagrangian formulation. Then, to include the notion of deformation, the Green-Saint-Venant strain tensor E is employed and its components are expressed as

E ij (X,t) = 1 2 u i,j + u j,i + u i,k u k,j (1)
In addition, following the framework of the mixture theory (see e.g. [START_REF] Valentín | A finite element-based constrained mixture implementation for arterial growth, remodeling, and adaptation: Theory and numerical verification[END_REF]; [START_REF] Placidi | On the role of grain growth, recrystallization and polygonization in a continuum theory for anisotropic ice sheets[END_REF]; Placidi andHutter (2005, 2006); [START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF]; Tomic et al ( 2014)), we introduce r b , r m and r f as apparent mass densities of bone, bio-material and marrow respectively, and define the mixture mass density by:

r = Â i=b,m,f r i = Â i=b,m,f ri V i = rb V b + rm V m + rf (1 V b V m ) (2)
in which V b and V m are the volume fraction respectively of bone and bio-material, while the Lagrangian porosity is

f = 1 (V b + V m ) (3)
and ri with {i = b,m,f} are the true mass densities of the phases.

Herein each material point at the position X in the reference configuration is also endowed with an internal degree of freedom that, in a sense, equips this point with a structure (see e.g. [START_REF] Altenbach | Acceleration waves and ellipticity in thermoelastic micropolar media[END_REF]; [START_REF] Altenbach | Analysis of the viscoelastic behavior of plates made of functionally graded materials[END_REF]; [START_REF] Federico | Possible approaches in modelling rearrangement in a microstructured material[END_REF]). Specifically, the geometric object representing the microstructure, in our context, is a scalar quantity assumed to be the change of the effective volume of the fluid content per unit volume of the body with respect to an equilibrium volume, z -i.e. the change of the Lagrangian porosity:

z (X,t) = f (c(X,t),t) f ⇤ (X,t) (4)
where f ⇤ stands for the porosity in the reference configuration. All quantities in the reference configuration are denoted by means of a superscript *. For nonhomogeneous porosity, the change of the porosity can be expressed as:

z = -• [f (U u)] = -• w (5)
where u is the displacement of the solid matrix and U is the average fluid-displacement vector. This vector is defined in such a way that the volume of fluid displaced through any inner surface S endowed with a normal n is given by

Q = Z S f U • n dS (6)
Therefore, w = f (U u) is the flow of the fluid relative to the solid but measured in terms of volume (of fluid displaced) per unit area (Biot, 1962b). For the purposes of this work we postulate that the stored energy E (E, z ; r ⇤ b , r ⇤ m ) is a second order polynomial function in the finite strain measure E and the fluid volume distortion z from the reference configuration, kinematic quantities that are objective1 , namely

E = 1 2 C ijhk (r ⇤ b , r ⇤ m )E ij E hk + 1 2 K 1 (r ⇤ b , r ⇤ m )z 2 + 1 2 K 2 z ,i z ,i K 3 (r ⇤ b , r ⇤ m ) z E ii (7)
where, assuming the hypothesis of isotropic material, the Piola stress tensor can be expressed in terms of Young's modulus and Poisson's ratio as

T e ij = C ijhk (r ⇤ b , r ⇤ m )E ij = 2 Y (r ⇤ b ,r ⇤ m ) 2(1+n) E ij + Y (r ⇤ b ,r ⇤ m )n (1 2n)(1+n) E kk d ij (8)
It is worth noting that the energy also depends on the parameters r ⇤ b , r ⇤ m evaluated in a reference -i.e. free-stress-configuration which vary with the time in accordance with evolution rules driven by a mechanical stimulus. The stored energy density ( 7) is non-negative if the following inequalities are satisfied [START_REF] Cowin | Linear elastic materials with voids[END_REF]:

Y > 0, 1 > n > 0.5, K 1 > 0, K 2 > 0, Y 3(1 2n) K 1 > 4K 2 3 (9)
Notice that the term of the energy in Eq. ( 7) depending on the gradient of z allows us to consider typical features of second gradient fluids as e.g. capillarity (see for more details [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF]; [START_REF] Seppecher | Second-gradient theory: application to Cahn-Hilliard fluids[END_REF]; [START_REF] Rosi | Propagation of linear compression waves through plane interfacial layers and mass adsorption in second gradient fluids[END_REF]). For a different view on second-gradient fluid that are known for exhibiting special concentrated interaction we refer to [START_REF] Giusteri | The multiple nature of concentrated interactions in second-gradient dissipative liquids[END_REF].

Material properties of the mixture are assumed inhomogeneous. The Young modulus Y of the mixture is assumed obeying the power-law model:

Y = Y Max b (V b ⇤ ) b b +Y Max m (V m ⇤ ) b m (10)
and varies with reference mass densities [START_REF] Currey | The effect of porosity and mineral content on the young's modulus of elasticity of compact bone[END_REF]Lekszycki and dell'Isola, 2012).

The material constants Y Max b and Y Max m are the maximal elastic moduli of the constituents of the mixture, and the parameters b b and b m are constants that specify the shape of the stiffness. As a first approximation, we assume that Poisson's ratio is constant. The material parameter K 1 can be interpreted as a compressibility coefficient related to the fluid phase inside the pores, and then evaluated by:

K 1 = ✓ f ⇤ K f + (a B f ⇤ )(1 a B ) K dr ◆ 1 (11)
where K f is the fluid bulk modulus, K dr = Y /(3(1 2n)) is the drained bulk modulus of the porous matrix and a B is the Biot-Willis coefficient that lie within the range f ⇤ 6 a B 6 1 (see e.g. [START_REF] Müller | Solid-phase bulk modulus and microinhomogeneity parameter from quasistatic compression experiments[END_REF]). The material parameter K 3 is related to the interaction between the solid and the microstructure which herein is assumed to be spherical. It can be easily demonstrated that an equivalent Lamé parameter l eq = l l v can be introduced to consider the presence of pores [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF].

We hypothesize that the correction due to the presence of pores l v is a fraction of l depending on the reference porosity and as a result K 3 is expressed by

l v = K 2 3 /K 1 = ĝ(f ⇤ )l ) l eq = (1 ĝ(f ⇤ ))l ) K 3 = p ĝ(f ⇤ ) l K 1 (12)
where the function ĝ(f ⇤ ) is non-negative and less than one. The function ĝ(f ⇤ ) maximizes the effect of pores when the level of porosity is high and nullifies this effect when porosity tends to zero; an example of such a function is:

ĝ(f ⇤ ) = A k 3 p ⇢ atan  s k 3 ✓ f ⇤ 1 2 ◆ + atan ⇣ s k 3 2 ⌘ ( 13
)
in which A k 3 2 (0, 1] and s k 3 are coefficients that define the basic shape of the ĝ function and therefore can be chosen to characterize the coupling law used. The coefficient K 2 related to the gradient porosity term, as a first approximation, is assumed to be constant.

Dissipation Sources

Since bone is a hierarchical composite characterized by multiple structures at different length scales, dissipation can arise from a variety of mechanisms. In this work dissipative phenomena occurring in bones reconstructed with bio-resorbable materials are assumed to belong to three main groups: i) Darcy-Brinkman dissipation resulting from fluid flow inside the pores of the mixture (PIT); ii) a viscosity related to the porous mixture under consideration (PLC/PCA); iii) a microstructure viscosity related to the marrow, considered as a very viscous fluid.

Darcy-Brinkman dissipation

In order to examine the mechanics of flow through the porous mixture, it is possible to write a dissipation function, D DB , which depend on the rate of flow of the fluid, i.e. the time derivative of the volume flow vector ∂ w/∂t -a Darcy's dissipation due to the drag at the pores-and on its gradient related to Brinkman dissipation due mainly to shear in the fluid bulk (see Brinkman (1949a,b) and more recently dell'Isola et al ( 2009); [START_REF] Srinivasan | A thermodynamic basis for the derivation of the Darcy, Forchheimer and Brinkman models for flows through porous media and their generalizations[END_REF]). In this case, the dissipation function can be given by

2D DB ( ẇ, -ẇ) = K D ẇ • ẇ + K B -ẇ : -ẇ (14)
where K D and K B are the Darcy permeability second order tensor and the fourth order Brinkmann tensor, respectively. For the particular case of an isotropic medium, this two tensors are reduced to two positive scalar constants. In particular, we can take K D = h/k, where h is the dynamic viscosity of the fluid and k is the coefficient of permeability of the porous mixture (Biot, 1962a).

It is worth noting that the related viscous force can be evaluated by means of the virtual work done by this kind of dissipation on any sub-body C of the whole body B:

d W diss DB = Z C (K D ẇk d w k + K B ẇk,i d w k,i ) dV (15)
Indeed, the viscous force exerted by the 'fluid' marrow on the inner porous surface comes from an integration by parts of the Eq. ( 15) and can be recognized as the dual quantity of the relative flow of the fluid w, which yields the Brinkmann force density:

F Diss DB = K D ẇ -• (K B -ẇ) (16)
This equation has the advantage of approximating the Darcy action (proportional to K D ) for low values of permeability (k) while the other term becomes progressively important when the permeability increases (Brinkman, 1949a) as in PIT. Therefore, in the case under study the generalization due to Brinkman is not negligible.

Mixture viscosity

In [START_REF] Garner | Viscoelastic dissipation in compact bone: implications for stress-induced fluid flow in bone[END_REF] a damping increase with frequency in both wet and dry bone was observed; as a result, it is reasonable to attribute the damping not only to the fluid flow but also to other sources. Indeed, the solid part of the bone is made mainly of an inorganic mineral in the form of small crystals of hydroxyapatite and collagen, that is a natural fibrous polymer. Thus, due to its composition, the bone of a trabecula presents dissipation related to phenomena similar to those which have been already observed to occur in polymers [START_REF] Garner | Viscoelastic dissipation in compact bone: implications for stress-induced fluid flow in bone[END_REF], and to an interstitial fluid flow in the bone canaliculi. Therefore, for the sake of simplicity, we consider a standard way to describe viscoelastic materials, the Kelvin-Voigt model, in which the elastic term is given by the Eq. ( 8), while the viscous part highlighting the shear contribution and the volumetric one, is expressed by

T v ij ( Ė) = 2µ v ✓ Ėij 1 3 Ėkk d ij ◆ + k v Ėkk d ij ( 17
)
where µ v is the shear viscosity and k v is bulk viscosity evaluated for the mixture of trabecular bone and bio-resorbable material (PLC/PCA). Consequently the Rayleigh dissipation function in the solid matrix is

2D s = T v ij Ėij (18)
It should be noted that dissipation in bone can also arise from the frictional sliding of the faces of the micro-cracks present in the solid matrix (see e.g. [START_REF] Scerrato | A simple non-linear model for internal friction in modified concrete[END_REF][START_REF] Scerrato | A micro-structural model for dissipation phenomena in the concrete[END_REF] where internal mechanisms of dissipation in solids are investigated when no relevant damage phenomena occur at the macroscopic level). In this paper, however we do not take account of this possible source of dissipation because the evolution rules of the process of adaptive bone remodeling seem to be related to a viscous effect (in [START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF], on the basis of experimental observations, evolution rules are assumed to be rate dependent) rather than to a Coulomb-like dissipation that is rate independent.

Microstructure viscosity Also the microdeformation related to the variable z can be associated with dissipative mechanisms (see e.g. [START_REF] Cowin | Linear elastic materials with voids[END_REF]). Thus, by motivating our hypothesis by the pragmatic reasoning of simplicity, we can assume the dissipative potential as follows:

2D z ( ż ) = K z ż 2 (19)
where K z is a positive damping coefficient. Given the complex nature of the system considered, more general models of dissipation can be employed to better describe experimental observations, in particular in relation to the dynamic behavior of bone resulting from the application of time-varying loads, as studied in [START_REF] Carcaterra | Dissipation in a finite-size bath[END_REF]; [START_REF] Carcaterra | Fractional dissipation generated by hidden wavefields[END_REF].

Mechanical governing equations

In order to describe the mechanical behavior of the body under study, we employ the Generalized Principle of Virtual Work [START_REF] Green | Multipolar continuum mechanics: functional theory I[END_REF]; [START_REF] Mindlin | On the equations of elastic materials with micro-structure[END_REF]; [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF] and more recently [START_REF] Lekszycki | Application of variational methods in analysis and synthesis of viscoelastic continuous systems[END_REF]; dell'Isola and Placidi ( 2012)) according to which the total virtual work done by the actions applied inside the body d W ins and those which are applied on its boundary d W boun during any arbitrary virtual displacement u i and virtual change of porosity z equals the virtual dissipated energy d W diss . Mathematically:

d W ins + d W boun = d W diss (20)
As a consequence of the assumption (7) on the stored energy and neglecting the inertial effect, the work due to inside interactions have the following form:

d W ins = Z B h T e ij d E ij + K 1 z d z + K 2 z ,i d z ,i K 3 (E ii d z + z d E ii ) b i d u i i dV , (21)
the work done by the contact interactions which can be exerted at the boundary of the considered body can be represented as:

d W boun = Z ∂ t B t i d u i dS + Z ∂ B X dz dS , ( 22
)
where t i is the surface traction on the boundary ∂ t B, and X is a microstructural action introduced to describe the local dilatant behavior of a porous material induced by pore opening, elastic and capillary interaction phenomena among neighboring pores. Interface conditions can be imposed at the boundary between trabecular bone and bio-material zones ∂ B int (see Fig. 1), by assuming the energy density [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF]:

E int = 1 2 K 4 z + z 2 = 1 2 K 4 [[z ]] 2 K 4 > 0 (23)
Where appropriate, similar interactions can be also used for the displacement field to simulate a not perfect gluing, instead of the continuity conditions used here for that field. As a result, the virtual boundary work ( 22) becomes:

d W boun = Z ∂ t B t i d u i dS + Z ∂ B X dz dS Z ∂ B int K 4 [[z ]] d [[z ]]dS (24)
The total virtual dissipated energy can be expressed as

d W diss = Z B h T v ij d E ij + K z ż d z + K D ẇk d w k + K B ẇk,i d w k,i i dV (25)
In Eq. ( 20) an additional kinematic descriptor, i.e. the flow w of the fluid relative to the solid was introduced to describe the state of fluid-solid interaction. Albeit it is an unknown of the problem, this variable is not an independent one as it is clear from Eq. ( 5) and therefore some further discussions are needed to solve our mechanical problem (20). Specifically, we assume this vector w as the gradient of a scalar function, y. In this case the newly introduced potential (y) satisfies Poisson's equation:

z = -• w = -• (-y) (26)
that should be solved together with the Eq. ( 20) and Neumann boundary condition for y, i.e. boundary condition for the relative fluid flow of the kind:

w • n = a(X,t) on ∂ B ( 27
)
where n denotes the unit normal vector to the boundary ∂ B and therefore a is the normal component of the w imposed on the same boundary. With this assumption the virtual work done by the Darcy-Brinkman dissipative action becomes:

d W diss DB = Z B [K D ẇk d w k + K B ẇk,i d w k,i ] dV = = Z B [K D ẏ,k d (y ,k ) + K B ẏ,ki d (y ,ki )] dV (28)
It should be noted that the gradient of ẇ with the assumption of potential flow and, by virtue of Schwarz's theorem is a symmetric second order tensor; as a result, the dissipative function ( 14), and in particular the second term, is correctly-set, that is any relative rigid motion does not entail dissipation. Without this assumption, a more convenient dissipative function should be written in term of the symmetric part of such a gradient Sym(-ẇ) as employed in [START_REF] Srinivasan | A thermodynamic basis for the derivation of the Darcy, Forchheimer and Brinkman models for flows through porous media and their generalizations[END_REF].

As a note, an investigation of stability and bifurcation problems is also important for a possible loss of stability during growth or resorption process, and deserves attention for future investigations (see e.g. [START_REF] Yeremeyev | The stability of the equilibrium of twophase elastic solids[END_REF]; [START_REF] Rizzi | Initial postbuckling behavior of thin-walled frames under mode interaction[END_REF]; [START_REF] Ruta | A beam model for the flexural-torsional buckling of thin-walled members with some applications[END_REF]; Di Egidio et al ( 2007)).

Evolution equations for the mass densities of the two solid phases

According to the current understanding of the biological phenomena involved in osteogenesis, bone cellular populations are organized in basic multicellular units, i.e. temporary structures that operate as smart remodeling actuators [START_REF] Van Bezooijen | Control of bone formation by osteocytes? lessons from the rare skeletal disorders sclerosteosis and van Buchem disease[END_REF][START_REF] Turner | Homeostatic control of bone structure: an application of feedback theory[END_REF]. Osteocytes (i.e. sensor cells) transmit to the effector cells (osteoclasts, which remove bone from unloaded sites and osteoblasts, which add bone to highly loaded regions) a 'mechanical stimulus' that is transduced by forming a dense network throughout the whole bone tissue [START_REF] Klein-Nulend | Mechanical loading and how it affects bone cells: the role of the osteocyte cytoskeleton in maintaining our skeleton[END_REF][START_REF] Baïotto | Theoretical and numerical study of a bone remodeling model: the effect of osteocyte cells distribution[END_REF]. In this way, effector cells are able to optimize the bone mass, shape and architecture according to the loading conditions in a process called bone adaptation [START_REF] Lekszycki | Modelling of bone adaptation based on an optimal response hypothesis[END_REF]. Evolution has created bones that are as light as possible to accomplish their functions of support and protection (see [START_REF] Ambrosi | Perspectives on biological growth and remodeling[END_REF] for a more general reference on biological growth and remodeling). In order to make effective use of amount of material, bones adapt their shape and architecture to the loads applied. Indeed, a balance between bone resorption and bone formation determines the current distribution of bone mass density.

Below the basic assumptions made in order to formulate a proper description of the remodeling process in the considered mixture (constituted by bone tissue and bio-resorbable material) are listed (Lekszycki and dell'Isola, 2012;Andreaus et al, 2014b):

sensor cells are situated only inside the living bone tissue and not in the artificial graft; -the number of sensor cells present in a given material particle of the considered mixture is proportional to the apparent density of the bone tissue r ⇤ b evaluated for the free-stress configuration;

since osteocytes constitute an interconnected network of sensor cells, they can modulate the signal sent to actor cells in order to compensate any loss of signal due for instance to their lower density;

the signal sent from osteocytes to the actor cells decays with increasing source distance, such as an exponential function;

actor cells elaborate and integrate all signals sent by surrounding sensor cells which reach them;

osteoblasts and osteoclasts activities are proportional to the stimulus from the osteocytes;

the number of the actor cells in a given place depends on the porosity of the composite material;

the bone tissue can be resorbed or synthesized but the graft material can be only resorbed.

The differential equations describing the processes of bone and graft material remodeling involve the apparent mass density in the reference configuration that characterizes the evolution of both constituents of the considered mixture as mentioned above. The evolution equations are of the first order and are given by

( ṙ⇤ b = A b (S) H (f ) with 0 < r ⇤ b 6 rb ṙ⇤ m = A m (S) H (f ) with 0 < r ⇤ m 6 r 0 m ( 29
)
where rb and r 0 m are the true mass density of the bone and the initial mass density of the bio-material, respectively. The two constituents in the mixture are different and thus their resorption and synthesis rates will be different. Specifically, the dependence on the stimulus, S, is assumed piece-wise linear as follows:

A b (S) = ⇢ s b S for S 0 r b S for S < 0 (30) A m (S) = ⇢ 0 for S 0 r m S for S < 0 (31)
where s b , r b and r m are positive constants and they represent synthesis and resorption rate of the bone and the resorption rate of the graft material, respectively.

The function H depends on porosity and is evaluated on the basis of a normalized specific surface (S s ) i.e. the inner surface area of the pores per unit volume of bone [START_REF] Martin | Porosity and specific surface of bone[END_REF] (see next section and Fig. 2(b)). This function is conceived in order to account for the influence of 'effective' porosity (i.e. effective for the deposit of actor cells) on the biological activity of actor cells: when effective porosity is too large there is not enough material on which actor cells may deposit; when it is too small, there is not enough free space in the pores to allow their mobility and deposit (Lekszycki and dell'Isola, 2012). For the above reasons, and in particular since the actor cells operate on the internal surface of the pores, the function H was related to the aforesaid normalized specific surface, as it will be seen in more detail in the following section.

In the present context the explicit expression of the stimulus is (a) Comparison among different elementary cell geometries and some reference curves [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF][START_REF] Martin | Porosity and specific surface of bone[END_REF] porosity ? where P s ref and P r ref are two activation thresholds, the first for synthesis and the other for resorption, that define the width of the 'lazy zone' associated with a biological equilibrium state while the signal generated by osteocytes can be expressed in the non-local integral form as

S(X,t) = 8 > < > : P(X,
P (X,t) = R B E s (X 0 ,t) v ⇥ r ⇤ b (X 0 ,t) ⇤ e kX X 0 k 2 2D 2 dX 0 R B e kX X 0 k 2 2D 2 dX 0 (33)
in which the characteristic length D represents the range of action of sensor cells and the signal coming from osteocytes is assumed depending on the density of strain energy of solid matrix 2 E s = (CE) : E in the region where they are located. Since sensor cells are interconnected through dendritic processes inside the canalicular network, we assume that they can exchange information each other through these processes on the current state of the stimulus. To mimic this non-local behavior, we employed a spatially averaged stimulus over a spherical 'zone of influence' as in [START_REF] Kumar | Dissipation energy as a stimulus for cortical bone adaptation[END_REF]. The function v can be interpreted, in the framework of the feedback control theory [START_REF] Turner | Homeostatic control of bone structure: an application of feedback theory[END_REF], as a proper gain for the strain energy that takes the real activities of the osteocytes into account, and therefore we set it as

v (r ⇤ b ) = h tanh(x r ⇤ b ) 0 < h 6 1 (34)
The key idea here is to reach a saturation when a critical number of sensor cells, tuned with the positive parameter x , is exceeded in order to have the best gain possible when there are enough cells and then a good estimate of the strain energy. On the contrary, if there are only few sensor cells, they compensate their insufficient number with a greater activity compatibly with the hypothesis that the number of sensor cells is proportional to the apparent mass density of the bone tissue.

It is a matter of common knowledge that the outside layer of bones consists of a compact tissue, called cortical bone, while the inner core is constituted by a more porous tissue, called cancellous, spongy or trabecular bone. The overall porosity of the bone increases when approaching its inner part. Keeping in mind this observation, the expression of the stimulus proposed in the present paper, i.e. the Eq. ( 33), seems to better fit the distribution of porosity experimentally observed than the stimulus employed in Lekszycki and dell'Isola (2012); Andreaus et al ( 2015); [START_REF] Giorgio | The influence of different loads on the remodeling process of a bone and bioresorbable material mixture with voids[END_REF] and shown below for comparison purposes.

P (X,t) = Z B E s (X 0 ,t) v [r ⇤ b (X 0 ,t)] e kX X 0 k 2 2D 2 dX 0 ( 35
)
In Fig. 3 indeed the distributions over the sample of the two stimuli considered, Eqs. ( 33) and ( 35), in a non-dimensional form (see section 2.6) and in the case of a pure bending deformation as sketched in Fig. 1 are compared at the initial stage. For the stimulus (35) a border effect resulting in a lower production of bone right on the borders parallel to the x-axis in the bone zone can be noticed in contrast to what happens for the stimulus (33).

Finally, we remark that because of the complexity of the model, a large number of constitutive parameters should be taken into account. Therefore, a specific procedure is needed to identify all these parameters. In that regard, the method outlined in [START_REF] Placidi | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF]; [START_REF] Dietrich | Problems of identification of mechanical characteristics of viscoelastic composites[END_REF]; [START_REF] Lekszycki | Modelling and identification of viscoelastic properties of vibrating sandwich beams[END_REF] can be profitably employed.

Evaluation of the normalized specific surface

To generate the curves shown in Fig. 2(a), we proceeded as sketched in the following way. First, some geometric shapes have been appropriately selected which were particularly elemental and which, periodically repeated in space, enabled to fill it completely and reproducing in a schematic way the arrangement of trabeculae (see [START_REF] Wang | Scale and boundary conditions effects on the apparent elastic moduli of trabecular bone modeled as a periodic cellular solid[END_REF] for a similar approach). Indeed, we interpret the bone as a cellular solid, i.e. a particular solid made up of an assembly of 'cells' packed together so that they fill space, each of them is constituted by an interconnected network of solid struts or plates which form the edge or the faces of the cells [START_REF] Gibson | Cellular solids: structure and properties[END_REF]. In a three dimensional environment a great variety of cell shapes is possible, which can be packed together to fill space. The considered forms (see Fig. 4 and [START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF]; [START_REF] Wang | Scale and boundary conditions effects on the apparent elastic moduli of trabecular bone modeled as a periodic cellular solid[END_REF]) are: a) tetrahedron, b) cube, c) honeycomb, d) truncated octahedron. The choice of the cellular shape is important because the properties of this kind of solid depend directly on the shape and structure of the cells. Therefore, it is necessary to characterize their size, shape and topology.

It is noted that each module is delimited by cylinders which represent a simplified shape for trabeculae. By using the same radius for all cylinders of the module, the surface area, S i , of the cylinder face of the inner part of the module was calculated. Then, it was calculated the volume of the void, V v , obtained by subtracting the volume of the full cylinders to the volume of the whole module, V m (e.g. the volume of the cube in the case b) so as to evaluate the porosity -this procedure takes into account only the inter-trabecular space porosity (PIT)-as the ratio f IT = V v /V m . By varying the radius, the curve which has the porosity f IT in abscissas and the specific surface in ordinates was built point by point. At the end, the specific surface S i was normalizing by dividing it by the maximum value of the S i itself, obtaining in such a way the normalized specific surface S s . Figure 2(a) shows the curves S s related to the various modules (see Fig. 4), in addition to the parabola [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF] and to the polynomial interpolation [START_REF] Martin | Porosity and specific surface of bone[END_REF].

As it can be seen from Fig. 2(a) the S s curves are concentrated in a rather narrow band; in particular the curves related to the cubic and the tetrahedral modules coincide and therefore this has been chosen as the reference curve in the numerical simulation which follows.

The need to introduce a not-effective porosity in close proximity of both the total filling level and the total vacuum, as already observed in [START_REF] Martin | Porosity and specific surface of bone[END_REF], has led to the curve of Fig. 2(b) which was obtained by an appropriate contraction of the original curve along the f -axis. Specifically, the porosity can be decomposed into two contributions, namely f IT related to the inter-trabecular space and f R related to the vascular, lacunar-canalicular and collagen-apatite space:

f = f IT + f R (36)
In a first approximation, the not-effective porosity f R can be assumed as a linear function of being n 0 = 0.166 the porosity related to the cortical bone [START_REF] Martin | Porosity and specific surface of bone[END_REF]; indeed when f IT is zero there is no inter-trabecular space but only the void space present in the cortical bone.

f IT : f R = n 0 f IT + n 0 (37)
On the other hand, it is reasonable that decreasing the trabecular solid, the not-effective porosity f R decreases. With these assumptions and defining from Eqs. ( 36) and (37):

f = f (f IT ) = (1 n 0 )f IT + n 0
we can define the function H as:

S s (f IT ) = S s [ f 1 (f )] := H(f ) (38)
In Tab. 1, we set the level of H, near to zero porosity in the range [0, n 0 ], a little larger than zero to avoid that the material will not have a chance in future to evolve in this domain and change its porosity any more, what contrasts with a real situation. Indeed, we adopt this low threshold not only to allow a future evolution, as done in Lekszycki and dell'Isola (2012), but also to penalize the remodeling process where no effective porosity is available; on the other hand, we are not interested in considering the effects that happen at cellular level, and lower.

Stimulus driven by Strain and Dissipated Energy

Laboratory experiments [START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF] have confirmed that both loading frequency and strain rate are important determinants of bone adaptation. These results have been summarized in the following succinct statements [START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF]: (a) dynamic strains drive bone adaptation; (b) the strain stimulus is increased if the magnitude or frequency of the dynamic signal is increased; and (c) increasing strain rate enhances the strain stimulus. To incorporate these facts into a mathematical formula, [START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF] considered that peak strain rate is proportional to the frequency of the loading waveform and strain magnitude and assumed that the strain stimulus is proportional to strain rate. From these considerations it is clear that the strain rate, and then the dissipated energy, which depends on it, play the role of ideal candidates to represent an important component of the mechanical stimulus. Already, other researchers [START_REF] Kumar | Dissipation energy as a stimulus for cortical bone adaptation[END_REF][START_REF] Kumar | Modeling of cortical bone adaptation in a rat ulna: Effect of frequency[END_REF][START_REF] Pereira | The influence of load repetition in bone mechanotransduction using poroelastic finite-element models: the impact of permeability[END_REF] have studied both idealized and realistic situations, taking as stimulus just the density of energy dissipation. They assumed that the load induced fluid flow can occur at the lacuno-canalicular porosities, and generated a 3D finite-element model of the cortical bone as a continuous, homogenized linear poroelastic material. The intramedullary canal, which is at a higher scale, is treated as a source for fluid to move freely into and out of the cortical bone.

In the present work, in order to take all these features into account, a further assumption was added to the above mentioned ones: the stimulus is the linear combination of the straindependent energy, E s , and the dissipation power, D s (with coefficients a and b , respectively).

Particularly

P (X,t) = R B (a E s + b D s ) v r ⇤ b e kX X 0 k 2 2D 2 dX 0 R B e kX X 0 k 2 2D 2 dX 0 (39)
The advantage to propose such an expression for the stimulus is twofold. On the one hand, it allows to simultaneously address both static and dynamic situations; at the two extremes there are indeed instantaneous impacts, when the bone is called to perform the function of shock absorber, and permanent static loads, where carrying capacity of bone ensures the homeostatic equilibrium. A further advantage of using the expression (39) is to guarantee a greater generality to processing and the possibility of optimizing the coefficients (weights) of the linear combination in order to approximate to the best experimental results. This opens the way for the optimal adaptive control, particularly effective in the preparation of programs of even extreme training, and of post-traumatic and/or post-operative rehabilitation.

Besides, this kind of stimulus involving the actual spatially averaged values of the elastic strain and its rate is close related to a Proportional-Derivative (PD) control rule as employed in [START_REF] Andreaus | Optimal-tuning PID control of adaptive materials for structural efficiency[END_REF][START_REF] Andreaus | An optimal control procedure for bone adaptation under mechanical stimulus[END_REF]Andreaus et al ( , 2013Andreaus et al ( , 2014a) ) where the actual mechanical stimulus is defined by comparing the Strain Energy Density (SED) with a reference Strain Energy Density (SED*).

Non-dimensional form

In numerical simulations, a non-dimensional form of the considered problem is achieved normalizing the variables involved in the analysis by reference quantities. In this regard, quantities in non-dimensional form are denoted with a superimposed tilde. In particular, the stored energy density E , Eq. ( 7), is normalized with respect to the maximal bone stiffness (40)

As a consequence of this normalization, the material parameters involved in Eq. ( 7) become in a dimensionless form:

Ỹ = (V ⇤ b ) b b + Ỹ Max m (V ⇤ m ) b m , (41)
= K f /Y Max b , K2 = K 2 /(Y Max b L 2 0 ), l = l /Y Max b , K4 = K 4 /(Y Max b L 0 ), (42)
with L 0 a characteristic length. Herein L 0 is assumed to be equal to the length of the sample.

KD = K D L 2 0 /(t ref Y Max b ), KB = K B /(t ref Y Max b ), Kz = K z /(t ref Y Max b ), μv = µ v /(t ref Y Max b ), kv = k v /(t ref Y Max b ). (44)
The external applied force t i is normalized following the same criterion with respect to the stiffness Y Max b . Non-dimensional mass densities of bone and material are normalized with respect to rb = 1800 kg/m 3 and can be defined as:

mass density of bone, rb = r b / rb = V b , -mass density of material, rm = r m / rb = ( rm / rb ) V m , and therefore the non-dimensional total mass density is r = rb + rm .

The non-dimensional stimulus can be defined as: 33) and Eq. ( 39). As a result, from Eqs. (29) the non-dimension form of the evolution equations for the mass densities is:

S = P Pi ref with i = {r, s}, (45)
8 > > < > > : ∂ r⇤ b ∂ t = Ãb ( S)H(f ) ∂ r⇤ m ∂ t = Ãm ( S)H(f ) , (46)
where with some algebra from Eqs. ( 30) and (31) it follows:

Ãb S = ⇢ sb S for S > 0 rb S for S < 0 , (47)
; rb = Y Max b t ref rb r b ; rm = Y Max b t ref rb r m . ( 49
)
3 Numerical simulations and discussion of the results

This section presents the results obtained through a numerical investigation conducted by the FEM code COMSOL Multiphysics R . As for code description, the evolutionary problem in our case study involves time-dependent differential equations of first order for mass densities of bone tissue and bio-material coupled with PDE governing mechanical equilibrium. The solver employs a discrete time step increment algorithm based on the backward differentiation formula method and uses a Newton-type iterative method to solve this nonlinear systems of PDE's. Since the integro-differential equations we are integrating are strongly non-linear, a very delicate issue in the performed numerical simulations regards the choice of time step that in our case is variable and optimized. The current numerical study is performed by means of a commercial FEM code as mentioned above, however, the objective of achieving greater accuracy in the solution of a problem as complex as that addressed herein might justify the adoption of more sophisticated and original techniques based on finite elements such as those proposed in Cazzani et al (2014b); Turco and Aristodemo (1998); Cazzani et al (2014a[START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF]; Greco andCuomo (2014, 2016); [START_REF] Solari | A numerical algorithm for the aerodynamic identification of structures[END_REF]; [START_REF] Cuomo | A variational model based on isogeometric interpolation for the analysis of cracked bodies[END_REF]. The object of the study is a rectangular sample of bicomponent material initially characterized by two distinct areas: the host bone, and the biomaterial graft (see Fig. 1). The mass density distributions of bone and of graft material are assumed to be at the initial stage constant, with the same value for the two phases ( r 0 b = r 0 m = 0.5). One edge of the specimen is assumed fixed in the longitudinal direction, while on the opposite edge is applied the bending traction that has a fixed component and a variable one in time with sinusoidal law:

f x (y,t) = ✓ 2y W 1 ◆ [ f 0 + f 1 sin(Wt)] (50)
where f 1 = f 0 /2; according to the non-dimentional form introduced in the Subsect. 2.6, the following values have been given to the parameters of Eq. ( 50 Figure 5 shows the time evolutions of bone mass density at the probe point P m1 (Fig. 1). In particular, Fig. 5(a) shows the curves obtained by using as a stimulus that provided by Eq. ( 33), i.e. the energy density of elastic deformation. It is noted the difference between static ( f 0 6 = 0, f 1 = 0) and dynamic ( f 0 6 = 0, f 1 6 = 0) application, Eq. ( 50), of the load and the substantial independence of the frequency on the response. In Fig. 5(b) the adoption of the mixed stimulus given by Eq. ( 39) allows to highlight the dependence of the response from the frequency of the applied load, as experimentally observed. In particular, it is noted that the growth process is favored by the increase in frequency, even if the level of saturation of the mass density is reduced -albeit slightly-with increasing frequency, because of the interaction with the bioresorbable material; in fact, the accelerated growth of the bone does not allow time to the biomaterial to be reabsorbed and then does not leave space to the bone to fill it.

Figures 6 and7 show the distributions of the mass density of bone (Fig. 6) and biomaterial (Fig. 7) over the sample at the end of the process for the different frequencies. The two initially distinct areas of bone and artificial material are highlighted by the labels 'B' and 'M'. The stimulus combination of strain energy and power dissipation, Eq. ( 39), was used in this analysis. With increasing frequency, there is a diversification of the density distribution of bone mass, which reaches saturation in areas increasingly wide near the edges parallel to the x-axis in the region of the bone (i.e. in the most stressed zones). In the region of the bio-material there is a tendency to greater uniformity in the distribution of mass density, with a groove located at the neutral axis, always less accentuated as the frequency increases. In general, the increase in the frequency favors a greater uniformity. Figure 7 shows in a complementary manner the same effects described above seen from the point of view of the bio-material, that is a lower resorption of the material occurs where bone growth is faster; this phenomenon is the more marked the higher the frequency. At low frequencies and in the static case, the bone grows less quickly, and the bio-material has more time to be reabsorbed.

Figure 8 shows the distribution of the porosity of the mixture at the end of the process at the different frequencies examined. We can observe that the overall porosity decreases with increasing frequency, and in particular, this decrease is more evident in the region of the bio-material and in areas close to the edges parallel to the x-axis. Comparing Fig. 8 with the corresponding Figs. 6 and 7 it can be seen that the low porosity in the zones of the neutral axis is mainly due to the fact that, as the frequency increases, resorption of both bone and the material is less significant.

Conclusions

The subject of this work is the study of biological and mechanical interaction of a graft of bio-resorbable material with a portion of bone tissue under the action of a bending load; in this regard some aspects of the phenomena related to the remodeling under mechanical stimulus have been described. The two-component system is considered as a porous mixture. In particular, the behavior of the material is described according to the model used for the poro-visco-elastic materials with fully saturated voids. The fluid flow in inter trabecular porosity (PIT) is evaluated by means of a Darcy-Brinkman dissipative model. At the level of Collagen-Apatite porosity (PCA) and lacunar-canalicular porosity (PLC) the dissipation is taken into account by means of the Kelvin-Voigt rheological model. At the level of porosity of PIT, it is also considered a dissipation associated with the micromorphic variable that is the change of porosity. A non-linear behavior has been assumed to describe the relationship between strain and displacement. As a stimulus, a linear combination of strain energy density and power dissipation, duly weighted, has been proposed. The key idea of considering such a form of the stimulus has been suggested by the knowledge that the bone is called to play a role in carrying capacity of permanent static loads, fatigue resistance against excitation varying in time and dynamic resilience against traumatic-impulsive events. The non-local character of osteocytes's sensitivity has been schematically accounted for by spatially averaging the stimulus over a suitable zone of influence here supposed to be spherical. The expression used for the stimulus lends itself to the idea of optimizing the values of the coefficients of the linear combination, in order to obtain a numerical simulation of the experimental results as much as possible close to real situations. Moreover, this optimization can lead to interpret the phenomenon of functional adaptation in view of adaptive optimal control, especially in relation to the preparation of programs of extreme training and of posttraumatic and/or post-operative rehabilitation. Finally, the definition of the stimulus adopted in this work allows to put in evidence the dependence of the phenomenon of remodeling on the frequency of the applied load, as already observed in previous works. Another noteworthy result achieved in the present work was to have identified the "normalized specific surface" useful to the characterization of the effect of porosity on remodeling; this identification was made on the basis of a geometric model based on the assumption of some modular forms that simulate the trabecular network. The values obtained are very close to those obtained by similar interpolated curves from experimental results.

Fig. 1

 1 Fig. 1 Sample used in numerical simulations. The labels 'B' and 'M' stand for Bone and graft Material, respectively.

Fig. 2

 2 Fig. 2 Normalized Specific Surface

Fig. 3

 3 Fig. 3 Initial normalized stimulus in the case of a pure bending deformation on a specimen as sketched in Fig. 1. The labels 'B' and 'M' stand for Bone and graft Material, respectively.

Fig. 4

 4 Fig. 4 Unit module of idealized periodic model of trabecular bone regarded as cellular solid.

  and therefore takes the form:Ẽ = E /Y Max b .

  from the assumptions Ẽs = E s /Y Max b , Ds = (t ref /Y Max b )D s and b = b /t ref , and from Eq. (

  Fig. 5 Time evolutions of bone mass density in the probe point P m1 . Parametric analysis varying load frequencies: {0, 5, 10, 20} cycles per unit of time (a week).

  Fig. 6 Bone mass density at the end of the process. Parametric analysis varying load frequencies: {0, 5, 10, 20} cycles per unit of time (a week). The labels 'B' and 'M' stand for Bone and graft Material, respectively.

  Fig. 7 Biomaterial mass density at the end of the process. Parametric analysis varying load frequencies: {0, 5, 10, 20} cycles per unit of time (a week). The labels 'B' and 'M' stand for Bone and graft Material, respectively.

Table 1

 1 Values of the normalized Specific surface S s vs f IT and H vs f .

	Effective porosity f IT	Specific surface S s	Porosity f	Function H
	0.0000	0.0000	0.0000	0.0300
	0.0002	0.0068	0.1662	0.0300
	0.0011	0.0218	0.1669	0.0300
	0.0032	0.0459	0.1687	0.0321
	0.0073	0.0799	0.1721	0.0722
	0.0139	0.1247	0.1776	0.1233
	0.0239	0.1823	0.1859	0.1817
	0.0381	0.2564	0.1978	0.2565
	0.0580	0.3591	0.2144	0.3590
	0.0861	0.4952	0.2378	0.4942
	0.1225	0.6148	0.2682	0.6147
	0.1662	0.7188	0.3046	0.7186
	0.2161	0.8059	0.3462	0.8058
	0.2713	0.8773	0.3923	0.8773
	0.3306	0.9324	0.4417	0.9323
	0.3929	0.9710	0.4937	0.9709
	0.4573	0.9935	0.5474	0.9935
	0.5226	1.0000	0.6018	1.0000
	0.5878	0.9901	0.6562	0.9900
	0.6518	0.9641	0.7096	0.9640
	0.7136	0.9217	0.7611	0.9216
	0.7720	0.8632	0.8098	0.8631
	0.8262	0.7887	0.8551	0.7886
	0.8749	0.6977	0.8957	0.6976
	0.9171	0.5904	0.9309	0.5902
	0.9518	0.4673	0.9598	0.4669
	0.9779	0.3277	0.9816	0.3265
	0.9943	0.1719	0.9952	0.1279
	0.9991	0.0706	0.9992	0.0202
	1.0000	0.0000	1.0000	0.0000

  By considering a reference time t ref defined as the time duration in which the physiological processes that lead to fulfillment of the phenomenon of material resorption and bone synthesis are accomplished, the non-dimensional time t = t/t ref is introduced. In this paper the reference time is assumed t ref = 6.048⇥10 5 s in order to have a reference time of a week. Analogously, by defining the non-dimensional form of the total dissipation, D = D DB + D s + D

			z ,
	as follows:
	D = (t ref /Y Max b	)D	(43)
	one gets:

Table 2

 2 Material coefficients in non-dimensional form.

		Ỹ Max b		Ỹ Max m		sb	rb	rm
		1			0.8	7.26⇥10 5	6.05⇥10 5	9.07⇥10 5
		r 0 b			r 0 m		b b = b m	A k3	s k3
		0.5			0.5		2	0.9	15
		D			a		b	Pr ref	Ps ref
		0.1			1		1	4.68⇥10 8	5.72⇥10 8
		Kf			KD		KB	μv	kv
		0.1	1.32⇥10 2	1.32⇥10 4	2.50⇥10 5	2.00⇥10 5
		K2			K4		Kz	x	h
		1.0⇥10 3		0.1	1.67⇥10 6	10	0.2
				Ãm S =	⇢	0 rm S	for S > 0 for S < 0	,	(48)
	with	sb =	Y Max b rb	t ref	s b

i.e. they are invariant under superimposition of a rigid body motion.

Acknowledgements The research was funded by the Italian Ministry of University and Research, under the Scientific Research Program of Relevant National Interest: Year 2010-2011, Protocol 2010MBJK5B-005, Title "Dynamics, Stability and Control of Flexible Structures".