**1. Introduction**

368 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

Pan S et al. (2007). Theoretical approach to evaluate thermodiffusion in aqueous alkanol

Parola, A., Piazza, R. (2004). Particle thermophoresis in liquids. *The European Physical Journal*,

Ross, S. and Morrison, I. D. (1988) *Colloidal Systems and Interfaces,* John Wiley and Sons, New

Ruckenstein, E. (1981). Can phoretic motion be treated as interfacial tension gradient driven

Schimpf, M. E., Semenov, S. N. (2004). Thermophoresis of Dissolved Molecules and

Gradient. *Physical Review E*, Vol. 69, No. 1 (October 2004), 011201 (8 pages). Semenov, S. N., Schimpf, M. E. (2005). Molecular thermodiffusion (thermophoresis) in liquid mixtures. *Physical Review E*, Vol. 72, No. 4 (October *2005)* 041202 (9 pages). Semenov, S. N., Schimpf, M. E. (2000) . Mechanism of Polymer Thermophoresis in

Semenov, S. N., Schimpf, M. E. (2009). Mass Transport Thermodynamics in Nonisothermal

Semenov, S. N. (2010). Statistical thermodynamic expression for the Soret coefficient

Weinert, F. M., Braun, D. (2008). Observation of Slip Flow in Thermophoresis. *Physical* 

Semenov, S. N., Schimpf, M. E. (2011). Internal degrees of freedom, molecular symmetry and thermodiffusion. *Comptes Rendus* Mecanique, doi:10.1016/j.crme.2011.03.011. Semenov, S. N., Schimpf, M. E. (2011).Thermodynamics of mass transport in diluted

Semenov, S. N. (2011). Statistical thermodynamics of material transport in non-isothermal

Wiegand, S., Kohler, W. (2002). Measurements of transport coefficients by an optical grating

Wittko, G., Köhler, W. (2005). Universal isotope effect in thermal diffusion of mixtures

*Europhysics Letters*, Vol. 90, No. 5 (June 2010), 56002 (6 pages).

*Review Letters.* Vol. 101 (October 2008), 168301(4 pages).

*Rendus* Mecanique, doi:10.1016/j.crme.2011.03.002.

Springer, Berlin, Germany.

123, No. 6 (June 2005), 014506 (6 pages).

binary molecular systems. Submitted to *Europhysics Letters.* 

pages).

York, USA.

1981), 77-82.

200),9935-9942.

1045-1054.

Vol.15, No. 11(November2004), 255-263.

