**F FP** : .

The tissue element under grow is submitted to a surface force field *S***f** (surface density) and to line densities *p* , *p* defined as the projections onto the unit vectors , *g g* **τ ν** resp. along the contour of the open growing surface *Sg* (Figure 5); Hence, those line densities are respectively tangential and normal to the surface *Sg* (forces acting in the tangent plane).

The variation of the previously built potential energy of the growing tissue element*V* is next evaluated, assuming applied forces act as dead loads, using the fact that the variation is performed over a changing domain (Petryk and Mroz, 1986), here the growing surface *Sg* . We refer to the recent work in (Ganghoffer, 2010a) giving the detailed calculation of the material forces for surface growth, very similar to present developments.

The variation of the volumetric term (first term on the right hand side of *V* ) can be developed from the equalities (A2.1) through (A2.3) given in **(**Ganghoffer, 2010a, Appendix 2):

$$\delta \left( \int\_{\Omega\_{\mathcal{g}}} \mathcal{W}\_0 \Big( \mathbf{F}, \mathbf{X}\_{\mathcal{g}} \Big) d\mathbf{x}\_{\mathcal{g}} \right) = \int\_{\partial\Omega\_{\mathcal{g}}} \left( \mathbf{E} \, \delta \mathbf{X}\_{\mathcal{g}} + \mathbf{p} \, \delta \mathbf{x} \right) \mathbf{N} d\left( \partial \Omega\_{\mathcal{g}} \right) + \upsilon \, t. \tag{5.2}$$

with volumetric terms denoted as 'v.t.' that will not be expressed here, as we are mostly interested in surface growth. The r.h.s. in previous identity is a pure surface contribution involving the volumetric Eshelby stress built from the volumetric strain energy density and the so-called canonical momentum

$$\mathbf{E} \coloneqq \mathcal{W}\_0 \mathbf{I} - \mathbf{F}^t. \mathbf{p} \quad \mathbf{p} \coloneqq \frac{\partial \mathcal{W}\_0}{\partial \nabla \mathbf{x}} \tag{5.3}$$

As we perform material variations over an assumed fixed actual configuration, the contribution of the canonical momentum vanishes ( **x 0** ). Observe that the volumetric Eshelby stress **Σ** triggers surface growth in the sense of the boundary values taken by the normal Eshelby-like traction . **Σ N** . The variation of the surface energy contribution *<sup>S</sup>* can be expanded using the surface divergence theorem (equality (3.15) in Ganghoffer, 2010a) as

$$\delta \left( \int\_{\mathcal{S}\_{\mathcal{S}}} \nu^S \left( \tilde{\mathbf{F}}, \mathbf{N}; \mathbf{X}\_{\mathcal{S}} \right) d\sigma\_{\mathcal{S}} \right) = \int\_{\mathcal{S}} \left[ \nabla\_{\mathcal{S}} \tilde{\mathbf{Z}} - \mathbf{H} \mathbf{K}^t \, \partial\_N \nu^S + \left( \partial\_{X\_{\mathcal{S}}} \nu^s \right)\_{\text{exp}I} + \tilde{\mathbf{F}}^T \, \mathbf{f}\_{\mathcal{S}} \right] \delta \mathbf{X}\_{\mathcal{S}} d\sigma\_{\mathcal{S}} \tag{5.4}$$

The surface energy momentum tensor (of Eshelby type) is then defined as the second order tensor

$$\tilde{\mathbf{T}} \coloneqq \hat{\boldsymbol{\sigma}}\_{\tilde{F}} \boldsymbol{\mu}^{S} \to \tilde{\boldsymbol{\Sigma}} \coloneqq \tilde{\mathbf{F}}^{T} \, \tilde{\mathbf{T}} - \boldsymbol{\mu}^{s} \mathbf{I}\_{s} \tag{5.5}$$

basing on the *surface stress* **T** . The Lagrangian curvature tensor is defined as : **K N** *<sup>R</sup>* . The chemical potential as the partial derivative of the surface energy density with respect to the superficial concentration

$$\mu\_k \coloneqq \frac{\partial \boldsymbol{\upmu}^S}{\partial \boldsymbol{n}\_k^\sigma} \bigg|\_{\boldsymbol{X}, \boldsymbol{F}, \boldsymbol{N}} \equiv \mu\_k \left( \boldsymbol{n}\_k^\sigma \right) \tag{5.6}$$

The contributions arising from the domain variation due to surface growth are considered as irreversible.

The material surface driving force (for surface growth) triggers the motion of the surface of the growing solid; it is identified from the material variation of *V* as the vector acting on the variation of the surface position

$$\tilde{\mathbf{Y}}\_{\mathcal{S}} \coloneqq \mathbf{Z}\mathbf{N} + \nabla\_{\mathcal{S}}\tilde{\mathbf{Z}} - \mathbf{P}\mathbf{K}^{t}\,\partial\_{N}\boldsymbol{\mu}^{\mathcal{S}} + \mu\_{k}\nabla\_{\mathcal{S}}\mathbf{n}\_{k}^{\sigma} - \mathbf{f}\_{\mathcal{S}}\tag{5.7}$$

itself built from *the surface stress* : *<sup>S</sup> <sup>F</sup>* **<sup>T</sup>** , and on the curvature tensor : **K N** *<sup>R</sup>* in the referential configuration.

#### **5.2 Bone remodeling**

392 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

**X***S* on *Sg* (no tilda notation is adopted here since the support of **X***S* is strictly restricted to

defined as the projections onto the unit vectors , *g g* **τ ν** resp. along the

*<sup>k</sup>* the chemical potential of the surface

maps material lengths (or material

*V* ) can be

can

 

. The surface gradient **F**

material forces for surface growth, very similar to present developments. The variation of the volumetric term (first term on the right hand side of

*g g*

contribution of the canonical momentum vanishes (

 

*g g*

*S S*

tensor

, with

tangent vectors) onto the deformed surface; it is elaborated as the surface projection of **F**

**F FP** : .

The tissue element under grow is submitted to a surface force field *S***f** (surface density) and

contour of the open growing surface *Sg* (Figure 5); Hence, those line densities are respectively tangential and normal to the surface *Sg* (forces acting in the tangent plane). The variation of the previously built potential energy of the growing tissue element*V* is next evaluated, assuming applied forces act as dead loads, using the fact that the variation is performed over a changing domain (Petryk and Mroz, 1986), here the growing surface *Sg* . We refer to the recent work in (Ganghoffer, 2010a) giving the detailed calculation of the

developed from the equalities (A2.1) through (A2.3) given in **(**Ganghoffer, 2010a,

*W dx g g <sup>g</sup> <sup>g</sup>*

with volumetric terms denoted as 'v.t.' that will not be expressed here, as we are mostly interested in surface growth. The r.h.s. in previous identity is a pure surface contribution involving the volumetric Eshelby stress built from the volumetric strain energy density and

<sup>0</sup> : . *W <sup>t</sup>* **Σ I Fp** <sup>0</sup> :

As we perform material variations over an assumed fixed actual configuration, the

Eshelby stress **Σ** triggers surface growth in the sense of the boundary values taken by the normal Eshelby-like traction . **Σ N** . The variation of the surface energy contribution *<sup>S</sup>*

be expanded using the surface divergence theorem (equality (3.15) in Ganghoffer, 2010a) as

, ; . .. . . *<sup>S</sup>*

The surface energy momentum tensor (of Eshelby type) is then defined as the second order

: :. *S Ts* **<sup>T</sup>** *<sup>F</sup>*

exp

*S tS s T*

 

<sup>0</sup> , . .. . .

 

*W x*

*<sup>S</sup> <sup>g</sup> S NX S S <sup>g</sup> <sup>l</sup>*

 

 

**<sup>Σ</sup> FT I***<sup>s</sup>* (5.5)

*d d* 

**FNX Σ Π <sup>K</sup> Ff X** (5.4)

**F X <sup>Σ</sup> X pxN** (5.2)

*d vt*

**<sup>p</sup>** (5.3)

**x 0** ). Observe that the volumetric

 

the surface *Sg* ), and chemical energy *k k n*

concentration of species *nk*

to line densities *p* , *p*

Appendix 2):

(onto the tangent plane to *<sup>a</sup>* ), viz

 

the so-called canonical momentum

Bone is considered as a homogeneous single phase continuum material; from a microstructural viewpoint, bone consists mainly of hydroxyapatite, a type-I collagen, providing the structural rigidity. The collageneous fraction will be discarded, as the mineral carries most of the strain energy (Silva and Ulm, 2002). The ultrastructure may be considered as a continuum, subjected to a portion of its boundary to the chemical activity generated by osteoclasts, generating an overall change of mass of the solid (the mineral fraction) given by