solutions. *The Journal of Chemical Physics,* Vol. *126,* No. 1 (January 2007), 014502 (12

phenomena? *The Journal of Colloid and Interface Science,* Vol. 83No. 1 (September

Polymers: Consideration of the Temperature-Induced Macroscopic Pressure

Nonaqueous Solvents. *The Journal of Physical Chemistry B,* Vol. 104, No. 42 (July

Molecular Liquid Mixtures. *Physics – Uspekhi*, Vol. 52, No. 11 (November 2090),

suspensions of charged particles in non-isothermal liquid electrolytes. *Comptes* 

technique. In: *Thermal Nonequilibrium Phenomena in Fluid Mixtures* (Lecture Notes in Physics, Vol. 584, W. Kohler, S. Wiegand (Eds.), 189-210, ISBN 3-540-43231-0,

containing cyclohexane and cyclohexane-d12. *The Journal of Chemical Physics,* Vol.

In physics, surface growth classically refers to processes where material reorganize on the substrate onto which it is deposited (like epitaxial growth), but principally to phenomena associated to phase transition, whereby the evolution of the interface separating the phases produces a crystal (Kessler, 1990; Langer, 1980). From a biological perspective, *surface growth* refers to mechanisms tied to accretion and deposition occurring mostly in hard tissues, and is active in the formation of teeth, seashells, horns, nails, or bones (Thompson, 1992). A landmark in this field is Skalak (Skalak et al., 1982, 1997) who describe the growth or atrophy of part of a biological body by the accretion or resorption of biological tissue lying on the surface of the body. Surface growth of biological tissues is a widespread situation, with may be classified as either fixed growth surface (e.g. nails and horns) or moving growing surface (e.g. seashells, antlers). Models for the kinematics of surface growth have been developed in (Skalak et al., 1997), with a clear distinction between cases of fixed and moving growth surfaces, see (Ganghoffer et al., 2010a,b; Garikipati, 2009) for a recent exhaustive literature review.

Following the pioneering mechanical treatments of elastic material surfaces and surface tension by (Gurtin and Murdoch, 1975; Mindlin, 1965), and considering that the boundary of a continuum displays a specific behavior (distinct from the bulk behavior), subsequent contributions in this direction have been developed in the literature (Gurtin and Struthers, 1990; Gurtin, 1995, Leo and Sekerka, 1989) for a thermodynamical approach of the surface stresses in crystals; configurational forces acting on interfaces have been considered e.g. in (Maugin, 1993; Maugin and Trimarco, 1995) – however not considering surface stress -, and (Gurtin, 1995; 2000) considering specific balance laws of configurational forces localized at interfaces.

Biological evolution has entered into the realm of continuum mechanics in the 1990's, with attempts to incorporate into a continuum description time-dependent phenomena, basically consisting of a variation of material properties, mass and shape of the solid body. One outstanding problem in developmental biology is indeed the understanding of the factors that may promote the generation of biological form, involving the processes of growth (change of mass), remodeling (change of properties), and morphogenesis (shape changes), a classification suggested by Taber (Taber, 1995).

The main focus in this chapter is the setting up of a modeling platform relying on the thermodynamics of surfaces (Linford, 1973) and configurational mechanics (Maugin, 1993)

Thermodynamics of Surface Growth with Application to Bone Remodeling 371

*a t* ( ,) **x** . The particular form of the flux and source depend on the nature of the considered extensive quantity, as shall appear in the forthcoming balance laws. We consider a system including n constituents undergoing r chemical reactions; the local variations of the partial

**x** *t* the local production (or destruction) of

*<sup>k</sup>* , obey the local balance law (Vidal et al., 1994)

, such that the variation of the mass *<sup>k</sup> dm* of the

(3)

(4)

1.. *k k r M J* 

1 1

*j j*

*k kk*

**<sup>J</sup> u 0** , as

 

**u u** (5)

**<sup>J</sup>** (6)

its production term, given by De Donder definition

*n n*

 *<sup>k</sup>* the

. The molar masses *Mk*

*<sup>k</sup>***u** and a diffusive flux

1

 

*n k k* 

. In this viewpoint,

1..

*r*

 

 

**u J** (2)

*kk k k*

*<sup>M</sup> <sup>J</sup> <sup>t</sup>*

**u u** the local barycentric velocity, *Mk* the molar mass, and

1.. *kk k* , k=1..n *r*

*k k*

the system is in fact closed, since the balance law satisfied by the global density

 

0, 1..

*M r*

 

*k k*

 

1 1.. . .

This balance law does not involve any source term for the total density. Instead of using the partial densities of the system constituents, one can write balance equations for the number of moles of constituent k, *nmM kk k* / , with *mk* the mass of the same constituent. The molar concentration is defined as *c nV k k* / , its inverse being called the partial molar

> *k i k k n n div t t*

*k r M J <sup>t</sup>*

 

*n*

 

. *<sup>k</sup>*

*dm M*

denotes the degree of advancement of reaction

1

*k*

volume. The partial mole number *nk* satisfies the balance equation

*t* 

Observe that the total flux of mass is the sum of a convective flux

*<sup>k</sup>* **J** ; the mass production is identified as the contribution

writes (Vidal et al., 1994) accounting for the relation

with *<sup>k</sup>* **<sup>J</sup>** the flux of species k and *i k <sup>n</sup>*

of the rate of progress of the jth chemical reaction

*n*

with ( ,) *<sup>a</sup>* **J x** *t* the flux density of *a t* ( ,) **x** and ( ,) *<sup>a</sup>*

density of a given constituent k, quantity

1

*k* 

*k k*

stoechiometric coefficients in the reaction

species k due to chemical reactions expresses as

satisfy the global conservation law (due to Lavoisier)

*n*

1 :

 

with

wherein

for the treatment of surface growth phenomena in a biomechanical context. A typical situation is the external remodeling in long bones, which is induced by genetic and epigenetic factors, such as mechanical and chemical stimulations. The content of the chapter is the following: the thermodynamics of coupled irreversible phenomena is briefly reviewed, and balance laws accounting for the mass flux and the mass source associated to growth are expressed (section 2). Evolution laws for a growth tensor (the kinematic multiplicative decomposition of the transformation gradient into a growth tensor and an accommodation tensor is adopted) in the context of volumetric growth are formulated, considering the interactions between the transport of nutrients and the mechanical forces responsible for growth. As growth deals with a modification of the internal structure of the body in a changing referential configuration, the language and technique of Eshelbian mechanics (Eshelby, 1951) are adopted and the driving forces for growth are identified in terms of suitable Eshelby stresses (Ganghoffer and Haussy, 2005; Ganghoffer, 2010a). Considering next surface growth, the thermodynamics of surfaces is first exposed as a basis for a consistent treatment of phenomena occurring at a growing surface (section 3), corresponding to the set of generating cells in a physiological context. Material forces for surface growth are identified (section 4), in relation to a surface Eshelby stress and to the curvature of the growing surface. Considering with special emphasis bone remodeling (Cowin, 2001), a system of coupled field equations is written for the superficial density of minerals, their concentration and the surface velocity, which is expressed versus a surface material driving force in the referential configuration. The model is able to describe both bone growth and resorption, according to the respective magnitude of the chemical and mechanical contributions to the surface driving force for growth (Ganghoffer, 2010a). Simulations show the shape evolution of the diaphysis of the human femur. Finally, some perspectives in the field of growth of biological tissues are mentioned.

As to notations, vectors and tensors are denoted by boldface symbols. The inner product of two second order tensors is denoted . *ik kj ij* **A B** *A B* . The material derivative of any function is denoted by a superposed dot.

## **2. Thermodynamics of irreversible coupled phenomena: a survey**

We consider multicomponent systems, mutually interacting by chemical reactions. Two alternative viewpoints shall be considered: in the first viewpoint, the system is closed, which in consideration of growth phenomena means that the nutrients are included into the overall system. The second point of view is based on the analysis of a solid body as an open system exchanging nutrients with its surrounding; hence growth shall be accounted for by additional source terms and convective fluxes.

#### **2.1 Multiconstituents irreversible thermodynamics**

We adopt the thermodynamic framework of open systems irreversible thermodynamics, which shall first be exposed in a general setting, and particularized thereafter for growing continuum solid bodies. Recall first that any extensive quantity *A* with volumetric density *aa t* ( ,) **x** satisfies a prototype balance law of the form

$$\frac{\partial \mathfrak{a}(\mathbf{x},t)}{\partial t} = -\nabla \cdot \mathbf{J}\_a(\mathbf{x},t) + \sigma\_a(\mathbf{x},t) \tag{1}$$

for the treatment of surface growth phenomena in a biomechanical context. A typical situation is the external remodeling in long bones, which is induced by genetic and epigenetic factors, such as mechanical and chemical stimulations. The content of the chapter is the following: the thermodynamics of coupled irreversible phenomena is briefly reviewed, and balance laws accounting for the mass flux and the mass source associated to growth are expressed (section 2). Evolution laws for a growth tensor (the kinematic multiplicative decomposition of the transformation gradient into a growth tensor and an accommodation tensor is adopted) in the context of volumetric growth are formulated, considering the interactions between the transport of nutrients and the mechanical forces responsible for growth. As growth deals with a modification of the internal structure of the body in a changing referential configuration, the language and technique of Eshelbian mechanics (Eshelby, 1951) are adopted and the driving forces for growth are identified in terms of suitable Eshelby stresses (Ganghoffer and Haussy, 2005; Ganghoffer, 2010a). Considering next surface growth, the thermodynamics of surfaces is first exposed as a basis for a consistent treatment of phenomena occurring at a growing surface (section 3), corresponding to the set of generating cells in a physiological context. Material forces for surface growth are identified (section 4), in relation to a surface Eshelby stress and to the curvature of the growing surface. Considering with special emphasis bone remodeling (Cowin, 2001), a system of coupled field equations is written for the superficial density of minerals, their concentration and the surface velocity, which is expressed versus a surface material driving force in the referential configuration. The model is able to describe both bone growth and resorption, according to the respective magnitude of the chemical and mechanical contributions to the surface driving force for growth (Ganghoffer, 2010a). Simulations show the shape evolution of the diaphysis of the human femur. Finally, some

perspectives in the field of growth of biological tissues are mentioned.

**2. Thermodynamics of irreversible coupled phenomena: a survey** 

is denoted by a superposed dot.

additional source terms and convective fluxes.

**2.1 Multiconstituents irreversible thermodynamics** 

*aa t* ( ,) **x** satisfies a prototype balance law of the form

*t*

As to notations, vectors and tensors are denoted by boldface symbols. The inner product of two second order tensors is denoted . *ik kj ij* **A B** *A B* . The material derivative of any function

We consider multicomponent systems, mutually interacting by chemical reactions. Two alternative viewpoints shall be considered: in the first viewpoint, the system is closed, which in consideration of growth phenomena means that the nutrients are included into the overall system. The second point of view is based on the analysis of a solid body as an open system exchanging nutrients with its surrounding; hence growth shall be accounted for by

We adopt the thermodynamic framework of open systems irreversible thermodynamics, which shall first be exposed in a general setting, and particularized thereafter for growing continuum solid bodies. Recall first that any extensive quantity *A* with volumetric density

> ( ,) . ( ,) ( ,) *a a a t t t*

**<sup>x</sup> <sup>J</sup> x x** (1)

with ( ,) *<sup>a</sup>* **J x** *t* the flux density of *a t* ( ,) **x** and ( ,) *<sup>a</sup>* **x** *t* the local production (or destruction) of *a t* ( ,) **x** . The particular form of the flux and source depend on the nature of the considered extensive quantity, as shall appear in the forthcoming balance laws. We consider a system including n constituents undergoing r chemical reactions; the local variations of the partial density of a given constituent k, quantity *<sup>k</sup>* , obey the local balance law (Vidal et al., 1994)

$$\frac{\partial \rho\_k}{\partial t} = -\nabla \cdot \left(\rho\_k \mathbf{u} + \mathbf{J}\_k\right) + M\_k \sum\_{a=1..r} \nu\_{ak} l\_a \tag{2}$$

with 1 1 : *n k k k* **u u** the local barycentric velocity, *Mk* the molar mass, and *<sup>k</sup>* the stoechiometric coefficients in the reaction , such that the variation of the mass *<sup>k</sup> dm* of the species k due to chemical reactions expresses as

$$d\mathfrak{m}\_k = \mathcal{M}\_k \sum\_{\alpha=1\dots r} \nu\_{\alpha k} \xi\_{\alpha \prime} \quad \text{k} = \mathbf{1} \dots \mathbf{n} \tag{3}$$

wherein denotes the degree of advancement of reaction . The molar masses *Mk* satisfy the global conservation law (due to Lavoisier)

$$\sum\_{k=1}^{n} \nu\_{\alpha k} M\_k = 0, \quad \alpha = 1..r \tag{4}$$

Observe that the total flux of mass is the sum of a convective flux *<sup>k</sup>***u** and a diffusive flux *<sup>k</sup>* **J** ; the mass production is identified as the contribution 1.. *k k r M J* . In this viewpoint,

the system is in fact closed, since the balance law satisfied by the global density 1 *n k k* 

writes (Vidal et al., 1994) accounting for the relation 1 1 *n n k kk j j* **<sup>J</sup> u 0** , as

$$\frac{\partial \rho}{\partial t} = -\nabla. \left(\rho \mathbf{u}\right) + \sum\_{k=1}^{n} \sum\_{\alpha=1\dots r} M\_k \nu\_{\alpha k} l\_\alpha \equiv -\nabla. \left(\rho \mathbf{u}\right) \tag{5}$$

This balance law does not involve any source term for the total density. Instead of using the partial densities of the system constituents, one can write balance equations for the number of moles of constituent k, *nmM kk k* / , with *mk* the mass of the same constituent. The molar concentration is defined as *c nV k k* / , its inverse being called the partial molar volume. The partial mole number *nk* satisfies the balance equation

$$\frac{\partial \mathbf{n}\_k}{\partial t} = -\text{div}\mathbf{J}\_k + \frac{\partial\_i \mathbf{n}\_k}{\partial t} \tag{6}$$

with *<sup>k</sup>* **<sup>J</sup>** the flux of species k and *i k <sup>n</sup> t* its production term, given by De Donder definition of the rate of progress of the jth chemical reaction

$$\frac{\partial \hat{\jmath}\_i n\_k}{\partial t} = \sum\_{j=1}^r \nu\_{kj} \dot{\tilde{\varphi}}\_j \tag{7}$$

Thermodynamics of Surface Growth with Application to Bone Remodeling 373

.. : *<sup>k</sup>*

One has assumed in this alternative that the mechanical power : *eq*

which is conveniently rewritten in terms of Helmholtz free energy density :

**2.2 General balance laws accounting for mass production due to growth** 

*k* **σ ε**

equality combined with the second principle, equality . *<sup>q</sup>*

 

irreversible mechanical power spent and to chemical reactions.

energy. The contribution : / *k k*

accounted for in a global manner as a source term.

across the closed surface *<sup>t</sup>* expresses then as

 

> 

 

Hence, the internal entropy production is identified as

1 11 <sup>1</sup>

.: / *<sup>q</sup> k k k*

 

*<sup>q</sup> <sup>i</sup>* : / *k k*

*V n* is identified to the term

*k*

*k*

*k*

**<sup>J</sup> σ ε** (17)

**<sup>J</sup> σ ε** (16)

**<sup>J</sup> σ ε** (15)

*<sup>u</sup>* **<sup>J</sup> σ ε** *V n* (14)

**J J σ ε** (13)

 

 

*s s*

 

*i*

**<sup>J</sup>** (the entropy flux

*u Ts* as

*w* **σ u** does not

*w* . Previous

*i q k j j k j*

which is due to the gradient of intensive variables (temperature, chemical potential), to the

An alternative to the previous writing of the internal entropy production bearing the name of Clausius-Duhem inequality is frequently used; as a starting point, the first principle is

include a flux contribution, hence only the heat diffusion contributes to the flux of internal

resumes to the sole heat flux), delivers after a few manipulations the variation of the internal

*u s Ts V n*

*<sup>i</sup> <sup>q</sup>* : / *k k*

*<sup>i</sup> <sup>q</sup>* .: / *k k*

*s s V n* 

This is at variant with the point of view adopted next, which consists in insulating a growing solid body from the external nutrients, identified as one the chemical species, but

In the case of mass being created / resorbed within a solid body considered as an open system from a general thermodynamic point of view, one has to account for a source term

 being produced (by a set of generating cells) at each point within the time varying volume *<sup>t</sup>* ; a convective term is also added, corresponding to the transport of nutrients by the velocity field of the underlying continuum. For any quantity a, the convective flux is locally defined in terms of its surface density as **F v** ( ) *a a* ; the overall convective flux of a

*s u Ts V n* 

*s w A V* 

and of the internal entropy production

written as

energy as

The two previous equalities enter into Gibbs relation as

$$
\rho \dot{\mu} = \theta \rho \dot{\mathbf{s}}\_{\varepsilon} + \theta \rho \dot{\mathbf{s}}\_{i} + \mathbf{o} \, : \, \dot{\mathbf{c}} - \sum\_{k} \mu\_{k} \frac{\mathcal{P}}{M} \text{div} \mathbf{J}\_{k} + \sum\_{k} \mu\_{k} \frac{\mathcal{P}}{M} \sum\_{j} \nu\_{kj} \dot{\tilde{\xi}}\_{j} \tag{8}
$$

with the temperature and *<sup>k</sup>* the chemical potential of constituent k. The chemical affinity in the sense of De Donder is defined as the force conjugated to the rate *<sup>j</sup>* 

$$A\_j = -\sum\_k \left(\mu\_k \frac{\rho}{M}\right) \nu\_{kj} = -\sum\_k \left(\mu\_k \nmid V\right) \nu\_{kj} \tag{9}$$

Hence, Gibbs relation can be rewritten in order to highlight the variation of entropy

$$\rho \dot{\mathbf{s}} = \frac{1}{\theta} \rho \dot{\mathbf{u}} - \frac{1}{\theta} \mathbf{o} : \dot{\mathbf{c}} + \sum\_{k} \frac{1}{\theta} \left( \mu\_{k} \frac{1}{V} \right) \text{div} \mathbf{J}\_{k} + \sum\_{j} \frac{1}{\theta} A\_{j} \dot{\check{\mathbf{s}}}\_{j} \tag{10}$$

The local balance of internal energy traduces the first principle of thermodynamics as

$$
\rho \dot{u} = -\nabla\_\cdot \mathbf{J}\_q + \rho \dot{w}
$$

with *<sup>q</sup>* **J** the heat flux, and the term *w* is relative to all forms of work. One shall isolate the flux-like contributions in the entropy variation, which after a few transformations writes

$$\rho \dot{\mathbf{s}} = \rho \dot{\mathbf{s}}\_{\varepsilon} + \rho \dot{\mathbf{s}}\_{i} = \frac{1}{\theta} \nabla \cdot \mathbf{J}\_{q} + \sum\_{k} \nabla \cdot \left( \mathbf{J}\_{k} \frac{\mu\_{k}}{\theta} \frac{1}{V} \right) - \sum\_{k} \mathbf{J}\_{k} \cdot \nabla \left( \frac{\mu\_{k}}{\theta} \frac{1}{V} \right) + \frac{1}{\theta} (\rho \dot{w} - \mathbf{o} : \dot{\mathbf{e}}) + \sum\_{j} \frac{1}{\theta} A\_{j} \dot{\xi}\_{j} \dot{\xi}\_{j}$$

The contribution **σ ε**: / (involving the virtual power of internal forces) is further decomposed into

$$\frac{1}{\theta}\mathbf{o} : \dot{\mathbf{c}} = \frac{1}{\theta}\mathbf{o} : \nabla \mathbf{u} = \nabla. \left(\frac{1}{\theta}\mathbf{o} : \mathbf{u}\right) - \mathbf{u}. \nabla. \left(\frac{1}{\theta}\mathbf{o}\right).$$

Hence, the rate of the entropy density decomposes into

$$\begin{aligned} \rho \dot{\mathbf{s}} &= \rho \dot{\mathbf{s}}\_c + \rho \dot{\mathbf{s}}\_i = -\nabla \cdot \left( \mathbf{J}\_q + \sum\_k \mathbf{J}\_k \frac{\mu\_k}{\theta} \frac{\mathbf{1}}{V} \right) + \mathbf{J}\_q \cdot \nabla \left( \frac{\mathbf{1}}{\theta} \right) - \\ \sum\_k \mathbf{J}\_k \cdot \nabla \left( \frac{\mu\_k}{\theta} \frac{\mathbf{1}}{V} \right) + \frac{1}{\theta} (\rho \dot{w} - \mathbf{o} : \dot{\mathbf{e}}) + \sum\_j \frac{1}{\theta} A\_j \dot{\xi}\_j \end{aligned} \tag{11}$$

This writing allows the identification of the divergential contribution to the exchange entropy, hence to the entropy flux

$$\mathbf{J}\_s = \mathbf{J}\_q + \sum\_k \mathbf{J}\_k \frac{\mu\_k}{\theta} \frac{1}{V} \tag{12}$$

and of the internal entropy production

372 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

*<sup>r</sup> i k*

*n t*

*us s div*

affinity in the sense of De Donder is defined as the force conjugated to the rate *<sup>j</sup>*

Hence, Gibbs relation can be rewritten in order to highlight the variation of entropy

The local balance of internal energy traduces the first principle of thermodynamics as

 

The two previous equalities enter into Gibbs relation as

 

Hence, the rate of the entropy density decomposes into

*ss s*

*k*

  

*V*

with *<sup>q</sup>* **J** the heat flux, and the term

 

entropy, hence to the entropy flux

decomposed into

The contribution **σ ε**: /

the temperature and

with  1

 

: *ei k k k kj j*

 **σ ε J**

*<sup>j</sup> k kj k k* / *<sup>j</sup> k k A V M* 

1 1 11 1

.

*u w* **J***<sup>q</sup>*

*s u div A*

 

flux-like contributions in the entropy variation, which after a few transformations writes

 **JJ J σ ε**

*ss s w A V V*

> 11 1 1 : : . : ..

This writing allows the identification of the divergential contribution to the exchange

*k j j*

*sq k*

. :

**J σ ε**

*k j*

1 1 1

*ei q k q k*

 

 **σε σ <sup>u</sup> <sup>σ</sup> u u <sup>σ</sup>**

1 1 . .

*V*

 

*k*

**JJ J**

*w A*

1 *<sup>k</sup>*

*<sup>k</sup> V* 

*k kj*

 

: *k k j j k j*

 <sup>1</sup> 1 11 <sup>1</sup> .. . : *k k e i q k k j j*

*k k j*

 

(involving the virtual power of internal forces) is further

 

> 

**JJ J** (12)

 

 

**σ ε <sup>J</sup>** (10)

*V*

*M M*

(8)

*j*

*kj j*

(7)

 

*<sup>k</sup>* the chemical potential of constituent k. The chemical

(9)

*w* is relative to all forms of work. One shall isolate the

 

 

 

(11)

$$\rho \dot{\mathbf{s}}\_i = \mathbf{J}\_q \cdot \nabla \left(\frac{1}{\theta}\right) - \sum\_k \mathbf{J}\_k \cdot \nabla \left(\frac{\mu\_k}{\theta} \frac{1}{V}\right) - \frac{1}{\theta} (\mathbf{o} : \dot{\mathbf{e}} - \rho \dot{w}) + \sum\_j \frac{1}{\theta} A\_j \dot{\xi}\_j \tag{13}$$

which is due to the gradient of intensive variables (temperature, chemical potential), to the irreversible mechanical power spent and to chemical reactions.

An alternative to the previous writing of the internal entropy production bearing the name of Clausius-Duhem inequality is frequently used; as a starting point, the first principle is written as

$$
\rho \dot{\boldsymbol{\mu}} = -\nabla \cdot \mathbf{J}\_q + \mathbf{o} : \dot{\mathbf{e}} + \sum\_k \left( \mu\_k \, / \, V \right) \dot{\boldsymbol{n}}\_k \tag{14}
$$

One has assumed in this alternative that the mechanical power : *eq w* **σ u** does not include a flux contribution, hence only the heat diffusion contributes to the flux of internal energy. The contribution : / *k k k* **σ ε** *V n* is identified to the term *w* . Previous

equality combined with the second principle, equality . *<sup>q</sup> i s s* **<sup>J</sup>** (the entropy flux

resumes to the sole heat flux), delivers after a few manipulations the variation of the internal energy as

$$\rho \dot{\boldsymbol{\mu}} = \left( \theta \rho \dot{\mathbf{s}} - \mathbf{J}\_q \frac{\nabla \theta}{\theta} - T \,\rho \dot{\mathbf{s}}\_i \right) + \mathbf{o} \,\mathbf{r} \,\dot{\mathbf{z}} + \sum\_k \left( \mu\_k \,/\, V \right) \dot{\boldsymbol{n}}\_k \tag{15}$$

Hence, the internal entropy production is identified as

$$\rho \partial \rho \dot{\mathbf{s}}\_i = -\rho \left(\dot{\mu} - T\dot{\mathbf{s}}\right) - \mathbf{J}\_q \frac{\nabla \theta}{\theta} + \mathbf{o} : \dot{\mathbf{c}} + \sum\_k \left(\mu\_k \,/\, V\right) \dot{n}\_k \tag{16}$$

which is conveniently rewritten in terms of Helmholtz free energy density : *u Ts* as

$$\boldsymbol{\dot{\theta}} \cdot \boldsymbol{\theta} \dot{\boldsymbol{\varphi}}\_i = -\boldsymbol{\rho} \left( \dot{\boldsymbol{\nu}} + \boldsymbol{s} \dot{\boldsymbol{\theta}} \right) - \mathbf{J}\_q \cdot \frac{\nabla \boldsymbol{\theta}}{\boldsymbol{\theta}} + \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} + \sum\_k \left( \mu\_k \;/\; V \right) \dot{\boldsymbol{n}}\_k \tag{17}$$

This is at variant with the point of view adopted next, which consists in insulating a growing solid body from the external nutrients, identified as one the chemical species, but accounted for in a global manner as a source term.

#### **2.2 General balance laws accounting for mass production due to growth**

In the case of mass being created / resorbed within a solid body considered as an open system from a general thermodynamic point of view, one has to account for a source term being produced (by a set of generating cells) at each point within the time varying volume *<sup>t</sup>* ; a convective term is also added, corresponding to the transport of nutrients by the velocity field of the underlying continuum. For any quantity a, the convective flux is locally defined in terms of its surface density as **F v** ( ) *a a* ; the overall convective flux of a across the closed surface *<sup>t</sup>* expresses then as

Thermodynamics of Surface Growth with Application to Bone Remodeling 375

*<sup>D</sup> adx dx a a dA*

The first term on the r.h.s. corresponds to mass production, the second contribution to convection of the produced mass through the boundary *<sup>t</sup>* , and the third contribution to diffusion through the boundary of the moving volume *<sup>t</sup>* . One can see that only the relative velocity of particles w.r. to the surface velocity matters. Combining this identity

. .

**w n vw J n** (24)

.

**<sup>v</sup>** (28)

.**v** , and for the source and flux of mass reflected by the right hand side

 

**vJ v J** (25)

**J n** (26)

**<sup>J</sup> <sup>v</sup>** (27)

the physical source of mass, and *m* : . **m n** the scalar

 

**m v** (29)

, the

by the

*<sup>a</sup> dx a dA dx a a dA*

The corresponding local balance law is obtained after elimination of the velocity **w** , hence

**Mass balance**: the mass balance equation is deduced from the identification *a*

*tt t*

The **mass balance in Eulerian format** is given in terms of the actual density

.

*va <sup>a</sup> adiv div a div a div a*

*<sup>D</sup> dx dx dA*

. *<sup>v</sup> div*

following reasoning: we first write the general form of the balance of mass in physical space as

*t t t tt D D dx dx dx mds dx*

physical mass flux across the boundary, projection of the flux (vector) **m** . The previous balance law is quite general, as we account for both the variation of the integration volume

. . *<sup>D</sup>*

law has been obtained in (Epstein and Maugin, 2000) starting form its Lagrangian

 

**v x** the Eulerian velocity, which proves identical to (27); the same balance

 

*tt t*

*Dt*

*t t tt*

*t t*

*t* 

*Dt* 

*Dt*

The strong form of the balance law of mass writes finally

*Dt Dt*

of (28). Localization of previous integral equation gives

**x**,*t* the actual density,

*X*

*t* **x**

actual density. Hence, (23) gives

*t*

with (22) gives

with 

through the term

with ( , ):

counterpart.

*t*

.

**vw J n** (23)

$$\phi\_{conv}(a) = -\int\_{\partial\Omega\_t} a(\mathbf{v} - \mathbf{w}) \, \mathbf{n} dA \tag{18}$$

(the minus accounts for the unit exterior normal **n** ). This diffusive flux corresponds to a macroscopic flux

$$\phi\_{\rm diff}(a) = -\int\_{\hat{c}\Omega\_{t}} \mathbf{J}(a) . \mathbf{n} dA \tag{19}$$

The density of microscopic flux **J** *a* is associated to an invisible motion of molecules within a continuum description, hence must be described by a specific constitutive law. It does not depend on the velocity of the points of *<sup>t</sup>* .

The convective derivative along the vector field **w** of the field *aa t* **x**, writes

$$\frac{\partial \mathcal{L}\_w a}{\partial t} = \left(\frac{\partial a}{\partial t}\right)\_x + \nabla a.\mathbf{w} \tag{20}$$

In the case **w** coincides with the velocity of the material particles, previous relation delivers the definition of the material (or particular) derivative

$$\frac{da}{dt} \equiv \frac{\mathcal{S}\_v a}{\delta t} = \frac{\mathcal{S}\_w a}{\delta t} + \nabla a. (\mathbf{v} - \mathbf{w}) \tag{21}$$

The derivative of the volume integral : *t A adx* is next calculated, according to Leibniz rule:

$$\frac{D}{Dt} \int\_{\Omega\_t} ad\mathbf{x} = \int\_{\Omega\_t} \frac{\partial a}{\partial t} d\mathbf{x} + \int\_{\partial\Omega\_t} a (\mathbf{w}.\mathbf{n}) dA \tag{22}$$

with **w** the velocity field of the points on *<sup>t</sup>* , which is associated to a variation of the domain occupied by the material points of the growing solid body (Figure 1).

Fig. 1. Domain variation due to the virtual velocity field **w**

A global balance equation can next be written, according to the natural physical rule: the balance of any quantity is the sum of the production / destruction term and of the flux; this yields

( ) . *t*

(the minus accounts for the unit exterior normal **n** ). This diffusive flux corresponds to a

() . *t*

 *a a dA* 

The density of microscopic flux **J** *a* is associated to an invisible motion of molecules within a continuum description, hence must be described by a specific constitutive law. It

> . *<sup>w</sup> x*

 

In the case **w** coincides with the velocity of the material particles, previous relation delivers

. *v w da a a*

*D a adx dx a dA*

with **w** the velocity field of the points on *<sup>t</sup>* , which is associated to a variation of the

A global balance equation can next be written, according to the natural physical rule: the balance of any quantity is the sum of the production / destruction term and of the flux; this

 

> 

*t A adx* 

*tt t*

*Dt t* 

domain occupied by the material points of the growing solid body (Figure 1).

*a*

*a*

.

*a a*

*t t*

**v wn** (18)

**J n** (19)

**v w** (21)

**w n** (22)

is next calculated, according to Leibniz rule:

**w** (20)

 *a a dA* 

*conv* 

> *diff*

The convective derivative along the vector field **w** of the field *aa t* **x**, writes

*dt t t* 

does not depend on the velocity of the points of *<sup>t</sup>* .

the definition of the material (or particular) derivative

Fig. 1. Domain variation due to the virtual velocity field **w**

yields

The derivative of the volume integral :

macroscopic flux

$$\frac{D}{Dt} \int\_{\Omega\_t} ad\mathbf{x} = \int\_{\Omega\_t} \Pi d\mathbf{x} - \int\_{\partial\Omega\_t} \left\{ a(\mathbf{v} - \mathbf{w}) + \mathbf{J}(a) \right\} . \mathbf{n} dA \tag{23}$$

The first term on the r.h.s. corresponds to mass production, the second contribution to convection of the produced mass through the boundary *<sup>t</sup>* , and the third contribution to diffusion through the boundary of the moving volume *<sup>t</sup>* . One can see that only the relative velocity of particles w.r. to the surface velocity matters. Combining this identity with (22) gives

$$\int\_{\Omega\_t} \frac{\partial a}{\partial t} d\mathbf{x} + \int\_{\partial \Omega\_t} a(\mathbf{w}. \mathbf{n}) dA = \int\_{\Omega\_t} \Pi d\mathbf{x} - \int\_{\partial \Omega\_t} \left\{ a(\mathbf{v} - \mathbf{w}) + \mathbf{J}(a) \right\} . \mathbf{n} dA \tag{24}$$

The corresponding local balance law is obtained after elimination of the velocity **w** , hence

$$\frac{\delta\_v a}{\delta t} + a \text{div} \mathbf{v} = \Pi - \text{div} \mathbf{J}(a) \Leftrightarrow \frac{\partial a}{\partial t} + \text{div}(a\mathbf{v}) = \Pi - \text{div} \mathbf{J}(a) \tag{25}$$

**Mass balance**: the mass balance equation is deduced from the identification *a* , the actual density. Hence, (23) gives

$$\frac{D}{Dt} \int\_{\Omega\_{\rm{t}}} \rho \, d\mathbf{x} = \int\_{\Omega\_{\rm{t}}} \pi \, d\mathbf{x} - \int\_{\partial\Omega\_{\rm{t}}} \mathbf{J}(\rho) \cdot \mathbf{n} \, dA \tag{26}$$

The strong form of the balance law of mass writes finally

$$\frac{\delta \mathcal{S}\_{\boldsymbol{\nu}} \rho}{\delta t} = \boldsymbol{\Pi} - \operatorname{div} \mathbf{J}(\rho) - \rho \nabla \cdot \mathbf{v} \tag{27}$$

The **mass balance in Eulerian format** is given in terms of the actual density by the following reasoning: we first write the general form of the balance of mass in physical space as

$$\frac{D}{Dt} \int\_{\Omega\_{\!\!\!}} \rho \mathbf{dx} = \int\_{\Omega\_{\!\!\!}} \left( \frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{v} \right) d\mathbf{x} = \int\_{\Omega\_{\!\!\!}} \pi d\mathbf{x} + \int\_{\partial\Omega\_{\!\!\!}} m d\mathbf{s} \equiv \int\_{\Omega\_{\!\!\!}} \Gamma \rho d\mathbf{x} \tag{28}$$

with **x**,*t* the actual density, the physical source of mass, and *m* : . **m n** the scalar physical mass flux across the boundary, projection of the flux (vector) **m** . The previous balance law is quite general, as we account for both the variation of the integration volume through the term .**v** , and for the source and flux of mass reflected by the right hand side of (28). Localization of previous integral equation gives

$$\frac{D\rho}{Dt} = \pi + \nabla.\mathbf{m} - \rho \nabla.\mathbf{v} \tag{29}$$

with ( , ): *X t t* **x v x** the Eulerian velocity, which proves identical to (27); the same balance law has been obtained in (Epstein and Maugin, 2000) starting form its Lagrangian counterpart.

Thermodynamics of Surface Growth with Application to Bone Remodeling 377

*<sup>D</sup> dx dx d dx <sup>d</sup>*

with **σ** Cauchy stress and **f** body forces per unit physical volume. Localizing (32) gives

. *<sup>D</sup> div*

**Balance of kinetic and internal energy**: the first law of thermodynamics for an open system has to account for the contributions to kinetic and internal energies due to the incoming material. Denoting *u* the specific internal energy density, one may write the energy balance

2

**n σ vm v q x**

2 2

**v v fv v <sup>σ</sup> <sup>m</sup>**

*<sup>D</sup> dx div dx*

*<sup>D</sup> <sup>k</sup> div div*

<sup>1</sup> .. . 2 2

. : . . 2 2

 

**v v f v <sup>σ</sup> v n <sup>σ</sup> v m**

*t t*

 

2 2

*dx d*

*u d*

with *r* the volumetric heat supply (generated by growth), and **q** the heat flux across *<sup>t</sup>* . This writing of the energy balance can be simplified using the balance of kinetic energy with volumetric density *k* , obtained by multiplying (33) by the velocity and integrating over *<sup>t</sup>* ,

> 2 2

**v v fv v <sup>σ</sup> v m v fv v <sup>σ</sup> <sup>m</sup>**

2 2 2

**<sup>v</sup> v v m fv v <sup>σ</sup> <sup>m</sup>**

 

<sup>2</sup> . .. . <sup>2</sup> 2 2

*DK dx div dx*

Using again (30) delivers similarly the material derivative of the total energy (left-hand side

: . . .. . . . 2 2

The left hand side of previous equality can be expressed versus the material derivative of the total kinetic energy of the growing body, using the general equality (30) with <sup>1</sup> <sup>2</sup>

**v fv v**

1 1 . 2 2

*<sup>D</sup> u dx r u dx*

2 2

 

*t tt t t*

*Dt* 

<sup>1</sup> . . <sup>2</sup> *t t*

*t*

1

*t t*

*Dt*

2

**v**

*D*

*Dt Dt*

*t t*

in (34)) as (the total internal energy is denoted *U* )

*Dt*

*Dt*

*Dt*

using the mass balance (29)

in the actual configuration as

hence

hence (35)

 . .

**vf n <sup>σ</sup> v nm v** (32)

 

*t t*

**<sup>v</sup> <sup>f</sup> <sup>σ</sup> m v** (33)

(34)

2 *a* **v** ,

(35)

In the sequel, we shall extensively use the following expression of the material derivative of integrals of specific quantities (defined per unit mass) *aa t* **x**, , obtained using the mass balance (28)

$$\frac{D}{Dt} \int\_{\Omega\_{\rm t}} \rho \, d\mathbf{x} = \int\_{\Omega\_{\rm t}} \left\{ \rho \frac{Da}{Dt} + a \left( \boldsymbol{\pi} + \nabla \cdot \mathbf{m} \right) \right\} d\mathbf{x} \tag{30}$$

The comparison of (27) with (29) gives the identification of fluxes **J m** ; the balance law is further consistent with (and equivalent to) the writing (Ganghoffer and Haussy, 2005)

$$
\dot{\rho} + \rho \text{div}(\mathbf{v}) = \Phi\_{\rho} + \sigma\_{\rho},
$$

with . **m** the total flux of conduction and the volumetric source of mass. Observe the difference with the treatment of section 1 considering overall a closed system with no internal sources, reflected by equation (1.5): this first point of view considers the nutrients responsible for growth as part of the system, whereas they appear as external sources in the second viewpoint.

Expressing the total mass of the domain *t* as ( ) () *t m x dx <sup>t</sup>* , the mass variation due to

the transport phenomena is written as the following integral accounting for source terms, allowing the identification of the production term

$$\left(\frac{dm}{dt}\right)\_{source} \coloneqq \int\_{\Omega\_t} \pi d\mathbf{x} = \int\_{\Omega\_t} d\mathbf{x} \Longrightarrow \pi = \Gamma \rho^2$$

The time variation of chemical concentration of nutrients is due to exchange through the boundary accounted for by a flux . *t <sup>k</sup> ds* **<sup>j</sup> <sup>n</sup>** , and to a source term due to growth

$$\int\_{\Omega\_t} \Gamma \rho d\mathbf{x} = \int\_{\Omega\_t} Tr\left(\dot{\mathbf{F}}\_{\mathcal{S}} \mathbf{F}\_{\mathcal{S}}^{-1}\right) \rho d\mathbf{x} \text{ \textquotedbl{} hence \textquotedbl{}} \text{(see (28))}$$

$$\frac{D}{Dt} \int\_{\Omega\_{\ell}} \rho \mathbf{dx} = \int\_{\Omega\_{\ell}} \text{Tr} \left( \dot{\mathbf{F}}\_{\mathcal{S}} \mathbf{F}\_{\mathcal{S}}^{-1} \right) \rho \mathbf{dx} = \int\_{\Omega\_{\ell}} \rho \dot{m}\_{k} \mathbf{dx} + \int\_{\partial \Omega\_{\ell}} \dot{\mathbf{j}}\_{k} \cdot \mathbf{N} \mathbf{ds} \to \rho \dot{m}\_{k} + d \dot{v} \mathbf{j}\_{k} = \rho \left( \dot{\mathbf{F}}\_{\mathcal{S}} \mathbf{F}\_{\mathcal{S}}^{-1} : \mathbf{I} \right) = \rho \Gamma \tag{31}$$

The last equality is nothing else than . / **m** - a consequence of (28) - expressed in material format, with the identifications : *Tr***D***<sup>g</sup>* , ; **m j** *ii k n* . The global form *D D dx JdX Dt Dt* also fits within the general balance law for an open system, relation

#### (30) with 1 *a* .

*R R*

**Balance of momentum:** the Eulerian version of the balance of momentum writes (Epstein and Maugin, 2000)

In the sequel, we shall extensively use the following expression of the material derivative of integrals of specific quantities (defined per unit mass) *aa t* **x**, , obtained using the mass

*<sup>D</sup> Da adx a dx*

law is further consistent with (and equivalent to) the writing (Ganghoffer and Haussy,

*div*( )

 

Observe the difference with the treatment of section 1 considering overall a closed system with no internal sources, reflected by equation (1.5): this first point of view considers the nutrients responsible for growth as part of the system, whereas they appear as external

the transport phenomena is written as the following integral accounting for source terms,

The time variation of chemical concentration of nutrients is due to exchange through the

*<sup>k</sup> ds*

<sup>1</sup> <sup>1</sup> . . . :

*g g k k k k gg*

also fits within the general balance law for an open system, relation

**Balance of momentum:** the Eulerian version of the balance of momentum writes (Epstein

**F F j N j FF I** (31)

: *t t source dm dx dx*

*t*

. / **m**

 **v** 

*t t*

*Dt Dt* 

**m** the total flux of conduction and

Expressing the total mass of the domain *t* as ( ) ()

*dt*

, hence (see (28))

*<sup>D</sup> dx Tr dx n dx ds n div*

 

material format, with the identifications : *Tr***D***<sup>g</sup>* , ; **m j** *ii k*

*t t t t*

allowing the identification of the production term

boundary accounted for by a flux .

The last equality is nothing else than

<sup>1</sup> .

*g g*

 

*dx Tr dx*

The comparison of (27) with (29) gives the identification of fluxes **J m**

.

 

 

*t m x dx <sup>t</sup>* 

> 

**<sup>j</sup> <sup>n</sup>** , and to a source term due to growth

 


**<sup>m</sup>** (30)

, the mass variation due to

the volumetric source of mass.

; the balance

 

*n* . The global form

balance (28)

2005)

with . 

*t t*

**F F**

*R R D D dx JdX Dt Dt*

and Maugin, 2000)

*Dt*

(30) with 1 *a* .

sources in the second viewpoint.

$$\frac{D}{Dt} \int\_{\Omega\_{\natural}} \rho \mathbf{v} d\mathbf{x} = \int\_{\Omega\_{\natural}} \mathbf{f} d\mathbf{x} + \int\_{\partial \Omega\_{\natural}} \mathbf{n} \, d\sigma\_{\natural} + \int\_{\Omega\_{\natural}} \mathbf{n} \mathbf{v} d\mathbf{x} + \int\_{\partial \Omega\_{\natural}} \mathbf{n} (\mathbf{m} \otimes \mathbf{v}) d\sigma\_{\natural} \tag{32}$$

with **σ** Cauchy stress and **f** body forces per unit physical volume. Localizing (32) gives using the mass balance (29)

$$
\rho \frac{D\mathbf{v}}{Dt} = \mathbf{f} + \operatorname{div} \mathbf{o} + \left(\mathbf{m} \nabla\right) \mathbf{v} \tag{33}
$$

**Balance of kinetic and internal energy**: the first law of thermodynamics for an open system has to account for the contributions to kinetic and internal energies due to the incoming material. Denoting *u* the specific internal energy density, one may write the energy balance in the actual configuration as

$$\begin{split} \frac{D}{Dt} \int\_{\Omega\_{t}} \rho \left( u + \frac{1}{2} \mathbf{v}^{2} \right) d\mathbf{x} &= \int\_{\Omega\_{t}} \left\{ \left( \mathbf{f}, \mathbf{v} + r \right) + \pi \left( u + \frac{1}{2} \mathbf{v}^{2} \right) \right\} d\mathbf{x} \\ + \int\_{\partial \Omega\_{t}} \mathbf{n} \cdot \left\{ \mathbf{o}, \mathbf{v} + \mathbf{m} \left( u + \frac{1}{2} \mathbf{v}^{2} \right) - \mathbf{q} \right\} d\sigma(\mathbf{x}) \end{split} \tag{34}$$

with *r* the volumetric heat supply (generated by growth), and **q** the heat flux across *<sup>t</sup>* . This writing of the energy balance can be simplified using the balance of kinetic energy with volumetric density *k* , obtained by multiplying (33) by the velocity and integrating over *<sup>t</sup>* , hence

$$\begin{aligned} k &:= \frac{1}{2}\rho \frac{D\mathbf{v}^2}{Dt} = \mathbf{f}. \mathbf{v} + \mathbf{v}. div\mathbf{o} + \mathbf{v}. \left(\mathbf{m}. \nabla\right)\mathbf{v} \equiv \mathbf{f}. \mathbf{v} + \mathbf{v}. div\mathbf{o} + \mathbf{m}. \nabla\frac{\mathbf{v}^2}{2} \Rightarrow\\ \int\_{\Omega\_t} \frac{1}{2}\rho \frac{D\mathbf{v}^2}{Dt} dx &= \int\_{\Omega\_t} \left(\mathbf{f}. \mathbf{v} + \mathbf{v}. div\mathbf{o} + \mathbf{m}. \nabla\left(\frac{\mathbf{v}^2}{2}\right)\right) dx \end{aligned}$$

The left hand side of previous equality can be expressed versus the material derivative of the total kinetic energy of the growing body, using the general equality (30) with <sup>1</sup> <sup>2</sup> 2 *a* **v** , hence (35)

$$\begin{split} \frac{DK}{Dt} &= \int\_{\Omega\_{i}} \left\{ \rho \frac{D\left(\frac{\mathbf{v}^{2}}{2}\right)}{Dt} + \frac{\mathbf{v}^{2}}{2} \left(\boldsymbol{\pi} + \nabla.\mathbf{m}\right) \right\} d\mathbf{x} = \int\_{\Omega\_{i}} \left(\mathbf{f}.\mathbf{v} + \mathbf{v}.d\mathbf{i}v\mathbf{0} + \nabla.\left(\mathbf{m}\frac{\mathbf{v}^{2}}{2}\right) + \boldsymbol{\pi}\frac{\mathbf{v}^{2}}{2}\right) d\mathbf{x} = \\ & \int\_{\Omega\_{i}} \left\{ \mathbf{f}.\mathbf{v} - \mathbf{o} : \nabla\mathbf{v} + \boldsymbol{\pi}\frac{\mathbf{v}^{2}}{2} \right\} d\mathbf{x} + \int\_{\partial\Omega\_{i}} \mathbf{n}.\left\{ \mathbf{o}.\mathbf{v} + \mathbf{m}\frac{\mathbf{v}^{2}}{2} \right\} d\boldsymbol{\sigma} \end{split} \tag{35}$$

Using again (30) delivers similarly the material derivative of the total energy (left-hand side in (34)) as (the total internal energy is denoted *U* )

growth:

mass : *u s* 

Thermodynamics of Surface Growth with Application to Bone Remodeling 379

The transformation gradients , , **FFF** *a g* define the mappings of the tangent spaces to the various configurations. The Jacobean of the growth mapping informs about the nature of

: det( ) *g g J* **F** (3.3)

Hence 1 *gJ* describes growth, whereas 1 *gJ* represents resorption. Growth essentially

Adopting the framework of hyperelasticity, the first Piola-Kirchhoff stress **P** expresses from the strain energy density per unit volume in the reference configuration *W* **F X***<sup>a</sup>* ; with argument the reversible part of the transformation gradient (a possible explicit dependence

A more explicit (compared to (38)) expression of the dissipation accounting for heat and matter exchanges is obtained by considering the general form of the balance of energy and entropy: let denote *u* and *s* the density of internal energy and entropy per unit mass respectively; the first and second principles of thermodynamics write (Munster, 1970)

*u p <sup>q</sup> i kk* **J J F** ; . *s s*

production, always positive (it is dissipated). Introducing the free energy density per unit

 *s p* **J J JF**

*pi kk* **J** . . **F JJ** *<sup>q</sup>*

power of external and internal forces respectively), leads to the global form of previous

the boundary of . Previous inequality traduces the fact that the flux of mechanical work

*i i* , with *s* the entropy density, we then immediately obtain the rate of

*s* **J**

**P FX** : ; **F***W <sup>a</sup>* (3.4)

**<sup>J</sup> J J** the total entropy flux, *<sup>k</sup>* **<sup>J</sup>** the diffusion flux

*<sup>s</sup>* in previous inequality then expresses as

*i i* (3.7)

*dt* ( *<sup>K</sup>* is the kinetic energy, , *P P e i* being the virtual

(3.8)

**<sup>J</sup> <sup>n</sup>** respectively the flux of heat and mass through

. *e m <sup>k</sup> <sup>q</sup> <sup>k</sup>*

(3.5)

(3.6)

*<sup>s</sup>* the entropy

occurs between the referential and the actual configurations.

 .. 

with *<sup>q</sup>* **<sup>J</sup>** the heat diffusion flux, <sup>1</sup>

variation of the free energy density

The positivity of the entropy production

The principle of virtual power *e i*

**J n** and : . *m ii*

inequality in Eulerian format:

with : . *q q d*

upon the Lagrangian variable is included for heterogeneous media) as

: *<sup>s</sup> <sup>q</sup>*

of the k-specie, **F x** *<sup>k</sup>* ,*t* an external force acting on the k-specie, and

 . .. *<sup>q</sup> s kk i s* 

*dK P P*

*dK dx P J F dt* 

 *d*

$$\begin{split} \frac{D}{Dt}(\boldsymbol{L} + \boldsymbol{K}) &= \int\_{\Omega\_{t}} \rho \frac{D}{Dt} \Big( \boldsymbol{u} + \frac{1}{2} \mathbf{v}^{2} \Big) d\mathbf{x} + \int\_{\Omega\_{t}} \left( \boldsymbol{u} + \frac{1}{2} \mathbf{v}^{2} \right) (\boldsymbol{\pi} + \nabla \cdot \mathbf{m}) d\mathbf{x} = \int\_{\Omega\_{t}} \Big[ (\mathbf{f}, \mathbf{v} + \boldsymbol{r}) + \boldsymbol{\pi} \Big( \boldsymbol{u} + \frac{1}{2} \mathbf{v}^{2} \Big) \Big] d\mathbf{x} \\ &+ \int\_{\partial \Omega\_{t}} \mathbf{n} \cdot \Big[ \mathbf{o}, \mathbf{v} + \mathbf{m} \Big( \boldsymbol{u} + \frac{1}{2} \mathbf{v}^{2} \Big) - \mathbf{q} \Big] d\sigma(\mathbf{x}) \end{split}$$

Using the balance of kinetic energy (35) allows isolating the material derivative of the internal energy

$$\frac{DII}{Dt} = \int\_{\Omega\_t} (-\mathbf{o} : \nabla \mathbf{v} + r + \pi \mu) d\mathbf{x} + \int\_{\partial \Omega\_t} \mathbf{n} \, (\mathbf{m} \boldsymbol{\mu} - \mathbf{q}) d\sigma(\mathbf{x}) \, \mathbf{n}$$

Its strong form is given by localization using the general equality (30) with the identification *a u*

$$
\rho \frac{Du}{Dt} = -\mathbf{o} : \nabla \mathbf{v} + r + \mathbf{m} \nabla . \mu - \nabla . \mathbf{q} \tag{36}
$$

The Lagrangian counterpart of previous balance laws has been expressed in (Epstein and Maugin, 2000).

**Dissipation and second principle**: the dissipation inequality writes in global form as

$$\frac{D}{Dt} \int\_{\Omega\_t} \rho \mathbf{s} d\mathbf{x} \ge \int\_{\Omega\_t} (\pi \mathbf{s} + \theta^{-1} r) d\mathbf{x} - \int\_{\partial \Omega\_t} \mathbf{n} \cdot \frac{\mathbf{q}}{\theta} d\sigma(\mathbf{x}) \\ \Rightarrow \int\_{\Omega\_t} (\rho \frac{Ds}{Dt} + s \nabla \cdot \mathbf{m}) d\mathbf{x} \ge \int\_{\Omega\_t} (\theta^{-1} r) d\mathbf{x} - \int\_{\partial \Omega\_t} \mathbf{n} \cdot \frac{\mathbf{q}}{\theta} d\sigma(\mathbf{x}) \\ \Rightarrow \int\_{\Omega\_t} (\mathbf{n} \cdot \mathbf{p}) d\mathbf{x} - \int\_{\Omega\_t} \mathbf{n} \cdot \nabla \cdot \mathbf{q} d\sigma(\mathbf{x}) \\ \Rightarrow \int\_{\Omega\_t} (\mathbf{n} \cdot \mathbf{p}) d\mathbf{x} - \int\_{\Omega\_t} \mathbf{n} \cdot \nabla \cdot \mathbf{q} d\sigma(\mathbf{x})$$

Hence, the local dissipation inequality localizes as Clausius-Duhem inequality

$$
\rho \frac{Ds}{Dt} \ge \theta^{-1} r - \operatorname{div} \left( \frac{\mathbf{q}}{\theta} \right) - s \nabla. \mathbf{m} \tag{37}
$$

The previous balance laws are general balance laws in the framework of open systems irreversible thermodynamics; we shall in the next section make the fluxes and source terms involved in those balance laws more specific, in order to identify an evolution law for the volumetric growth of solid bodies.

## **3. Volumetric growth**

The kinematics of growth is elaborated from the classical multiplicative decomposition (Rodriguez et al., 1994) of the transformation gradient

$$\mathbf{F} = \nabla\_{\lambda} \mathbf{x}(\mathbf{X}, t) \to \mathbf{J} \; \coloneqq \det(\mathbf{F}) \tag{3.1}$$

with **X x**, the Lagrangian end Eulerian positions in the referential and actual configurations denoted , *R t* respectively, as the product of the growth deformation gradient **F***g* and the growth accommodation mapping **F***<sup>a</sup>*

$$\mathbf{F} = \mathbf{F}\_a.\mathbf{F}\_\mathcal{g}\tag{3.2}$$

*D D U K u dx u dx r u dx*

Using the balance of kinetic energy (35) allows isolating the material derivative of the

*DU r u dx u d x*

**<sup>σ</sup> <sup>v</sup> nm q**

Its strong form is given by localization using the general equality (30) with the identification

The Lagrangian counterpart of previous balance laws has been expressed in (Epstein and

**Dissipation and second principle**: the dissipation inequality writes in global form as

Hence, the local dissipation inequality localizes as Clausius-Duhem inequality

 

*Dt* 

*t t t t t t <sup>D</sup> Ds sdx s r dx d s dx r dx d*

 **q q nx m n x**

<sup>1</sup> . *Ds r div s*

The previous balance laws are general balance laws in the framework of open systems irreversible thermodynamics; we shall in the next section make the fluxes and source terms involved in those balance laws more specific, in order to identify an evolution law for the

The kinematics of growth is elaborated from the classical multiplicative decomposition

with **X x**, the Lagrangian end Eulerian positions in the referential and actual configurations denoted , *R t* respectively, as the product of the growth deformation gradient **F***g* and the

 

*t t*

: .

: ..

*r u*

1 1 .. .

2 2 2

**v v m fv v**

 

> 

**<sup>q</sup> <sup>m</sup>** (37)

**σ vm q** (36)

( , ) : det( ) *<sup>X</sup>* **F xX** *t J* **F** (3.1)

. **F FF** *<sup>a</sup> <sup>g</sup>* (3.2)

 

 

1 1 <sup>1</sup> . . 2 2 <sup>2</sup>

*tt t*

*Du*

*Dt* 

*Dt Dt*

(Rodriguez et al., 1994) of the transformation gradient

2

*Dt*

*u dx*

<sup>1</sup> . . <sup>2</sup>

**n σ vm v q**

*Dt Dt*

*t*

internal energy

Maugin, 2000).

 

volumetric growth of solid bodies.

growth accommodation mapping **F***<sup>a</sup>*

**3. Volumetric growth** 

*a u*

The transformation gradients , , **FFF** *a g* define the mappings of the tangent spaces to the various configurations. The Jacobean of the growth mapping informs about the nature of growth:

$$J\_{\mathcal{g}} \coloneqq \det(\mathbf{F}\_{\mathcal{g}}) \tag{3.3}$$

Hence 1 *gJ* describes growth, whereas 1 *gJ* represents resorption. Growth essentially occurs between the referential and the actual configurations.

Adopting the framework of hyperelasticity, the first Piola-Kirchhoff stress **P** expresses from the strain energy density per unit volume in the reference configuration *W* **F X***<sup>a</sup>* ; with argument the reversible part of the transformation gradient (a possible explicit dependence upon the Lagrangian variable is included for heterogeneous media) as

$$\mathbf{P} \coloneqq \partial\_{\mathbf{F}} \mathcal{W} \left( \mathbf{F}\_{a}; \mathbf{X} \right) \tag{3.4}$$

A more explicit (compared to (38)) expression of the dissipation accounting for heat and matter exchanges is obtained by considering the general form of the balance of energy and entropy: let denote *u* and *s* the density of internal energy and entropy per unit mass respectively; the first and second principles of thermodynamics write (Munster, 1970)

$$
\rho \dot{u} = -\nabla\_\cdot \mathbf{J}\_q - p\_i + \mathbf{J}\_k \, \mathbf{E}\_k \; ; \; \rho \dot{s} = -\nabla\_\cdot \mathbf{J}\_s + \sigma\_s \tag{3.5}
$$

with *<sup>q</sup>* **<sup>J</sup>** the heat diffusion flux, <sup>1</sup> : *<sup>s</sup> <sup>q</sup> i i* **<sup>J</sup> J J** the total entropy flux, *<sup>k</sup>* **<sup>J</sup>** the diffusion flux of the k-specie, **F x** *<sup>k</sup>* ,*t* an external force acting on the k-specie, and *<sup>s</sup>* the entropy production, always positive (it is dissipated). Introducing the free energy density per unit mass : *u s* , with *s* the entropy density, we then immediately obtain the rate of variation of the free energy density

$$\rho \dot{\boldsymbol{\varphi}} = \rho \mathbf{s} \dot{\boldsymbol{\theta}} - \nabla \cdot \mathbf{J}\_q + \theta \nabla \cdot \mathbf{J}\_s + \mathbf{J}\_k \, \mathbf{E}\_k - p\_i - \sigma\_s \tag{3.6}$$

The positivity of the entropy production *<sup>s</sup>* in previous inequality then expresses as

$$
\rho \rho \dot{\nu} \le -p\_i + \mathbf{J}\_k \mathbf{E}\_k + \nabla \cdot \left( \mathbf{J}\_q - \mu\_i \mathbf{J}\_i \right) \tag{3.7}
$$

The principle of virtual power *e i dK P P dt* ( *<sup>K</sup>* is the kinetic energy, , *P P e i* being the virtual power of external and internal forces respectively), leads to the global form of previous inequality in Eulerian format:

$$\frac{dK}{dt} + \int\_{\Omega} \rho \eta \dot{\nu} d\mathbf{x} \le P\_e + \underline{I}\_k \underline{F}\_k + \Phi\_m + \Phi\_q \tag{3.8}$$

with : . *q q d* **J n** and : . *m ii d* **<sup>J</sup> <sup>n</sup>** respectively the flux of heat and mass through the boundary of . Previous inequality traduces the fact that the flux of mechanical work

Thermodynamics of Surface Growth with Application to Bone Remodeling 381

*ii a k g a*

**J F I L F**

 

*a kg*

 **F I L**

 

(3.13)

**F** (3.14)

**L***<sup>g</sup> f* **Σ***<sup>a</sup>* (3.15)1

in the sequel), as pictured in figure 2. The

in the sequel.

(3.15)2

0 .: *<sup>t</sup>*

The dissipation splits into the sum of the thermal and chemical dissipation and the intrinsic

*Js Grad* **J** ; . :0 *<sup>t</sup>*

From (3.14), and as a generalization of the growth models initially written in a purely mechanical context, relations (3.3) and (3.4), one is entitled to write a general growth model

with the Eshelby stress accounting for both mechanical and chemical energy contributions

 **Σ F I F**

: . *<sup>t</sup> a a k a* 

Thereby, the Eshelby stress accounts for the change of domain induced by growth; this is further reflected in the material driving force for growth, including the (material) divergence of Eshehlby stress (Ganghoffer, 2010a, b). The exchange of matter is accounted for by the number of moles (with corresponding driving forces the chemical potentials), which may obey specific kinetic equations, of evolution diffusion type in a general setting

Simulations of volumetric growth based on this formalism have been done for academic situations in (Ganghoffer, 2010b). The objective of the present contribution is rather to unify volumetric and surface growth under a common umbrella, basing on the framework of

Surface thermodynamics is clearly a pluridisciplinary topic, which has its origins in the study of liquids, and touches various disciplines, such as metallurgy (grain boundary energy), fracture mechanics (fracture energy, mechanics (surface stress), physics of fluids (surface tension) and of solids (surface stress). Surface thermodynamic data are important parameters

The thermodynamics of surfaces has a long history, tracing back to Gibbs; an interface exists when a thin inhomogeneous element of material forms a transition zone separating two

The aim of this section is not to give a detailed account in each of those fields, but rather to provide the reader with a broad overview of the basic surface thermodynamics and to

Different viewpoints have been considered in the literature as to the geometry of the surface (this coinage used in Linford refers to the surface, as opposed to the bulk phases): the

for specialists in each of those fields, with however a different acceptance of the term.

 

transition zone between the bulk phases will be denoted by the Greek letter

review the major underlying parameters and their possible source of variation.

 

*a* 

*Js Grad J*

0 *i i*

   

(mechanical) dissipation

(Ganghoffer, 2010a, b).

Eshelby stress and material forces.

phases of different materials (denoted ,

**4. Surface growth: A review of the thermodynamics** 

according to

and mass increases the kinetic and internal free energy of the system, the difference being dissipated.

The second principle may be rewritten after a few manipulations in terms of a dynamical Eshelby stress accounting for all sources of energies (mechanical, chemical, thermal): the free energy density is taken to depend on the elastic part of the transformation gradient **F***<sup>a</sup>* , the concentration of chemical specie *nk* and the temperature , so that Clausius-Duhem inequality (3.7) becomes in material format:

$$\begin{split} & \frac{D}{Dt} (\rho \| \boldsymbol{\nu}) \leq \mathbf{T} : \dot{\mathbf{F}} + \nabla \left( \mathbf{J}\_{q} - \mu\_{i} \mathbf{J}\_{i} \right) \to \\ & \rho \| \rho \mathbf{I} : \dot{\mathbf{F}}\_{\mathcal{S}} \mathbf{F}\_{\mathcal{S}}^{-1} + \rho \mathbf{J} \frac{\partial \mathcal{\boldsymbol{\nu}} \boldsymbol{\nu}}{\partial \mathbf{F}\_{a}} \dot{\mathbf{F}}\_{a} + \rho \mathbf{J} \frac{\partial \mathcal{\boldsymbol{\nu}} \boldsymbol{\nu}}{\partial n\_{k}} \dot{n}\_{k} \\ & + \rho \mathbf{J} \frac{\partial \mathcal{\boldsymbol{\nu}} \boldsymbol{\nu}}{\partial \boldsymbol{\theta}} \dot{\boldsymbol{\theta}} \leq \mathbf{T} : \left( \dot{\mathbf{F}}\_{a} \mathbf{F}\_{\mathcal{S}} + \mathbf{F}\_{a} \dot{\mathbf{F}}\_{\mathcal{S}} \right) + \operatorname{Div} \mathbf{J}\_{q} - \mu\_{i} \operatorname{Div} \mathbf{J}\_{i} - \mathbf{J}\_{i} \operatorname{Grad} \mu\_{i} \end{split} \tag{3.9}$$

The balance of biochemical energy expresses that the time variation of chemical concentration of nutrients is due to exchange through the boundary accounted for by the term . *R dA* **J N** and to a source term due to growth <sup>1</sup> . *R R g g JdX Tr JdX* **F F** , hence

$$\begin{split} \frac{D}{Dt} \int\_{\Omega\_{\ell}} \rho d\mathbf{x} &= \int\_{\Omega\_{k}} Tr \Big( \dot{\mathbf{F}}\_{\mathcal{S}}, \mathbf{F}\_{\mathcal{S}}^{-1} \Big) \rho \mathbf{J} dX = \int\_{\Omega\_{k}} \rho f \dot{n}\_{k} dX \\ &+ \int\_{\partial \Omega\_{k}} \mathbf{J}\_{k} \mathbf{N} dA \to \rho f \dot{n}\_{k} + Div \mathbf{J}\_{k} = \rho f \Big( \dot{\mathbf{F}}\_{\mathcal{S}}, \mathbf{F}\_{\mathcal{S}}^{-1} : \mathbf{I} \Big) \rho \mathbf{J} \end{split} \tag{3.10}$$

The last equality is nothing else than . / **m** expressed in material format, identifying : *Tr***D***<sup>g</sup>* , **M J** *i i* . Accordingly, (3.9) becomes

$$\begin{split} 0 \le & \left( \mathbf{T} \mathbf{F}\_{\mathcal{S}}^{t} - \rho \mathbf{J} \frac{\partial \boldsymbol{\nu}}{\partial \mathbf{F}\_{a}} \right) \ddot{\mathbf{F}}\_{a} - \left( \frac{\partial \boldsymbol{\nu}}{\partial \mathbf{n}\_{k}} - \mu\_{k} \right) \rho \operatorname{J} \dot{\boldsymbol{n}}\_{k} - \mathbf{J}\_{i} \operatorname{Grad} \boldsymbol{\mu}\_{i} \\ & + \left( \mathbf{F}\_{a}^{t} \cdot \mathbf{T} \, \mathbf{F}\_{\mathcal{S}}^{t} - \rho \operatorname{J} \boldsymbol{\nu} \mathbf{I} - \rho \operatorname{J} \frac{\partial \boldsymbol{\nu}}{\partial \mathbf{n}\_{k}} \right) : \dot{\mathbf{F}}\_{\mathcal{S}} \, \mathbf{F}\_{\mathcal{S}}^{-1} + \rho \operatorname{J} \dot{\boldsymbol{s}} \dot{\boldsymbol{\sigma}} \end{split} \tag{3.11}$$

Since previous equality must hold true for arbitrary variations , **<sup>F</sup>***a k <sup>n</sup>* , the following constitutive equations for the first Piola-Kirchhoff stress and the chemical potential are obtained

$$\mathbf{T} = \rho \mathbf{J} \frac{\partial \boldsymbol{\varphi}}{\partial \mathbf{F}\_a} \mathbf{F}\_g^{-t}; \quad \mu\_k = \frac{\partial \boldsymbol{\varphi}}{\partial n\_k} \tag{3.12}$$

Especially, (3.12)1 is an alternative to (3.4) using the specific free energy instead of a strain energy potential; observe that is expressed per unit mass, in contrast to *W* **F X***<sup>a</sup>* ; in (3.4), expressed per unit referential volume.

The residual dissipation then writes from (3.11)

and mass increases the kinetic and internal free energy of the system, the difference being

The second principle may be rewritten after a few manipulations in terms of a dynamical Eshelby stress accounting for all sources of energies (mechanical, chemical, thermal): the free energy density is taken to depend on the elastic part of the transformation gradient **F***<sup>a</sup>* , the

*q ii*

 

**J N** and to a source term due to growth <sup>1</sup> .

*<sup>D</sup> dx Tr JdX Jn dX Dt*

.

*t R R*

**F F**

. . : .

**F TF I F F**

*J*

**T F**

*a g g g*

 

*a k*

**TF F <sup>J</sup> <sup>F</sup>**

identifying : *Tr***D***<sup>g</sup>* , **M J** *i i* . Accordingly, (3.9) becomes

*t*

*t t*

1

*g g k*

*dA Jn Div J J*

. . :

. / **m**

 

*g a k ki i*

*J J Js n*

*J Jn Grad n*

*k*

Since previous equality must hold true for arbitrary variations , **<sup>F</sup>***a k <sup>n</sup>* , the following constitutive equations for the first Piola-Kirchhoff stress and the chemical potential are

> . ; *<sup>t</sup> g k a k*

Especially, (3.12)1 is an alternative to (3.4) using the specific free energy instead of a strain

 

**J N J FF I**

*k k k gg*

*a k*

 

*J Div Div Grad*

The balance of biochemical energy expresses that the time variation of chemical concentration of nutrients is due to exchange through the boundary accounted for by the

**T FF FF J J J**

*n*

 

*a g a g q i ii i*

*R R*

1

*n*

 

1

 

> 

**F** (3.12)

is expressed per unit mass, in contrast to *W* **F X***<sup>a</sup>* ; in (3.4),

**F F**

, so that Clausius-Duhem

(3.9)

 

*g g*

*JdX Tr JdX*

 

expressed in material format,

, hence

(3.10)

(3.11)

concentration of chemical specie *nk* and the temperature

1

**IFF F**

: .

 

*R*

The last equality is nothing else than

0 .

 

: .

*J J Jn*

*g g a k*

**TF J J**

**F**

:. .

inequality (3.7) becomes in material format:

*<sup>D</sup> <sup>J</sup> Dt*

 

dissipated.

term . *R dA*

obtained

energy potential; observe that

expressed per unit referential volume.

The residual dissipation then writes from (3.11)

$$0 \le \left(\rho \text{Js}\dot{\boldsymbol{\theta}} - \mathbf{J}\_i \text{Grad}\boldsymbol{\mu}\_i\right) + \rho \mathbf{J}\left(\mathbf{F}\_a^\dagger \cdot \frac{\partial \boldsymbol{\nu}}{\partial \mathbf{F}\_a} - \left(\boldsymbol{\nu} + \boldsymbol{\mu}\_k\right) \mathbf{I}\right) \colon \mathbf{L}\_g\tag{3.13}$$

The dissipation splits into the sum of the thermal and chemical dissipation and the intrinsic (mechanical) dissipation

$$\{\rho\} \mathbf{s}\dot{\boldsymbol{\theta}} - \mathbf{J}\_i \mathbf{G} \boldsymbol{\pi} d\mu\_i \ge 0 : \left( \mathbf{F}\_a^\dagger, \frac{\partial \boldsymbol{\nu}}{\partial \mathbf{F}\_a} - \left( \boldsymbol{\nu} + \mu\_k \right) \mathbf{I} \right) \colon \mathbf{L}\_\mathcal{g} \ge 0 \tag{3.14}$$

From (3.14), and as a generalization of the growth models initially written in a purely mechanical context, relations (3.3) and (3.4), one is entitled to write a general growth model according to

$$\mathbf{L}\_{\mathcal{g}} = f\left(\tilde{\mathbf{L}}\_a\right) \tag{3.15}$$

with the Eshelby stress accounting for both mechanical and chemical energy contributions

$$\tilde{\boldsymbol{\Sigma}}\_{a} \coloneqq \rho \mathbf{F}\_{a}^{\;t} \cdot \frac{\hat{\boldsymbol{\mathcal{O}}} \boldsymbol{\nu}}{\hat{\boldsymbol{\mathcal{O}}} \mathbf{F}\_{a}} - \rho \left(\boldsymbol{\nu} + \mu\_{k}\right) \mathbf{I} \tag{3.15}$$

Thereby, the Eshelby stress accounts for the change of domain induced by growth; this is further reflected in the material driving force for growth, including the (material) divergence of Eshehlby stress (Ganghoffer, 2010a, b). The exchange of matter is accounted for by the number of moles (with corresponding driving forces the chemical potentials), which may obey specific kinetic equations, of evolution diffusion type in a general setting (Ganghoffer, 2010a, b).

Simulations of volumetric growth based on this formalism have been done for academic situations in (Ganghoffer, 2010b). The objective of the present contribution is rather to unify volumetric and surface growth under a common umbrella, basing on the framework of Eshelby stress and material forces.

### **4. Surface growth: A review of the thermodynamics**

Surface thermodynamics is clearly a pluridisciplinary topic, which has its origins in the study of liquids, and touches various disciplines, such as metallurgy (grain boundary energy), fracture mechanics (fracture energy, mechanics (surface stress), physics of fluids (surface tension) and of solids (surface stress). Surface thermodynamic data are important parameters for specialists in each of those fields, with however a different acceptance of the term.

The thermodynamics of surfaces has a long history, tracing back to Gibbs; an interface exists when a thin inhomogeneous element of material forms a transition zone separating two phases of different materials (denoted , in the sequel), as pictured in figure 2. The transition zone between the bulk phases will be denoted by the Greek letter in the sequel. The aim of this section is not to give a detailed account in each of those fields, but rather to provide the reader with a broad overview of the basic surface thermodynamics and to review the major underlying parameters and their possible source of variation.

Different viewpoints have been considered in the literature as to the geometry of the surface (this coinage used in Linford refers to the surface, as opposed to the bulk phases): the

Thermodynamics of Surface Growth with Application to Bone Remodeling 383

*k k dF S dT gdA dN*

The quantity *g* accounts for both the creation of new surface (with a fixed number of atoms) and the elastic deformation of the surface (also with a fixed number of atoms). The addition of atoms (particles) on the surface is accounted for by the last term. Considering

*dN dN dN dN* 

The superficial excess or molar superficial concentrations are then defined as : / *n NA*

*Z Z Z Z zV zV z A*

The purely elastic variation of the surface area, expressed by a superficial stress **σ** ,

*k k dU dS gdA dN*

*ZZ Z Z* 

*k k dU dS pdV gdA dN*

*k k dF Sd pdV gdA dN*

 

 

> 

, , *TVNi*

 ; *F ANk k* 

 

> 

 

. Combining both relations

*F <sup>g</sup> <sup>A</sup>*

> 

 

0 *S dT N d g dA Ad k k*

with *A* the area of the interface. Any extensive quantity *Z* can be decomposed as

surface (irreversible phenomena), with a constant number of particles.

*k k dF S d gdA dN*

 

 

the *superficial excess* quantity. Regarding surface quantities, one makes a

 

 

 

<sup>2</sup> *J m*/ - a scalar - accounting for the creation of a new

 

 

in-between, one can write the differential of

 ,

with a separating interface

 

For the whole system, using the previous decomposition

two phases ,

with : / *z ZA* 

one has

gives

distinction between:

 

The superficial energy

dual to an elastic surface strain **ε** . The excess total internal energy writes

The variation of the free energy is

This leads to the differential

Hence, *g* is defined as the partial derivative

 

the total number of particles as

surface phase is considered as two-dimensional by Gibbs, and coined the mathematical dividing surface, as the neat separation between fluid and solid phases. Gibbs viewpoint may be called the *surface excess approach (at fixed volume)*, in which the composite system (bulk phases and the interface) is the sum of the reference system without the interface and a correction; the difference of any quantity between the actual and the reference system leads to an *interfacial excess quantity*.

Fig. 2. Formation of an interface from a fixed number of moles of α and β.

Important to this viewpoint is the fact that the reference and actual systems have the same volume.

Guggenheim considered the surface phase as a three dimensional body of finite small thickness, and is commonly coined the *surface phase approach*. A third approach has been introduced by Goodrich, relying on Guggenheim vision, but with the interfaces between the surface phase and the two bulk phases identified to the walls of a confining vessel. A last vision at variant with Gibbs treatment advocates that both the actual and reference systems have the same mass, but possibly different volumes: it bears the name *Surface excess approach (at fixed mass)*, and was hardly considered in the literature, although rapidly mentioned by Gibbs in 1878. One drawback of the Guggenheim model is that the volume of the interfacial region *V* is arbitrary, and has nothing to do with the volume change that occurs during the formation of the interface; this difficulty is not apparent in Gibbs approach, for which the excess volume *V*is always zero.

For liquids, the situation is simple, as a single scalar parameter, the surface tension, is sufficient. Three parameters are required to characterize the thermodynamics of surfaces: the reversible work to produce unit area of new surface, sometimes called the *specific surface work* (the counterpart of the surface tension in liquids), the specific surface Helmholtz energy, as the change of energy of the surface region (as opposed to the bulk phases), and the surface stress tensor, defined as the reversible work required to produced a unit area of new surface by deformation. In order to avoid some existing confusion in the early literature (this is due to the oversimplified situation that prevails for liquids), those three parameters are next introduced in a distinct manner.

The thermodynamics of surfaces is based on the setting up of *excess quantities*. The reader is referred to (Linford, 1973) and (Couchman and Linford, 1980) for more details on the topic. Hence, the excess (Helmholtz) free energy is defined through its differential

$$dF^{\sigma} = -S^{\sigma}dT + gdA + \mu\_k dN^{\sigma}\_{\ \ \ k}$$

The quantity *g* accounts for both the creation of new surface (with a fixed number of atoms) and the elastic deformation of the surface (also with a fixed number of atoms). The addition of atoms (particles) on the surface is accounted for by the last term. Considering two phases , with a separating interface in-between, one can write the differential of the total number of particles as

$$d\mathcal{N} = d\mathcal{N}^a + d\mathcal{N}^\beta + dN^\sigma$$

The superficial excess or molar superficial concentrations are then defined as : / *n NA* , with *A* the area of the interface. Any extensive quantity *Z* can be decomposed as

$$\mathbf{Z} = \mathbf{Z}^{\alpha} + \mathbf{Z}^{\beta} + \mathbf{Z}^{\sigma} = z^{\alpha}V^{\alpha} + z^{\beta}V^{\beta} + z^{\sigma}A$$

with : / *z ZA* the *superficial excess* quantity. Regarding surface quantities, one makes a distinction between:


The excess total internal energy writes

$$d\mathcal{U}^{\sigma} = \theta dS^{\sigma} + \mathcal{g}dA + \mu\_k dN^{\sigma}{}\_k$$

For the whole system, using the previous decomposition

$$Z = Z^{\alpha} + Z^{\beta} + Z^{\sigma}$$

one has

382 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

surface phase is considered as two-dimensional by Gibbs, and coined the mathematical dividing surface, as the neat separation between fluid and solid phases. Gibbs viewpoint may be called the *surface excess approach (at fixed volume)*, in which the composite system (bulk phases and the interface) is the sum of the reference system without the interface and a correction; the difference of any quantity between the actual and the reference system leads

Fig. 2. Formation of an interface from a fixed number of moles of α and β.

Important to this viewpoint is the fact that the reference and actual systems have the same

Guggenheim considered the surface phase as a three dimensional body of finite small thickness, and is commonly coined the *surface phase approach*. A third approach has been introduced by Goodrich, relying on Guggenheim vision, but with the interfaces between the surface phase and the two bulk phases identified to the walls of a confining vessel. A last vision at variant with Gibbs treatment advocates that both the actual and reference systems have the same mass, but possibly different volumes: it bears the name *Surface excess approach (at fixed mass)*, and was hardly considered in the literature, although rapidly mentioned by Gibbs in 1878. One drawback of the Guggenheim model is that the volume of the interfacial

the formation of the interface; this difficulty is not apparent in Gibbs approach, for which

For liquids, the situation is simple, as a single scalar parameter, the surface tension, is sufficient. Three parameters are required to characterize the thermodynamics of surfaces: the reversible work to produce unit area of new surface, sometimes called the *specific surface work* (the counterpart of the surface tension in liquids), the specific surface Helmholtz energy, as the change of energy of the surface region (as opposed to the bulk phases), and the surface stress tensor, defined as the reversible work required to produced a unit area of new surface by deformation. In order to avoid some existing confusion in the early literature (this is due to the oversimplified situation that prevails for liquids), those three parameters

The thermodynamics of surfaces is based on the setting up of *excess quantities*. The reader is referred to (Linford, 1973) and (Couchman and Linford, 1980) for more details on the topic.

Hence, the excess (Helmholtz) free energy is defined through its differential

is arbitrary, and has nothing to do with the volume change that occurs during

to an *interfacial excess quantity*.

volume.

region *V*

the excess volume *V*

are next introduced in a distinct manner.

is always zero.

$$d\mathcal{L}I = \theta dS - p dV + \mathcal{g} dA + \mu\_k dN\_k$$

The variation of the free energy is

$$dF = -S d\theta - p dV + \mathcal{g} dA + \mu\_k dN\_k$$

Hence, *g* is defined as the partial derivative , , *TVNi F <sup>g</sup> <sup>A</sup>* . Combining both relations

$$dF^{\sigma} = -S^{\sigma}d\theta + \mathcal{g}dA + \mu\_k dN^{\sigma}{}\_k ; F^{\sigma} = \mathcal{\gamma}A + \mu\_k N^{\sigma}{}\_k$$

gives

$$S^{\sigma}dT + N^{\sigma} \, \_k d\mu\_k - \left(g - \gamma\right) dA + A d\gamma = 0$$

This leads to the differential

Thermodynamics of Surface Growth with Application to Bone Remodeling 385

; the variation of free

from the bulk, and

*<sup>i</sup>* . Imagine a

, at fixed

*<sup>i</sup>* are

, and modify

be adjusted to maintain the chemical potentials at fixed specified values

 

returned to their initial values, but with the surface area increased by *dA*

Subtracting both variations of Helmholtz free energy by unit surface gives

, the specific surface entropy *s S S dA* :' /

But one can also express the specific surface Helmholtz energy as *i i f*

by

surface area; *<sup>i</sup> dn* particles from the bulk will enter the solid phase

thereafter the temperature by *dT* , and the ith chemical potential by *<sup>i</sup> d*

 

ithechange of Helmholtz free energy *F*

energy of this system of larger area is

*f* :' / *dF dF dA*

and thus finally

tensor, denoted

The second order tensor *ij*

induced by the component

*ij* . It is related to

excess *n d N N dA i ii* :' / , ,

the entropy of the phase

with *S*

> 

 

> 

 

 

leads to

modification of the temperature by *dT* , and of the ith chemical potential by *<sup>i</sup> d*

will be

*i ii* , , *i i dF S dT d dN S dT dN*

, '' ' *i i dF S dT dN*

, , ' ' ' *i i*

 

*dF dF S S dN N dT dA dA dA*

Introducing therein the definitions of the specific surface Helmholtz energy

*i i df s dT d*

*ii ii i i df d d d s dT d*

*i i d s dT d* 

The same identity was derived by (Goodrich, 1969) for a one-component system using the method of Lagrange multipliers. The reversible work needed to generate a unit area of new surface by stretching at constant pressure and temperature represents the surface stress

*ij ij*

 

 

 

 

  *i*

 

*ij* 

is the strain (a small perturbation scheme is presently adopted)

*ij* acting in the jth direction per unit length of the edge normal

 

> 

 

 , and the surface

, hence

. Consider next a new system for which *T* and

$$d\gamma = -\mathbf{s}^{\sigma}d\theta - \mathbf{n}^{\sigma}d\mu\_k + (\mathcal{g} - \mathcal{\gamma})d\mathbf{A}/A$$

expressing the variation of the superficial energy. In the case of an isothermal surface stretch with a constant chemical potential, one gets the Couchmann-Everett formula

$$\mathbf{g} = \mathcal{Y} + A \left( \frac{\partial \mathcal{Y}}{\partial \mathcal{A}} \right)\_{\mathcal{T}, \mu\_0}$$

In the case of a purely elastic stretch, previous formula specializes to the relation *el*

, *i T g A <sup>A</sup>* .

The reversible work needed to form a unit area of new surface is defined at constant temperature and pressure as the partial derivative of the Gibbs free energy of the entire system (bulk phases and surface), quantity *G PTn A* ,, ,*<sup>i</sup>* , with respect to the formed area *A* , at constant temperature *T* , pressure *P* , and number of moles of each component *ni* , viz

$$dG = -VdP - SdT + \mu\_i dN\_i + \gamma dA\_i$$

whereby the multiplicative factors of the differential elements on the right-hand side of *dG* are the partial derivatives

$$V = -\frac{\partial G}{\partial P};\ \mathbf{S} = -\frac{\partial G}{\partial T};\ \mu\_i = \frac{\partial G}{\partial N\_i};\ \mathcal{V} = \frac{\partial G}{\partial A}$$

The partial derivatives are evaluated with all three other variables being held fixed. The specific surface work includes two contributions, the change of Gibbs free energy per unit area for the surface region, denoted *g* , and the change per unit area of surface created from the surrounding bulk phases, evaluated as the sum *i i n* over all components,

$$\boldsymbol{\gamma} = \mathbf{g}^{\sigma} - \mu\_i \mathbf{n}\_i^{\sigma}$$

with : / *n NA i i* , the surface excess of the ith species. In terms of the Helmholtz energy of the whole system *F* , one has the similar relation involving the Helmholtz free energy per unit area *f* , viz

$$\gamma = f^{\sigma} + Pe - \mu\_i n\_i^{\sigma}$$

with *e* the thickness of the surface; in most cases, the parameter *e* is small, and one may neglect the contribution *Pe* , hence one has the identification *f g* . The last two formulas are expressions of Gibb's adsorption equation, with the derivation due to Mullins, which is next reproduced. We consider a system with n components consisting of a solid phase in contact with a fluid phase and a solid phase acting as a thermal bath at temperature *T* and as a chemical reservoir for each component; accordingly, the components concentrations can be adjusted to maintain the chemical potentials at fixed specified values *<sup>i</sup>* . Imagine a modification of the temperature by *dT* , and of the ith chemical potential by *<sup>i</sup> d* , at fixed surface area; *<sup>i</sup> dn* particles from the bulk will enter the solid phase from the bulk, and ithechange of Helmholtz free energy *F*will be

$$dF^{\sigma} = -S^{\sigma}dT + \left(\mu\_i + d\mu\_i\right)dN\_{i,\sigma} \approx -S^{\sigma}dT + \mu\_i dN\_{i,\sigma}$$

with *S* the entropy of the phase . Consider next a new system for which *T* and *<sup>i</sup>* are returned to their initial values, but with the surface area increased by *dA* , and modify thereafter the temperature by *dT* , and the ith chemical potential by *<sup>i</sup> d* ; the variation of free energy of this system of larger area is

$$dF^{\prime \sigma} \approx -S^{\prime \sigma} \, dT + \mu\_i dN^{\prime}\_{\,i,\sigma}$$

Subtracting both variations of Helmholtz free energy by unit surface gives

$$\frac{dF^{\sigma\sigma} - dF^{\sigma}}{dA^{\sigma}} \approx -\frac{\left(S^{\sigma\sigma} - S^{\sigma}\right)}{dA^{\sigma}} dT + \mu\_i \frac{d\left(N^{\prime}\_{\ i, \sigma} - N\_{i, \sigma}\right)}{dA^{\sigma}}$$

Introducing therein the definitions of the specific surface Helmholtz energy *f* :' / *dF dF dA* , the specific surface entropy *s S S dA* :' / , and the surface excess *n d N N dA i ii* :' / , , leads to

$$df^{\sigma} = -s^{\sigma}dT + \mu\_i d\Gamma\_i$$

But one can also express the specific surface Helmholtz energy as *i i f* , hence

$$df^{\sigma} = d\underline{\gamma} + \mu\_i d\Gamma\_i + \Gamma\_i d\mu\_i = -\mathbf{s}^{\sigma} dT + \mu\_i d\Gamma\_i$$

and thus finally

384 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

*d s d n d g dA A k k* /

expressing the variation of the superficial energy. In the case of an isothermal surface stretch

*g A <sup>A</sup>*

In the case of a purely elastic stretch, previous formula specializes to the relation

The reversible work needed to form a unit area of new surface is defined at constant temperature and pressure as the partial derivative of the Gibbs free energy of the entire system (bulk phases and surface), quantity *G PTn A* ,, ,*<sup>i</sup>* , with respect to the formed area *A* , at constant temperature *T* , pressure *P* , and number of moles of each component *ni* ,

> *i i dG VdP SdT dN dA*

> > ; S ; ; *<sup>i</sup>*

*i i g n* 

*i i f Pe n*

 

the whole system *F* , one has the similar relation involving the Helmholtz free energy per

with *e* the thickness of the surface; in most cases, the parameter *e* is small, and one may

are expressions of Gibb's adsorption equation, with the derivation due to Mullins, which is next reproduced. We consider a system with n components consisting of a solid phase

contact with a fluid phase and a solid phase acting as a thermal bath at temperature *T* and as a chemical reservoir for each component; accordingly, the components concentrations can

neglect the contribution *Pe* , hence one has the identification *f g*

 

, the surface excess of the ith species. In terms of the Helmholtz energy of

> >

. The last two formulas

in

*G G GG <sup>V</sup> P TNA* 

The partial derivatives are evaluated with all three other variables being held fixed. The

whereby the multiplicative factors of the differential elements on the right-hand side of *dG*

 

, *T <sup>i</sup>*

 

*i*

 

includes two contributions, the change of Gibbs free energy per unit

*n* 

, and the change per unit area of surface created

over all components,

 

with a constant chemical potential, one gets the Couchmann-Everett formula

, *i*

*T*

*el*

.

are the partial derivatives

specific surface work

with : / *n NA i i* 

> , viz

unit area *f*

from the surrounding bulk phases, evaluated as the sum *i i*

area for the surface region, denoted *g*

 

*g A <sup>A</sup>*

viz

$$d\gamma = -s^{\sigma}dT - \Gamma\_i d\mu\_i$$

The same identity was derived by (Goodrich, 1969) for a one-component system using the method of Lagrange multipliers. The reversible work needed to generate a unit area of new surface by stretching at constant pressure and temperature represents the surface stress tensor, denoted *ij* . It is related to by

$$
\tilde{\sigma}\_{ij} = \mathcal{\mathcal{Y}}\delta\_{ij} + \frac{\partial \mathcal{Y}}{\partial \tilde{\mathcal{E}}\_{ij}}
$$

The second order tensor *ij* is the strain (a small perturbation scheme is presently adopted) induced by the component *ij* acting in the jth direction per unit length of the edge normal

Thermodynamics of Surface Growth with Application to Bone Remodeling 387

Combining the stretch and shear processes then lead to the expression of the surface stress

*ij ij*

represented by a 2 by 2 symmetrical matrix (3 independent components). For an isotropic material or a crystal with a threefold (or greater axis of symmetry), it follows as shown by

 

Shuttleworth (using the principle of virtual work) the isotropic surface stress

Lastly, consider a square section in the xy plane of side

of surface energy, which expresses for a high symmetry isotropic crystal as

 /

thermodynamic quantities of the whole system are due to the new /

*excess values* of the corresponding quantities, denoted with a subscript

 

. Therefore, one has

1/2

1/2

*xx A*

*yy A*

due to the equality *A xx yy dA* 

overall additional energy contribution.

along a plane to generate the

of molecules)

 

*g*

, with an expanded work given by

<sup>2</sup> 1 *WA A yy xx yy* 

Assuming the deformation is reversible and isothermal, the total work spent is the variation

1 2 1 *xx xx yy xx yy dA W W A A gdA*

*<sup>d</sup> gdA dA A d g dA*

Note that the last term vanishes for liquids; as a corollary, liquid films can easily be stretched since atoms can more from the bulk to the surface without additional energy costs. The opposite situation prevails for solids, as they shear and their structure changes with an

The Gibbs approach towards interfacial excess quantities is as previously mentioned valid only at fixed volume; a parallel approach that is valid at fixed mass instead has been

initially separated and interface-free, and are in a thought experiment imagined to be joined

them from Gibbs approach at fixed volume. The differential of the Gibbs energy of the system before and after formation of the interface successively writes (for a constant number

developed in (Muller and Kern, 2001), which is next exposed. The bulk phases ,

 

*ij d d* 

1 0 0 1 *ij*

> 1/2 *A*

*WA A* <sup>1</sup> *xx xx* 

 

 

interface. Since mass is conserved, any change in the

 

 

> are

interface, coined

to distinguish

 

and imagine an extension

 

. Extend

; the required work is 1/2 1/2

1/2 1/2

 

 

tensor

of the x edge by

next the y edge by

to the ith direction, with both indices i, j lying in the plane of the surface. Previous equation is valid for an anisotropic solid, and reduces in the case of an isotropic surface to the previously written Couchmann-Everett formula, with *g* half the trace of the surface stress tensor. The proof of previous formula follows Mullins derivation: let imagine a unit cube with edges parallel to the axes , , *x y z* , and perform two distinct operations on it:


$$\mathcal{W}\_0 + \mathcal{Z} \{ \boldsymbol{\gamma} + d\boldsymbol{\gamma} \} (1 + d\mathbf{x}) = \mathcal{Z}\boldsymbol{\gamma} + \mathcal{W}\_1$$

The difference *W W* 1 0 is the work due to the stretching operations of both processes, and can be equalized to the x-component *xx* of a force in the newly formed surface times the distance 2*dx* through which this force acts, hence

$$2\tilde{\sigma}\_{\infty}d\mathbf{x} = \mathcal{W}\_1 - \mathcal{W}\_0$$

The strain *xx dx* (since the other side has unit length), hence

$$2\gamma + 2\gamma \Delta \varepsilon\_{\infty} + 2\Delta \gamma + 2\Delta \gamma \Delta \varepsilon\_{\infty} = 2\gamma + 2\tilde{\sigma}\_{\infty} \Delta \varepsilon\_{\infty}$$

Due to the equalities

$$
\Delta \mathfrak{P} \mathfrak{A} \varepsilon\_{\infty} \approx 0 \; ; \; \Delta \chi \; / \; \Delta \varepsilon\_{\infty} \approx d \gamma \; / \; d \varepsilon\_{\infty}.
$$

it finally results

$$
\tilde{\sigma}\_{\infty} = \gamma + \frac{d\gamma}{d\varepsilon\_{\infty}}.
$$

Similar analogous processes with the stretching replaced by shear lead to the relation (Linford, 1973)

$$
\tilde{\sigma}\_{xy} = \frac{d\mathcal{Y}}{d\varepsilon\_{xy}}
$$

to the ith direction, with both indices i, j lying in the plane of the surface. Previous equation is valid for an anisotropic solid, and reduces in the case of an isotropic surface to the previously written Couchmann-Everett formula, with *g* half the trace of the surface stress tensor. The proof of previous formula follows Mullins derivation: let imagine a unit cube

i. Stretch the cube reversibly along the x axis by an amount *dx* , with the y edge fixed, but allowing the edge z to vary its length. The surface in the xz plane may then change by an inflow (or outflow) of material from the bulk, increasing (or decreasing), the cube height; denote *W*0 the work expanded in this transformation. Let next separate the

ii. Separate the original unit cube into two parts along an xy plane, requiring the work

two surfaces are created, and the factor (1 ) *dx* since the specific surface work applies

edge. Let *W*1 be the work expanded in this operation. The final configuration is the same as that obtained in the first process, hence the same total work has been expanded,

*W d dx W* 0 1 2 (1 ) 2

The difference *W W* 1 0 is the work due to the stretching operations of both processes, and

1 0 2 *xx*

22 2 2 22 *xx xx xx xx*

 

 ; / / *xx xx* 

*dx W W*

 

*xx*

 

Similar analogous processes with the stretching replaced by shear lead to the relation

*xy*

*dx* (since the other side has unit length), hence

, and stretch each half by *dx* in the x direction, at fixed y edge, but varying z

 

> 

 *d d*

*xx d d* 

*xy d d* 

 

*xx* of a force in the newly formed surface times the

due to the stretch *dx* (factor 2 arises since

, with *d*

stretched cube along the xy plane, requiring the work *W d dx* <sup>2</sup> 2 (1 )

the variation of the specific surface work

per unit surface area).

hence *WWWW* <sup>0213</sup> , viz

can be equalized to the x-component

distance 2*dx* through which this force acts, hence

0 *xx* 

<sup>3</sup> *W* 2

The strain *xx* 

Due to the equalities

it finally results

(Linford, 1973)

with edges parallel to the axes , , *x y z* , and perform two distinct operations on it:

Combining the stretch and shear processes then lead to the expression of the surface stress tensor

$$
\tilde{\sigma}\_{ij} = \mathcal{\eta}\delta\_{ij} + \frac{d\mathcal{\eta}}{d\varepsilon\_{ij}}.
$$

represented by a 2 by 2 symmetrical matrix (3 independent components). For an isotropic material or a crystal with a threefold (or greater axis of symmetry), it follows as shown by Shuttleworth (using the principle of virtual work) the isotropic surface stress

$$
\tilde{\sigma}\_{ij} = \mathcal{g} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}
$$

Lastly, consider a square section in the xy plane of side 1/2 *A* and imagine an extension of the x edge by 1/2 *xx A* ; the required work is 1/2 1/2 *WA A* <sup>1</sup> *xx xx* . Extend next the y edge by 1/2 *yy A* , with an expanded work given by

$$\mathcal{W}\_2 = \tilde{\sigma}\_{yy} \left( A^{\sigma} \right)^{1/2} \left( 1 + \varepsilon\_{\infty} \right) \varepsilon\_{yy} \left( A^{\sigma} \right)^{1/2}$$

Assuming the deformation is reversible and isothermal, the total work spent is the variation of surface energy, which expresses for a high symmetry isotropic crystal as

$$d\left(A^{\sigma}\mathcal{Y}\right) = \mathcal{W}\_1 + \mathcal{W}\_2 = \tilde{\sigma}\_{xx}\varepsilon\_{xx}A^{\sigma} + \tilde{\sigma}\_{yy}\left(1 + \varepsilon\_{xx}\right)\varepsilon\_{yy}A^{\sigma} \approx gdA^{\sigma}$$

due to the equality *A xx yy dA* . Therefore, one has

$$
\mathcal{g}dA^{\sigma} = \mathcal{\mathcal{Y}}dA^{\sigma} + A^{\sigma}d\mathcal{Y} \Rightarrow \mathcal{g} = \mathcal{Y} + \frac{d\mathcal{Y}}{dA^{\sigma}}
$$

Note that the last term vanishes for liquids; as a corollary, liquid films can easily be stretched since atoms can more from the bulk to the surface without additional energy costs. The opposite situation prevails for solids, as they shear and their structure changes with an overall additional energy contribution.

The Gibbs approach towards interfacial excess quantities is as previously mentioned valid only at fixed volume; a parallel approach that is valid at fixed mass instead has been developed in (Muller and Kern, 2001), which is next exposed. The bulk phases , are initially separated and interface-free, and are in a thought experiment imagined to be joined along a plane to generate the / interface. Since mass is conserved, any change in the thermodynamic quantities of the whole system are due to the new / interface, coined *excess values* of the corresponding quantities, denoted with a subscript to distinguish them from Gibbs approach at fixed volume. The differential of the Gibbs energy of the system before and after formation of the interface successively writes (for a constant number of molecules)

Thermodynamics of Surface Growth with Application to Bone Remodeling 389

The specific interfacial excess energy is obtained by simply integrating the differential

can be determined from those of \* *G*

\* \*\* \* *<sup>n</sup> <sup>U</sup> TS PV*

, , *<sup>n</sup> PAn*

identified as the surface energy, which is the sum of the interfacial tension and the

 

*T*

heat supplied by the surrounding for an isothermal creation of new interface. The advantages of this last approach in comparison to Gibbs treatment is that it leads to non-nil interfacial volumes, analogues of the Maxwell relations for bulk phases can be derived, and the temperature and pressure dependence of the interfacial tension can be accessed from a

orientation dependent, although not for a liquid. This surface energy parameter has been up to now considered under the thermodynamic continuum viewpoint; we next examine two other viewpoints, the atomistic approach and Wulff plot. The atomistic approach considers the interaction between atoms to calculate the surface energy; arrangement of atoms in crystals are such that one can order atoms according to the energy required to remove atoms from the bulk: first nearest neighbours requiring more energy compared to second and third nearest neighbours. For a crystal lattice presenting dislocations, the number of broken bonds

\* \* \* \*\*

   

 

excess energy is obtained from a Legendre transform to \* *dG V dP S dT dA*

at constant pressure and temperature, introducing

: / . Last relation implies that the temperature

,

 

, is for a solid generally

. The specific interfacial

and

 

 ,

substitution of the previous interfacial Maxwell relations, thus

It immediately results the specific interfacial excess enthalpy

 ,

comparison between simple formulae and experiments.

Fig. 3. The broken bond model for surface energy

The reversible work to form new surface area, parameter

 

*<sup>n</sup> dG V dP S dT dA G*

 

the specific interfacial excess energy \* *G GA*

 

and pressure dependence of \*

with \* *H* \* \*\*

   ,

> 

*H TS T*

$$\mathbf{d} \cdot \mathbf{G}\_1 = \begin{pmatrix} V\_{\alpha} + V\_{\beta} \\ \end{pmatrix} \mathbf{d}P - \begin{pmatrix} S\_{\alpha} + S\_{\beta} \\ \end{pmatrix} \mathbf{d}T \;/; \mathbf{d} \mathbf{G}\_2 = \begin{pmatrix} V\_{\alpha} + V\_{\beta} + V\_{\gamma} \\ \end{pmatrix} \mathbf{d}P - \begin{pmatrix} S\_{\alpha} + S\_{\beta} + S\_{\gamma} \\ \end{pmatrix} \mathbf{d}T + \gamma^\* \mathbf{d}A$$

with \* the *reorganization surface energy*, although commonly referred to as the interfacial tension in the literature; it is a mechanical positive quantity, that may depend upon interface curvature. Note that the number of atoms is the same in the reference and final states, in contrast with Gibbs approach. Hence, the variation of the excess Gibbs free energy between states 1 and 2 for the fixed masses , *m m* is

$$d\mathbf{G}\_{\boldsymbol{\gamma}} = d\mathbf{G}\_2 - d\mathbf{G}\_1 = V\_{\boldsymbol{\gamma}}dP - \mathbf{S}\_{\boldsymbol{\gamma}}dT + \boldsymbol{\gamma}^\*d\mathbf{A}$$

which may be interpreted from an energetic point of view as follows: the term *V dP* is the mechanical work done against the external force field, the contribution *S dT* represents the heat of formation of the interface, and \* *dA* is the mechanical work done against the internal force field of both phases , by motion of the molecules from the bulk to generate a new interface. The *excess free energy* of formation of the interface, potential *G* , is the additional free energy required to form the interface from fixed masses of the preexisting bulk phases , . The above equations implicitly use the conservation of mass, equation

$$n\_{\text{total}} = n\_a + n\_\beta$$

and the definition of the excess interfacial volume *V* from the contributions to the total volume after interface formation (balance law for the volume)

$$V\_{total} = V\_{a} + V\_{\beta} + V\_{\gamma}$$

In contrast to this treatment, Gibbs assumes a conservation of the total volume as *V VV total* , but with addition of the new mass *n*such that

$$m\_{\text{total}} = n\_{\alpha} + n\_{\beta} + n^{\sigma}$$

As a compensation for the volume change accompanying the formation of the interface; hence, *n* is a supply of material from outside the system, with the sense that the Gibbs volume is an open thermodynamic volume.

Due to its status as a state function, the previous differential OF *G* allows writing relations between partial derivatives as the analogues of the bulk phase Maxwell relations

$$\begin{split} \left(\frac{\partial \boldsymbol{\hat{\gamma}}^{\star}}{\partial \boldsymbol{T}}\right)\_{\boldsymbol{P},\boldsymbol{A},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} &= -\left(\frac{\partial \boldsymbol{S}}{\partial \boldsymbol{A}}\right)\_{\boldsymbol{P},\boldsymbol{T},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} = -\boldsymbol{S}\_{\boldsymbol{\gamma}}^{\*} \; : \; \left(\frac{\partial \boldsymbol{\gamma}^{\star}}{\partial \boldsymbol{P}}\right)\_{\boldsymbol{T},\boldsymbol{A},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} = \left(\frac{\partial \boldsymbol{V}}{\partial \boldsymbol{A}}\right)\_{\boldsymbol{P},\boldsymbol{T},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} &= \boldsymbol{V}\_{\boldsymbol{\gamma}} \cdot \boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}} \\ & \quad \left(\frac{\partial \boldsymbol{V}\_{\boldsymbol{\gamma}}}{\partial \boldsymbol{T}}\right)\_{\boldsymbol{P},\boldsymbol{A},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} = -\left(\frac{\partial \boldsymbol{S}}{\partial \boldsymbol{P}}\right)\_{\boldsymbol{A},\boldsymbol{T},\boldsymbol{n}\_{\boldsymbol{a},\boldsymbol{\beta}}} \end{split}$$

The introduced quantities \* *S V*, are respectively the *interfacial excess entropy* and the *specific interfacial excess volume*; the compact notation , *<sup>n</sup>* stands for the two quantities *n n* , .

tension in the literature; it is a mechanical positive quantity, that may depend upon interface curvature. Note that the number of atoms is the same in the reference and final states, in contrast with Gibbs approach. Hence, the variation of the excess Gibbs free energy between

2 1 *dG dG dG V dP S dT dA*

which may be interpreted from an energetic point of view as follows: the term *V dP*

generate a new interface. The *excess free energy* of formation of the interface, potential *G*

the additional free energy required to form the interface from fixed masses of the pre-

*n nn total* 

*V VVV total* 

In contrast to this treatment, Gibbs assumes a conservation of the total volume as

*n nnn total*

As a compensation for the volume change accompanying the formation of the interface;

*V S T P* 

;

, , , , *PAn* , , *AT n*

 

 

 

 is

mechanical work done against the external force field, the contribution *S dT*

 

  

the *reorganization surface energy*, although commonly referred to as the interfacial

; \*

<sup>2</sup> *dG V V V dP S S S dT dA*

\*

. The above equations implicitly use the conservation of mass,

such that

\* \*

*P A* 

> 

, ,

stands for the two quantities *n n*

*<sup>V</sup> <sup>V</sup>*

 

*T An* , , *PTn* , ,

are respectively the *interfacial excess entropy* and the *specific* 

is a supply of material from outside the system, with the sense that the Gibbs

from the contributions to the total

allows writing relations

.

 , .

*dA* is the mechanical work done against the

by motion of the molecules from the bulk to

represents the

is the

, is

*dG V V dP S S dT* <sup>1</sup> 

states 1 and 2 for the fixed masses , *m m*

heat of formation of the interface, and \*

 

volume is an open thermodynamic volume.

and the definition of the excess interfacial volume *V*

volume after interface formation (balance law for the volume)

, but with addition of the new mass *n*

Due to its status as a state function, the previous differential OF *G*

, ,

, , , , *PTn PAn*

 

*T A*

*interfacial excess volume*; the compact notation , *<sup>n</sup>*

 

The introduced quantities \* *S V*,

\* \*

between partial derivatives as the analogues of the bulk phase Maxwell relations

*<sup>S</sup> <sup>S</sup>*

 

internal force field of both phases ,

existing bulk phases ,

equation

*V VV total* 

hence, *n*

with \* 

  The specific interfacial excess energy is obtained by simply integrating the differential , \* \*\* *<sup>n</sup> dG V dP S dT dA G* at constant pressure and temperature, introducing the specific interfacial excess energy \* *G GA* : / . Last relation implies that the temperature and pressure dependence of \* can be determined from those of \* *G* . The specific interfacial excess energy is obtained from a Legendre transform to \* *dG V dP S dT dA* and substitution of the previous interfacial Maxwell relations, thus

$$\left(\left(\boldsymbol{\mathcal{U}}\_{\boldsymbol{\mathcal{V}}}\right)\_{\boldsymbol{n}\_{\boldsymbol{\alpha},\boldsymbol{\beta}}} = \boldsymbol{\mathcal{V}}^{\*} + \boldsymbol{T}\boldsymbol{S}\_{\boldsymbol{\mathcal{V}}}^{\*} - \boldsymbol{P}\boldsymbol{V}\_{\boldsymbol{\mathcal{V}}}^{\*}\right)$$

It immediately results the specific interfacial excess enthalpy

$$\left(H\_{\mathcal{V}}^{\*}\right)\_{n\_{a,\beta}} = \mathcal{Y}^{\*} + TS\_{\mathcal{V}}^{\*} \equiv \mathcal{Y}^{\*} - T\left(\frac{\partial \mathcal{Y}^{\*}}{\partial T}\right)\_{P,A,n\_{a,\beta}}$$

with \* *H* identified as the surface energy, which is the sum of the interfacial tension and the heat supplied by the surrounding for an isothermal creation of new interface. The advantages of this last approach in comparison to Gibbs treatment is that it leads to non-nil interfacial volumes, analogues of the Maxwell relations for bulk phases can be derived, and the temperature and pressure dependence of the interfacial tension can be accessed from a comparison between simple formulae and experiments.

Fig. 3. The broken bond model for surface energy

The reversible work to form new surface area, parameter , is for a solid generally orientation dependent, although not for a liquid. This surface energy parameter has been up to now considered under the thermodynamic continuum viewpoint; we next examine two other viewpoints, the atomistic approach and Wulff plot. The atomistic approach considers the interaction between atoms to calculate the surface energy; arrangement of atoms in crystals are such that one can order atoms according to the energy required to remove atoms from the bulk: first nearest neighbours requiring more energy compared to second and third nearest neighbours. For a crystal lattice presenting dislocations, the number of broken bonds

Thermodynamics of Surface Growth with Application to Bone Remodeling 391

The present model aims at describing radial bone remodeling, accounting for chemical and mechanical influences from the surrounding. Our approach of bone growth typically follows the streamlines of continuum mechanical models of bone adaptation, including the time-dependent description of the external geometry of cortical bone surfaces in the spirit of free boundary value problems – a process sometimes called net 'surface remodeling' - and of the bone material properties, sometimes coined net 'internal remodeling' (Cowin, 2001).

In the sequel, the framework for surface growth elaborated in (Ganghoffer, 2010) will be applied to describe bone modeling and remodeling. As a prerequisite, we recall the identification of the driving forces for surface growth. We consider a tissue element under grow 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

Focusing on the surface behavior, the potential energy of the growing tissue element is set as

Thereby, the growing solid surface is supposed to be endowed with a volumetric density *W*<sup>0</sup> **F** depending upon the transformation gradient : **F x** *<sup>X</sup>* , a surface energy with

the unit normal vector **N** to *Sg* , and possibly explicitly upon the surface position vector

*g S g k k g S S*

 

**FNX***<sup>S</sup>* per unit reference surface, depending upon the surface gradient **F**

 

(5.1)

<sup>0</sup> , ;

**F FNX**

*S g g g g g*

*d p dl p dl*

*g g g*

*S*

*V W dx d nd*

...

**fx x τ x ν**

*gg g*

*SS S*

 

,

**n**

**5. Model of surface growth with application to bone remodeling** 

with anisotropic surface tension

the expression

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

**5.1 Material driving forces for surface growth** 

normal to the surface *Sg* (forces acting in the tangent plane).

Fig. 5. Tissue element under growth: elements of differential geometry.

is direction dependent, and is given by the expressions cos / *a* and sin / *a* in the x and y directions respectively (figure 3), with *a* the distance to nearest neighbour (function of the type of atomic packing) and the inclination of the overall crystal shape resulting from the total number of steps being created.

The surface energy is given by the expression

$$E\_{surf} = \left(\cos\theta + \sin\theta\right)\varepsilon\_b \not\begin{pmatrix} 2a^2 \end{pmatrix}$$

with *<sup>b</sup>* the energy per bond. The broken bond model can be used to determine the shape of a small crystal from the minimization of the sum of surface energies *<sup>i</sup>* over all crystal faces, a concept introduced in 1878 by J. W. Gibbs, considering constant pressure, volume, temperature and molar mass:

$$M \text{im} \sum\_{i} A\_{i} \mathcal{V}\_{i}$$

at constant energy, hence adding the constraint 0 *i i i dE dA* . The dependence of on orientation of the crystal's surface and its equilibrium shape are condensed into a *Wulff plot*; in 1901, George Wulff stated that the length of a vector normal to a crystal face is proportional to its surface energy in this orientation. This is known as the Gibbs-Wulff theorem, which was initially given without proof, and was proven in 1953 by Conyers Herring, who at the same time provided a two steps method to determine the equilibrium shape of a crystal: in a first step, a polar plot of the surface energy as a function of orientation is made, given as the so-called gamma plot denoted as **n** , with **n** the normal to the surface corresponding to a particular crystal face. The second step is Wulff construction, in which the gamma plot determines graphically which crystal faces will be present: Wulff construction of the equilibrium shape consists in drawing a plane through each point on the γ-plot perpendicular to the line connecting that point to the origin. The inner envelope of all planes is geometrically similar to the equilibrium shape (figure 4).

Fig. 4. Wulff's construction to calculate the minimizing surface for a fixed volume

with anisotropic surface tension **n**

390 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

directions respectively (figure 3), with *a* the distance to nearest neighbour (function of the

 <sup>2</sup> cos sin / 2 *E a surf* 

a concept introduced in 1878 by J. W. Gibbs, considering constant pressure, volume,

orientation of the crystal's surface and its equilibrium shape are condensed into a *Wulff plot*; in 1901, George Wulff stated that the length of a vector normal to a crystal face is proportional to its surface energy in this orientation. This is known as the Gibbs-Wulff theorem, which was initially given without proof, and was proven in 1953 by Conyers Herring, who at the same time provided a two steps method to determine the equilibrium shape of a crystal: in a first step, a polar plot of the surface energy as a function of

to the surface corresponding to a particular crystal face. The second step is Wulff construction, in which the gamma plot determines graphically which crystal faces will be present: Wulff construction of the equilibrium shape consists in drawing a plane through each point on the γ-plot perpendicular to the line connecting that point to the origin. The inner envelope of all planes is geometrically similar to the equilibrium shape (figure 4).

Fig. 4. Wulff's construction to calculate the minimizing surface for a fixed volume

 *b*

the energy per bond. The broken bond model can be used to determine the shape of

*i i i Min A*

> *i dE dA*

the inclination of the overall crystal shape resulting from the

*a* and sin /

. The dependence of

*a* in the x and y

over all crystal faces,

**n** , with **n** the normal

on

is direction dependent, and is given by the expressions cos /

a small crystal from the minimization of the sum of surface energies *<sup>i</sup>*

at constant energy, hence adding the constraint 0 *i i*

orientation is made, given as the so-called gamma plot denoted as

type of atomic packing) and

temperature and molar mass:

with *<sup>b</sup>* 

total number of steps being created.

The surface energy is given by the expression