$$\frac{d}{dt} \int\_{\Omega\_{\mathcal{g}}} \rho\_{\mathcal{g}} d\mathbf{x}\_{\mathcal{g}} = \int\_{\mathcal{S}\_{\mathcal{g}}} \rho\_{\mathcal{g}} \mathbf{V}\_{\mathcal{S}} \mathbf{N} d\sigma\_{\mathcal{g}}.$$

The quantity . *<sup>g</sup> <sup>S</sup> <sup>g</sup>* **V N***d* therein represents the molar flux of bone material being dissolved, hence

$$
\rho\_{\mathcal{X}} V\_N d\sigma\_{\mathcal{X}} = M \mathfrak{l} d\sigma\_{\mathcal{X}} \tag{5.8}
$$

with *VN* the normal surface velocity, *M* the bone mineral molar mass, and / *J VM g N* the molar influx of minerals (positive in case of bone apposition, and negative when resorption occurs). Clearly, the previous expression shows that the knowledge of the normal surface growth velocity determines the molar influx of minerals. Estimates of the order of magnitude of the dissolution rate given in (Christoffersen at al., 1997), for a pH of 7.2 (although much higher compared to the pH for which bone resorption takes place) and at a temperature of 310K, are indicative of values of the molar influx in the interval 9 8 12 *<sup>J</sup>* 10 ,1.8.10 . . *mol s m* . The osteoclasts responsible for bone resorption attach to the bone surface, remove the collageneous fraction of the material by transport phenomena, and diffuse within the material. This osteoclasts activity occurs at a typical scale of about 50*m* ,

Thermodynamics of Surface Growth with Application to Bone Remodeling 395

Fig. 6. Modeling occurring during growth of the proximal end of the femur. Frontal section of the original proxima tibia is indicated as the stippled area. The situation after a growth of 21 days is superimposed. Bone formation (+) and bone resorption zones indicated [Weiss,

with **ε P**. .. **ε I***S rr* **ε ε e e** the surface strain (induced by the existing volumetric strain), and *A*,*B* mechanical properties of the surface, expressing versus the surface density of minerals and the maximum value of the traction modulus as (the Poisson ratio is selected as

> 

(5.13)

; this applied stress generates a preexisting homogeneous

*<sup>S</sup> t* and its superficial concentration. We

(5.14)

(5.12)

 

 

3 3

 

*S*

 

stress, consisting of the superposition of an axial and a radial component, as

**σ** 

stress state within the bulk material, inducing a surface stress given by

(5.10). The surface stress results from (5.11), (5.12) as

surface*V t <sup>N</sup>* , the surface density of minerals

in the cylindrical basis **eee** *r z* , ,

max max ; 1 2 1 21 *Et Et S S A B* 

As the surface of bone undergoes resorption, its mechanical properties are continuously changing from the bulk behavior, due to the decrease of mineral density as reflected in

: 2

 *rr r r zz z z* **ee ee** 

**T σ εε I ε**

The unknowns of the remodeling problem are the normal velocity of the bone

shall herewith simulate the resorption of a hollow bone submitted to a composite applied

*mech A Btr <sup>S</sup>*

1988].

0.3 )

which is much larger compared to the characteristic size of the ultrastructure; the resorption phase takes typically 21 days (the complete remodeling cycle lasts 3 months). The osteoclasts, generate an acid environment causing simultaneously the dissolution of the mineral - hydroxyapatite, a strong basic mineral 3 4 2 2 <sup>3</sup> *Ca PO Ca OH* , abbreviated HA in the sequel - and the degradation of the collageneous fraction of the material. The metabolic processes behind bone remodeling are very complicated, with kinetics of various chemical substances, see (Petrtyl and Danesova, 1999).

The pure chemical driving force represents the difference of the chemical potential externally supplied *<sup>e</sup>* (biochemical activity generated by the osteoclasts) with the chemical potential of the mineral of the solid phase, denoted min ; it can be estimated from the change of activity of the *H* cation (Silva and Ulm, 2002):

$$
\Delta\mu \coloneqq \mu\_e - \mu\_{\rm min} = R\theta \ln \frac{\left[H^+\right]\_{eq}^2}{\left[H^+\right]\_{ex}^2} \tag{5.9}
$$

This chemical driving force is the affinity conjugated to the superficial concentration of minerals, denoted ( ) *n t* in the sequel. The conversion to mechanical units is done, considering a density of HA <sup>3</sup> 3000 / *kg m* (5.1), hence / 20 *M MPa* , according to (Silva and Ulm, 2002); the negative value means that the dissolution of HA is chemically more favorable (bone resorption occurs).

Relying on the biochemical description given thereabove, bone remodeling is considered as a pure surface growth process. In order to analyze the influence of mechanical stress on bone remodeling, a simple geometrical model of a long bone as a hollow homogeneous cylinder is introduced, endowed with a linear elastic isotropic behavior (the interstitial fluid phase in the bone is presently neglected). This situation is representative of the diaphysal region of long bones (Cowin and Firoozbakhsh, 1981), such as the human femur (figure 6). According to experiments performed by (Currey, 1988), the elastic modulus is assumed to

scale uniformly versus the bone density according to

$$E = E\_{\text{max}} \rho\_S(t)^p \tag{5.10}$$

with *<sup>S</sup> t* the surface density of mineral, max *E GPa* 15 (Reilly and Burstein, 1975) the maximum value of the tensile modulus, and *p* a constant exponent, here taken equal to 3 (Currey, 1988; Ruinerman et al., 2005).

Following the representation theorems for isotropic scalar valued functions of tensorial arguments, the surface strain energy density , ; *<sup>S</sup> mech* **FNX***<sup>S</sup>* of mechanical origin is selected as a function of the curvature tensor invariants, viz the mean and Gaussian curvatures, the invariants of the surface Cauchy-Green tensor : .*<sup>t</sup>* **C FF** and of its square. The following simple form depending on the second invariant of the linearized part of **C I** 2**ε** is selected, adopting the small strain framework, viz, hence

$$\left.\boldsymbol{\nu}\right|\_{\text{mech}}^{S}\left(\tilde{\mathbf{c}}\right) = \frac{A}{2}Tr(\tilde{\mathbf{c}})^2 + B\left(\tilde{\mathbf{c}}:\tilde{\mathbf{c}}\right) \tag{5.11}$$

which is much larger compared to the characteristic size of the ultrastructure; the resorption phase takes typically 21 days (the complete remodeling cycle lasts 3 months). The osteoclasts, generate an acid environment causing simultaneously the dissolution of the

the sequel - and the degradation of the collageneous fraction of the material. The metabolic processes behind bone remodeling are very complicated, with kinetics of various chemical

The pure chemical driving force represents the difference of the chemical potential

min <sup>2</sup> : ln *eq*

This chemical driving force is the affinity conjugated to the superficial concentration of

(Silva and Ulm, 2002); the negative value means that the dissolution of HA is chemically

Relying on the biochemical description given thereabove, bone remodeling is considered as a pure surface growth process. In order to analyze the influence of mechanical stress on bone remodeling, a simple geometrical model of a long bone as a hollow homogeneous cylinder is introduced, endowed with a linear elastic isotropic behavior (the interstitial fluid phase in the bone is presently neglected). This situation is representative of the diaphysal region of long bones (Cowin and Firoozbakhsh, 1981), such as the human femur (figure 6). According to experiments performed by (Currey, 1988), the elastic modulus is assumed to

> max *<sup>p</sup> EE t*

maximum value of the tensile modulus, and *p* a constant exponent, here taken equal to 3

Following the representation theorems for isotropic scalar valued functions of tensorial

as a function of the curvature tensor invariants, viz the mean and Gaussian curvatures, the invariants of the surface Cauchy-Green tensor : .*<sup>t</sup>* **C FF** and of its square. The following simple form depending on the second invariant of the linearized part of **C I** 2**ε** is

*A*

*<sup>S</sup> t* the surface density of mineral, max *E GPa* 15 (Reilly and Burstein, 1975) the

3000 / *kg m* (5.1), hence

*R*

 

*e*

*<sup>e</sup>* (biochemical activity generated by the osteoclasts) with the chemical

*H*

*H*

in the sequel. The conversion to mechanical units

2

*ex*

 

*<sup>S</sup>* (5.10)

*mech* **FNX***<sup>S</sup>* of mechanical origin is selected

**ε ε εε** *Tr B* (5.11)

*Ca PO Ca OH* , abbreviated HA in

min ; it can be estimated from the

(5.9)

/ 20 *M MPa* , according to

is done,

mineral - hydroxyapatite, a strong basic mineral 3 4 2 2 <sup>3</sup>

substances, see (Petrtyl and Danesova, 1999).

potential of the mineral of the solid phase, denoted

change of activity of the *H* cation (Silva and Ulm, 2002):

scale uniformly versus the bone density according to

arguments, the surface strain energy density , ; *<sup>S</sup>*

selected, adopting the small strain framework, viz, hence

 <sup>2</sup> () : <sup>2</sup> *S mech*

considering a density of HA <sup>3</sup>

more favorable (bone resorption occurs).

(Currey, 1988; Ruinerman et al., 2005).

externally supplied

minerals, denoted ( ) *n t*

with 

Fig. 6. Modeling occurring during growth of the proximal end of the femur. Frontal section of the original proxima tibia is indicated as the stippled area. The situation after a growth of 21 days is superimposed. Bone formation (+) and bone resorption zones indicated [Weiss, 1988].

with **ε P**. .. **ε I***S rr* **ε ε e e** the surface strain (induced by the existing volumetric strain), and *A*,*B* mechanical properties of the surface, expressing versus the surface density of minerals and the maximum value of the traction modulus as (the Poisson ratio is selected as 0.3 )

$$A = \frac{E\_{\text{max}} \rho\_{\text{S}} \left(t\right)^{\text{3}} \nu}{\left(1 - 2\nu\right) \left(1 + \nu\right)}; \; B = \frac{E\_{\text{max}} \rho\_{\text{S}} \left(t\right)^{\text{3}}}{2\left(1 + \nu\right)}\tag{5.12}$$

As the surface of bone undergoes resorption, its mechanical properties are continuously changing from the bulk behavior, due to the decrease of mineral density as reflected in (5.10). The surface stress results from (5.11), (5.12) as

$$\tilde{\mathbf{T}} \equiv \tilde{\mathbf{o}} := \frac{\partial \boldsymbol{\mu}^{\rm S}\_{mech}}{\partial \tilde{\mathbf{e}}} = A\tilde{\mathbf{e}} + 2Btr\left(\tilde{\mathbf{e}}\right) \mathbf{I}\_{\rm S} \tag{5.13}$$

The unknowns of the remodeling problem are the normal velocity of the bone surface*V t <sup>N</sup>* , the surface density of minerals *<sup>S</sup> t* and its superficial concentration. We shall herewith simulate the resorption of a hollow bone submitted to a composite applied stress, consisting of the superposition of an axial and a radial component, as

$$\mathbf{o} = \sigma\_{rr}\mathbf{e}\_r \otimes \mathbf{e}\_r + \sigma\_{zz}\mathbf{e}\_z \otimes \mathbf{e}\_z \tag{5.14}$$

in the cylindrical basis **eee** *r z* , , ; this applied stress generates a preexisting homogeneous stress state within the bulk material, inducing a surface stress given by

Thermodynamics of Surface Growth with Application to Bone Remodeling 397

reality a rather complex chemical reaction (Blair, 1998) that is here simply modeled as a

 2 2 3 4 2 2 4 2 <sup>3</sup> *Ca PO Ca OH H Ca HPO H O* 8 10 6 2

*s i*

*t r*

 incorporating the density of minerals. The rate coefficient of dissolution of HA, namely the

The present model involves a dependency of the triplet of variables *rt t n t i S* , ,

solution of the set of equations (5.15) through (5.19) on a set of parameters, arising from


The evolution versus time of some variables of interest is next shown, considering a time scale conveniently expressed in days. Numerical simulations of bone resorption are to be performed for three stress levels in the normal physiological range,

 1 ,2 ,5 *MPa MPa MPa* . The surface velocity (Figure 7) shows an acceleration of the resorption process with time, which is enhanced by the stress level, as expected from the

The density and concentration vanish over long durations, meaning that the bone has been

stress level of 1MPa (Cowin, 2001). The superficial density of minerals and its concentration are both weakly dependent upon stress; the density of minerals decreases by a factor two (for low stresses; the resorption is enhanced by the applied stress) over a period of one

Considering an imposed stress function of time, the surface driving force is seen to vanish

 <sup>10</sup> 3/2 1/2 ( ) 9.4.10 *crit zz <sup>S</sup> t t nt*

 

This expression gives an order of magnitude of the stress level above which bone apposition (growth) shall take place; when the critical stress is reached, the chemical and mechanical

*t* , depending upon the density and concentration, given from (5.18),

(5.16)

An order of magnitude of the simulated radial surface velocity is about 10 /

obtained by time integration of the normal velocity expressed in (5.16).

0

0

*tn t r t tn t*

, is taken at room temperature from literature values available for CHA

0

is taken as unity, viz <sup>3</sup>

 

<sup>0</sup> *n* 1 mol.m

.

*m day* for a

exp *<sup>s</sup> <sup>S</sup>*

(5.19)

2.2.10 *s* (Hankermeyer et al., 2002).

single first order kinetic reaction

The kinetic equation is chosen as:

parameter

**5.3 Simulation results** 

*n t*

(carbonated HA, similar to bone), viz 4 1

initial conditions satisfied by those variables: - The initial concentration of minerals *n*<sup>0</sup>

higher magnitude of the driving force.

*zz* 

driving forces do balance, and the bone microstructure is stable.

completely dissolved (Figure 8).

month resorption period.

for a critical stress ( ) *crit*

(5.19) as

. **σ P σ** *zz z z* **e e**

The radial component of Eshelby stress **Σ***rr* is then easily evaluated from the preexisting homogeneous stress state. Straightforward calculations deliver then the driving force for surface remodeling, as the sum of a chemical and a mechanical contribution due to the applied axial stress:

$$\tilde{\Upsilon}\_{\mathbf{g}^{\mathcal{N}}} = \frac{1}{r\_i(t)} \left| \frac{1}{8} \Delta \mu m^{\sigma} \left( t \right) + \frac{A + 2B}{8 \left( A + B \right) B} \sigma\_{zz} \right|^2 \right| \tag{5.15}$$

with the material coefficients , *A B* given in (5.12), and the axial stress *zz* possibly function of time. A simple linear relation of the velocity of the growing surface to the driving force is selected, viz

$$V\_N\left(t\right) = C\tilde{\Upsilon}\_{\mathfrak{g}^N}\left(t\right) \tag{5.16}$$

with *C* a positive parameter; the positive sign is due to the velocity direction being opposite to the outer normal (the inner radius is increasing). The chemical contribution leads by itself to resorption, hence the normal velocity has to be negative; the mechanical contribution in (5.15) brings a positive contribution to the driving force for bone growth, corresponding to apposition of new bone when the neat balance of energy is favorable to bone growth. An estimate of the amplitude of the normal velocity is given from the expression of the rate of dissolution of HA in (5.8) as

$$I = \rho\_{\mathcal{R}} V\_{\mathcal{N}} \;/\; M = 10^{-8} \,\text{mol} \,\text{s}^{-1} \,\text{m}^{-2} \Rightarrow V\_{\mathcal{N}} = \text{[M }/\; \rho\_{\mathcal{R}} \approx 3.3 \,\text{10}^{-12} \,\text{m} \;/\text{s} \approx 0.286 \,\mu\text{m} \;/\text{ day}$$

selecting a molar mass *M* 1.004 / *kg mol* , following (Silva and Ulm, 2002). This value is an initial condition for the radius evolution (its rate is prescribed), leading to 23 2 1 *C mk* 3.5.10 . . *g s* ; it is however much lower compared to typical values of the bulk growth velocity, about 10 / *m day* .

The mass balance equation for the surface density of minerals *<sup>S</sup>* writes

$$
\rho\_S + \rho\_S \nabla\_S \tilde{\mathbf{V}} = \Gamma^S \rho\_S \tag{5.17}
$$

expressing as the following conservation law

$$\frac{\dot{\rho\_S}}{\rho\_S} - \frac{V\_N}{r\_i(t)} = \Gamma\_0^S \Leftrightarrow \rho\_S(t) = \rho\_s^0 \frac{r\_0}{r\_i(t)} \exp\left(\Gamma\_0^S t\right) \tag{5.18}$$

The initial surface density of minerals <sup>0</sup> *s S* 0 *t* , is evaluated from the bulk density of HA, viz <sup>3</sup> 3000 / *kg m* , and the estimated thickness of the attachment region of osteoclasts, about 7*m* (Blair, 1998), hence 0 22 *<sup>s</sup>* 2.1.10 / *kg <sup>m</sup>* .

The surface growth rate of mass 0 *<sup>S</sup>* is here assumed to be constant (it represents a datum) and can be identified to the rate of dissolution of HA, adopting the chemical reaction model of (Blair, 1998): 0 *<sup>S</sup>* is estimated by considering that 80% of the superficial minerals have been dissolved in a 2 months period, hence 7 1 <sup>0</sup> 2.2.10 *<sup>S</sup> <sup>s</sup>* . The dissolution of HA is in reality a rather complex chemical reaction (Blair, 1998) that is here simply modeled as a single first order kinetic reaction

$$\left[\mathrm{Ca}\_3\mathrm{(PO}\_4\mathrm{)}\_2\right]\_3\mathrm{Ca}(\mathrm{OH})\_2 + 8H^+ \to 10\mathrm{Ca}^{2+} + 6\mathrm{HPO}\_4^{2-} + 2H\_2\mathrm{O}$$

The kinetic equation is chosen as:

396 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

. **σ P σ** *zz z z* **e e**

The radial component of Eshelby stress **Σ***rr* is then easily evaluated from the preexisting homogeneous stress state. Straightforward calculations deliver then the driving force for surface remodeling, as the sum of a chemical and a mechanical contribution due to the

> 1 1 2 <sup>2</sup> 8 8 *gN zz*

*r t A BB* 

of time. A simple linear relation of the velocity of the growing surface to the driving force is

with *C* a positive parameter; the positive sign is due to the velocity direction being opposite to the outer normal (the inner radius is increasing). The chemical contribution leads by itself to resorption, hence the normal velocity has to be negative; the mechanical contribution in (5.15) brings a positive contribution to the driving force for bone growth, corresponding to apposition of new bone when the neat balance of energy is favorable to bone growth. An estimate of the amplitude of the normal velocity is given from the expression of the rate of

*g N* / 10 . . *N g* / 3.3.10 / 0.286 /

selecting a molar mass *M* 1.004 / *kg mol* , following (Silva and Ulm, 2002). This value is an initial condition for the radius evolution (its rate is prescribed), leading to 23 2 1 *C mk* 3.5.10 . . *g s* ; it is however much lower compared to typical values of the bulk

*m s m day*

. *<sup>S</sup>*

0 0 exp ( ) *S N S S S s*

*s S* 0

HA, viz <sup>3</sup> 3000 / *kg m* , and the estimated thickness of the attachment region of osteoclasts,

and can be identified to the rate of dissolution of HA, adopting the chemical reaction model

 

*V r t t*

 

<sup>0</sup> <sup>0</sup>

*<sup>S</sup>* is estimated by considering that 80% of the superficial minerals have

(5.18)

 *S SS S* 

*A B n t*

 

*Vt C t N g gN* (5.16)

*<sup>S</sup>* writes

**V** (5.17)

*t* , is evaluated from the bulk density of

<sup>0</sup> 2.2.10 *<sup>S</sup> <sup>s</sup>* . The dissolution of HA is in

*<sup>S</sup>* is here assumed to be constant (it represents a datum)

*zz* possibly function

 

(5.15)

*i*

with the material coefficients , *A B* given in (5.12), and the axial stress

8 12 <sup>12</sup> *J V M mol s m V JM*

*S i i*

 

*<sup>s</sup>* 2.1.10 / *kg <sup>m</sup>* .

*r t r t*

applied axial stress:

selected, viz

dissolution of HA in (5.8) as

growth velocity, about 10 /

about 7

of (Blair, 1998): 0

The initial surface density of minerals <sup>0</sup>

*m* (Blair, 1998), hence 0 22 

been dissolved in a 2 months period, hence 7 1

expressing as the following conservation law

The surface growth rate of mass 0

*m day* . The mass balance equation for the surface density of minerals

$$\frac{\partial n^{\sigma}(t)}{\partial t} = -\tilde{\gamma}\rho\_{s}(t)n^{\sigma}(t) = -\tilde{\gamma}\eta\_{i}(t)\frac{\rho\_{s}^{0}}{r\_{0}}\exp\left(\Gamma\_{0}^{S}t\right)n^{\sigma}(t)\tag{5.19}$$

incorporating the density of minerals. The rate coefficient of dissolution of HA, namely the parameter , is taken at room temperature from literature values available for CHA (carbonated HA, similar to bone), viz 4 1 2.2.10 *s* (Hankermeyer et al., 2002).

### **5.3 Simulation results**

The present model involves a dependency of the triplet of variables *rt t n t i S* , , solution of the set of equations (5.15) through (5.19) on a set of parameters, arising from initial conditions satisfied by those variables:


The evolution versus time of some variables of interest is next shown, considering a time scale conveniently expressed in days. Numerical simulations of bone resorption are to be performed for three stress levels in the normal physiological range, 1 ,2 ,5 *MPa MPa MPa* . The surface velocity (Figure 7) shows an acceleration of the resorption process with time, which is enhanced by the stress level, as expected from the higher magnitude of the driving force.

The density and concentration vanish over long durations, meaning that the bone has been completely dissolved (Figure 8).

An order of magnitude of the simulated radial surface velocity is about 10 / *m day* for a stress level of 1MPa (Cowin, 2001). The superficial density of minerals and its concentration are both weakly dependent upon stress; the density of minerals decreases by a factor two (for low stresses; the resorption is enhanced by the applied stress) over a period of one month resorption period.

Considering an imposed stress function of time, the surface driving force is seen to vanish for a critical stress ( ) *crit zz t* , depending upon the density and concentration, given from (5.18), (5.19) as

$$
\sigma\_{zz}^{c\bar{\tau}t}(t) \approx 9.4.10^{10} \,\rho\_{\bar{\rm S}}(t)^{3/2} \, n^{\sigma} \left(t\right)^{1/2} \tag{5.16}
$$

This expression gives an order of magnitude of the stress level above which bone apposition (growth) shall take place; when the critical stress is reached, the chemical and mechanical driving forces do balance, and the bone microstructure is stable.

Thermodynamics of Surface Growth with Application to Bone Remodeling 399

growth will occur due to mineralization (the chemical driving force in (5.9) favors apposition of new bone on the surface), as reflected by the simulated decrease of the internal

Fig. 9. Evolution of the internal radius of the diaphysis of the human femur (in microns)

Apposition of new bone would occur in the absence of mechanical stimulus, under the influence of a pure chemical driving force; in that case, the internal radius will decrease very rapidly (Figure 9) and tends to an asymptotic value (for long times) after about two weeks growth. For a non vanishing axial stress above the critical stress in (5.16), the driving force is negative in the first growth period, and becomes thereafter positive due to the decrease of the surface density of minerals, indicating that growth takes over from bone resorption. Hence, the developed model is able to encompass both situations of growth and resorption, according to the level of applied stress (the nature of the stress, compressive or under traction, does not play a role according to (5.15)), which determine the mechanical

Surface growth is by essence a pluridisciplinary field, involving interactions between the physics and mechanics of surfaces and transport phenomena. The literature survey shows different strategies for treating superficial interactions, hence recognizing that no unitary viewpoint yet exists. The present contribution aims at providing a pluridisciplinary

*zz MPa* .

versus time. Applied stress above the critical stress level: 0.02

contribution of the overall driving force for growth.

**6. Concluding remarks** 

approach of surface growth focusing on

*zz MPa* lying slightly above the critical stress expressed in (5.16),

For an applied stress 0.2

radius over the first week (Figure 9).

Fig. 7. Evolution vs. time of the surface growth velocity for three stress levels: 1 *zz MPa* (thick line), 2 *zz MPa* (dashed line), 5 *zz MPa* (dash-dotted line).

Fig. 8. Evolution of the superficial density of HA versus time for three stress levels. 1 *zz MPa* (thick line), 2 *zz MPa* (dashed line), 5 *zz MPa* (dash-dotted line).

Fig. 7. Evolution vs. time of the surface growth velocity for three stress levels: 1

Fig. 8. Evolution of the superficial density of HA versus time for three stress levels.

*zz MPa* (dashed line), 5

*zz MPa* (dash-dotted line).

*zz MPa* (dash-dotted line).

*zz MPa* (dashed line), 5

(thick line), 2 

1

*zz MPa* (thick line), 2

*zz MPa*

For an applied stress 0.2 *zz MPa* lying slightly above the critical stress expressed in (5.16), growth will occur due to mineralization (the chemical driving force in (5.9) favors apposition of new bone on the surface), as reflected by the simulated decrease of the internal radius over the first week (Figure 9).

Fig. 9. Evolution of the internal radius of the diaphysis of the human femur (in microns) versus time. Applied stress above the critical stress level: 0.02 *zz MPa* .

Apposition of new bone would occur in the absence of mechanical stimulus, under the influence of a pure chemical driving force; in that case, the internal radius will decrease very rapidly (Figure 9) and tends to an asymptotic value (for long times) after about two weeks growth. For a non vanishing axial stress above the critical stress in (5.16), the driving force is negative in the first growth period, and becomes thereafter positive due to the decrease of the surface density of minerals, indicating that growth takes over from bone resorption.

Hence, the developed model is able to encompass both situations of growth and resorption, according to the level of applied stress (the nature of the stress, compressive or under traction, does not play a role according to (5.15)), which determine the mechanical contribution of the overall driving force for growth.
