Thermodynamics of Surface Growth with Application to Bone Remodeling

1. Introduction

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) 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$${\left(A.B\right)}_{ij}=A_{ik}{B}_{kj}$$. 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 $$a=a(x,t)$$ satisfies a prototype balance law of the form

$$\frac{\partial a(x,t)}{\partial t}=-\nabla .J_a(x,t)+{\sigma }_a(x,t)$$E1

with $$J_a(x,t)$$ the flux density of $$a(x,t)$$and $${\sigma }_a(x,t)$$ the local production (or destruction) of$$a(x,t)$$. 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$${\rho }_k$$, obey the local balance law (Vidal et al., 1994)

$$\frac{\partial {\rho }_k}{\partial t}=-\nabla .\left({\rho }_{k}u+J_k\right)+M_k{\displaystyle \sum _{\alpha =\mathrm{1..}r}{\nu }_{\alpha k}{J}_{\alpha }}$$E2

with $$u:=\frac{1}{\rho }{\displaystyle \sum _{k=1}^n{\rho }_k{u}_k}$$ the local barycentric velocity, $$M_k$$the molar mass, and $${\nu }_{\alpha k}$$ the stoechiometric coefficients in the reaction$$\alpha$$, such that the variation of the mass $$d{m}_k$$ of the species k due to chemical reactions expresses as

$$d{m}_k=M_k{\displaystyle \sum _{\alpha =\mathrm{1..}r}{\nu }_{\alpha k}{\xi }_{\alpha }},\text{k=1}\mathrm{..}\text{n}$$E3

wherein $${\xi }_{\alpha }$$ denotes the degree of advancement of reaction$$\alpha$$. The molar masses $$M_k$$ satisfy the global conservation law (due to Lavoisier)

$${\displaystyle \sum _{k=1}^n{\nu }_{\alpha k}}{M}_k=\mathrm{0,}\alpha =\mathrm{1..}r$$E4

Observe that the total flux of mass is the sum of a convective flux $${\rho }_{k}u$$ and a diffusive flux$$J_k$$ ; the mass production is identified as the contribution$$M_k{\displaystyle \sum _{\alpha =\mathrm{1..}r}{\nu }_{\alpha k}{J}_{\alpha }}$$. In this viewpoint, the system is in fact closed, since the balance law satisfied by the global density $$\rho ={\displaystyle \sum _{k=1}^n{\rho }_k}$$ writes (Vidal et al., 1994) accounting for the relation$${\displaystyle \sum _{j=1}^n{J}_k={\displaystyle \sum _{j=1}^n{\rho }_k{u}_k=0}}$$, as

$$\frac{\partial \rho }{\partial t}=-\nabla .\left(\rho u\right)+{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{\alpha =\mathrm{1..}r}{M}_k{\nu }_{\alpha k}{J}_{\alpha }}\equiv -\nabla .\left(\rho u\right)}$$E5

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, $$n_k=m_k/M_k$$, with $$m_k$$ the mass of the same constituent. The molar concentration is defined as$$c_k=n_k/V$$, its inverse being called the partial molar volume. The partial mole number $$n_k$$ satisfies the balance equation

$$\frac{\partial n_k}{\partial t}=-div{J}_k+\frac{{\partial }_i{n}_k}{\partial t}$$E6

with $$J_k$$ the flux of species k and $$\frac{{\partial }_i{n}_k}{\partial t}$$ its production term, given by De Donder definition of the rate of progress of the j^th chemical reaction

$$\frac{{\partial }_i{n}_k}{\partial t}={\displaystyle \sum _{j=1}^r{\nu }_{kj}{\dot{\xi }}_j}$$E7

The two previous equalities enter into Gibbs relation as

$$\rho \dot{u}=\theta \rho {\dot{s}}_e+\theta \rho {\dot{s}}_i+\sigma :\dot{\epsilon }-{\displaystyle \sum _k{\mu }_k\frac{\rho }{M}div{J}_k}+{\displaystyle \sum _k{\mu }_k\frac{\rho }{M}{\displaystyle \sum _j{\nu }_{kj}{\dot{\xi }}_j}}$$E8

with $$\theta$$ the temperature and $${\mu }_k$$ the chemical potential of constituent k. The chemical affinity in the sense of De Donder is defined as the force conjugated to the rate $${\dot{\xi }}_j$$

$$A_j=-{\displaystyle \sum _k\left({\mu }_k\frac{\rho }{M}\right){\nu }_{kj}}=-{\displaystyle \sum _k\left({\mu }_k/V\right){\nu }_{kj}}$$E9

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

$$\rho \dot{s}=\frac{1}{\theta }\rho \dot{u}-\frac{1}{\theta }\sigma :\dot{\epsilon }+{\displaystyle \sum _k\frac{1}{\theta }\left({\mu }_k\frac{1}{V}\right)div{J}_k}+{\displaystyle \sum _j\frac{1}{\theta }{A}_j{\dot{\xi }}_j}$$E10

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

$$\rho \dot{u}=-\nabla .J_q+\rho \dot{w}$$E11

with $$J_q$$ the heat flux, and the term $$\rho 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{s}=\rho {\dot{s}}_e+\rho {\dot{s}}_i=\frac{1}{\theta }\nabla .J_q+{\displaystyle \sum _k\nabla .\left(J_k\frac{{\mu }_k}{\theta }\frac{1}{V}\right)}-{\displaystyle \sum _k{J}_k.\nabla \left(\frac{{\mu }_k}{\theta }\frac{1}{V}\right)}+\frac{1}{\theta }\left(\rho \dot{w}-\sigma :\dot{\epsilon }\right)+{\displaystyle \sum _j\frac{1}{\theta }{A}_j{\dot{\xi }}_j}$$E12

The contribution $$-\sigma :\dot{\epsilon }/\theta$$ (involving the virtual power of internal forces) is further decomposed into

$$\frac{1}{\theta }\sigma :\dot{\epsilon }=\frac{1}{\theta }\sigma :\nabla u=\nabla .\left(\frac{1}{\theta }\sigma :u\right)-u.\nabla .\left(\frac{1}{\theta }\sigma \right)$$E13

Hence, the rate of the entropy density decomposes into

$$\begin{array}{l}\rho \dot{s}=\rho {\dot{s}}_e+\rho {\dot{s}}_i=-\nabla .\left(J_q+{\displaystyle \sum _k{J}_k}\frac{{\mu }_k}{\theta }\frac{1}{V}\right)+J_q.\nabla \left(\frac{1}{\theta }\right)-\\ {\displaystyle \sum _k{J}_k.\nabla \left(\frac{{\mu }_k}{\theta }\frac{1}{V}\right)}+\frac{1}{\theta }\left(\rho \dot{w}-\sigma :\dot{\epsilon }\right)+{\displaystyle \sum _j\frac{1}{\theta }{A}_j{\dot{\xi }}_j}\end{array}$$E14

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

$$J_s=J_q+{\displaystyle \sum _k{J}_k}\frac{{\mu }_k}{\theta }\frac{1}{V}$$E15

and of the internal entropy production

$$\rho {\dot{s}}_i=J_q.\nabla \left(\frac{1}{\theta }\right)-{\displaystyle \sum _k{J}_k.\nabla \left(\frac{{\mu }_k}{\theta }\frac{1}{V}\right)}-\frac{1}{\theta }\left(\sigma :\dot{\epsilon }-\rho \dot{w}\right)+{\displaystyle \sum _j\frac{1}{\theta }{A}_j{\dot{\xi }}_j}$$E16

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{u}=-\nabla .J_q+\sigma :\dot{\epsilon }+{\displaystyle \sum _k\left({\mu }_k/V\right){\dot{n}}_k}$$E17

One has assumed in this alternative that the mechanical power $$\rho \dot{w}={\sigma }^{eq}:\nabla u$$ does not include a flux contribution, hence only the heat diffusion contributes to the flux of internal energy. The contribution $$\sigma :\dot{\epsilon }+{\displaystyle \sum _k\left({\mu }_k/V\right){\dot{n}}_k}$$ is identified to the term$$\rho \dot{w}$$. Previous equality combined with the second principle, equality $$\rho \dot{s}=-\nabla .\left(\frac{J_q}{\theta }\right)+\rho {\dot{s}}_i$$ (the entropy flux resumes to the sole heat flux), delivers after a few manipulations the variation of the internal energy as

$$\rho \dot{u}=\left(\theta \rho \dot{s}-J_q\frac{\nabla \theta }{\theta }-T\rho {\dot{s}}_i\right)+\sigma :\dot{\epsilon }+{\displaystyle \sum _k\left({\mu }_k/V\right){\dot{n}}_k}$$E18

Hence, the internal entropy production is identified as

$$\theta \rho {\dot{s}}_i=-\rho \left(\dot{u}-T\dot{s}\right)-J_q\frac{\nabla \theta }{\theta }+\sigma :\dot{\epsilon }+{\displaystyle \sum _k\left({\mu }_k/V\right){\dot{n}}_k}$$E19

which is conveniently rewritten in terms of Helmholtz free energy density $$\psi :=u-Ts$$ as

$$\theta \rho {\dot{s}}_i=-\rho \left(\dot{\psi }+s\dot{\theta }\right)-J_q.\frac{\nabla \theta }{\theta }+\sigma :\dot{\epsilon }+{\displaystyle \sum _k\left({\mu }_k/V\right){\dot{n}}_k}$$E20

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 $$\pi$$ being produced (by a set of generating cells) at each point within the time varying volume$${\Omega }_t$$; 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(a)=av$$; the overall convective flux of a across the closed surface $$\partial {\Omega }_t$$ expresses then as

$${\varphi }_{conv}(a)=-{\displaystyle \underset{\partial {\Omega }_t}{\int }a\left(v-w\right).ndA}$$E21

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

$${\varphi }_{diff}(a)=-{\displaystyle \underset{\partial {\Omega }_t}{\int }J\left(a\right).ndA}$$E22

The density of microscopic flux $$J\left(a\right)$$ 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$$\partial {\Omega }_t$$.

The convective derivative along the vector field $$w$$of the field $$a=a\left(x,t\right)$$ writes

$$\frac{{\delta }_{w}a}{\delta t}={\left(\frac{\partial a}{\partial t}\right)}_x+\nabla a.w$$E23

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{{\delta }_{v}a}{\delta t}=\frac{{\delta }_{w}a}{\delta t}+\nabla a.\left(v-w\right)$$E24

The derivative of the volume integral $$A:={\displaystyle \underset{{\Omega }_t}{\int }adx}$$is next calculated, according to Leibniz rule:

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }adx}={\displaystyle \underset{{\Omega }_t}{\int }\frac{\partial a}{\partial t}dx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }a\left(w.n\right)dA}$$E25

with $$w$$ the velocity field of the points on$$\partial {\Omega }_t$$, which is associated to a variation of the domain occupied by the material points of the growing solid body (Figure 1).

Figure 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

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }adx}={\displaystyle \underset{{\Omega }_t}{\int }\Pi dx}-{\displaystyle \underset{\partial {\Omega }_t}{\int }\left\{a\left(v-w\right)+J\left(a\right)\right\}.ndA}$$E26

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

$${\displaystyle \underset{{\Omega }_t}{\int }\frac{\partial a}{\partial t}dx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }a\left(w.n\right)dA}={\displaystyle \underset{{\Omega }_t}{\int }\Pi dx}-{\displaystyle \underset{\partial {\Omega }_t}{\int }\left\{a\left(v-w\right)+J\left(a\right)\right\}.ndA}$$E27

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

$$\frac{{\delta }_{v}a}{\delta t}+adivv=\Pi -divJ\left(a\right)\iff \frac{\partial a}{\partial t}+div\left(av\right)=\Pi -divJ\left(a\right)$$E28

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

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho dx}={\displaystyle \underset{{\Omega }_t}{\int }\pi dx}-{\displaystyle \underset{\partial {\Omega }_t}{\int }J\left(\rho \right).ndA}$$E29

The strong form of the balance law of mass writes finally

$$\frac{{\delta }_v\rho }{\delta t}=\Pi -divJ\left(\rho \right)-\rho \nabla .v$$E30

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

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho dx}={\displaystyle \underset{{\Omega }_t}{\int }\left(\frac{D\rho }{Dt}+\rho \nabla .v\right)}dx={\displaystyle \underset{{\Omega }_t}{\int }\pi dx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }mds}\equiv {\displaystyle \underset{{\Omega }_t}{\int }\Gamma \rho dx}$$E31

with $$\rho \left(x,t\right)$$ the actual density, $$\pi$$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$$\rho \nabla .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 .m-\rho \nabla .v$$E32

with $$v(x,t):={\left(\frac{\partial x}{\partial t}\right)}_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.

In the sequel, we shall extensively use the following expression of the material derivative of integrals of specific quantities (defined per unit mass)$$a=a\left(x,t\right)$$, obtained using the mass balance (28)

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho adx=}{\displaystyle \underset{{\Omega }_t}{\int }\left\{\rho \frac{Da}{Dt}+a\left(\pi +\nabla .m\right)\right\}dx}$$E33

The comparison of (27) with (29) gives the identification of fluxes$$J\left(\rho \right)\equiv -m$$; the balance law is further consistent with (and equivalent to) the writing (Ganghoffer and Haussy, 2005)

$$\dot{\rho }+\rho div(v)={\Phi }_{\rho }+{\sigma }_{\rho }$$E34

with $${\Phi }_{\rho }\equiv \nabla .m$$ the total flux of conduction and $${\sigma }_{\rho }\equiv \pi$$ 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 $${\Omega }_t$$ as$$m({\Omega }_t)={\displaystyle \underset{{\Omega }_t}{\int }\rho (x)dx}$$, 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}:={\displaystyle \underset{{\Omega }_t}{\int }\pi dx}={\displaystyle \underset{{\Omega }_t}{\int }dx}\Rightarrow \pi =\Gamma \rho$$E35

The time variation of chemical concentration of nutrients is due to exchange through the boundary accounted for by a flux$${\displaystyle \underset{\partial {\Omega }_t}{\int }{j}_k.nds}$$, and to a source term due to growth$${\displaystyle \underset{{\Omega }_t}{\int }\Gamma \rho dx}={\displaystyle \underset{{\Omega }_t}{\int }Tr\left({\dot{F}}_g.F_g{}^{-1}\right)\rho dx}$$, hence (see (28))

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho dx}={\displaystyle \underset{{\Omega }_t}{\int }Tr\left({\dot{F}}_g.F_g{}^{-1}\right)\rho dx}={\displaystyle \underset{{\Omega }_t}{\int }\rho {\dot{n}}_{k}dx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }{j}_k.Nds}\to \rho {\dot{n}}_k+div{j}_k=\rho \left({\dot{F}}_g.F_g{}^{-1}:I\right)=\rho \Gamma$$E36

The last equality is nothing else than $$\Gamma =\left(\pi +\nabla .m\right)/\rho$$ - a consequence of (28) - expressed in material format, with the identifications$$\Gamma :=Tr\left(D_g\right)$$,$$m_i=j_i;\pi =\rho {\dot{n}}_k$$. The global form $$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_R}{\int }\rho dx}=\frac{D}{Dt}{\displaystyle \underset{{\Omega }_R}{\int }\rho JdX}$$ also fits within the general balance law for an open system, relation (30) with$$a\equiv 1$$.

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

$$\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho vdx}={\displaystyle \underset{{\Omega }_t}{\int }fdx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\sigma d{\sigma }_t}+{\displaystyle \underset{{\Omega }_t}{\int }\pi vdx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\left(m\otimes v\right)d{\sigma }_t}$$E37

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

$$\rho \frac{Dv}{Dt}=f+div\sigma +\left(m.\nabla \right)v$$E38

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{array}{l}\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho \left(u+\frac{1}{2}{v}^2\right)}dx={\displaystyle \underset{{\Omega }_t}{\int }\left\{\left(f.v+r\right)+\pi \left(u+\frac{1}{2}{v}^2\right)\right\}}dx\\ +{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\left\{\sigma .v+m\left(u+\frac{1}{2}{v}^2\right)-q\right\}}d\sigma \left(x\right)\end{array}$$E39

with $$r$$ the volumetric heat supply (generated by growth), and $$q$$ the heat flux across$$\partial {\Omega }_t$$. 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$${\Omega }_t$$, hence

$$\begin{array}{l}k:=\frac{1}{2}\rho \frac{D{v}^2}{Dt}=f.v+v.div\sigma +v.\left(m.\nabla \right)v\equiv f.v+v.div\sigma +m.\nabla \frac{v^2}{2}\Rightarrow \\ {\displaystyle \underset{{\Omega }_t}{\int }\frac{1}{2}\rho \frac{D{v}^2}{Dt}dx=}{\displaystyle \underset{{\Omega }_t}{\int }\left(f.v+v.div\sigma +m.\nabla \left(\frac{v^2}{2}\right)\right)dx}\end{array}$$E40

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$$a=\frac{1}{2}{v}^2$$, hence (35)

$$\begin{array}{l}\frac{DK}{Dt}={\displaystyle \underset{{\Omega }_t}{\int }\left\{\rho \frac{D\left({\frac{v}{2}}^2\right)}{Dt}+{\frac{v}{2}}^2\left(\pi +\nabla .m\right)\right\}dx}={\displaystyle \underset{{\Omega }_t}{\int }\left(f.v+v.div\sigma +\nabla .\left(m\frac{v^2}{2}\right)+\pi {\frac{v}{2}}^2\right)dx}=\\ {\displaystyle \underset{{\Omega }_t}{\int }\left(f.v-\sigma :\nabla v+\pi {\frac{v}{2}}^2\right)dx}+{\displaystyle \underset{\partial {\Omega }_t}{\int }n.}\left\{\sigma .v+m\frac{v^2}{2}\right\}d\sigma \end{array}$$E41

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$$)

$$\begin{array}{l}\frac{D}{Dt}\left(U+K\right)={\displaystyle \underset{{\Omega }_t}{\int }\rho \frac{D}{Dt}\left(u+\frac{1}{2}{v}^2\right)dx}+{\displaystyle \underset{{\Omega }_t}{\int }\left(u+\frac{1}{2}{v}^2\right)\left(\pi +\nabla .m\right)dx}={\displaystyle \underset{{\Omega }_t}{\int }\left\{\left(f.v+r\right)+\pi \left(u+\frac{1}{2}{v}^2\right)\right\}}dx\\ +{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\left\{\sigma .v+m\left(u+\frac{1}{2}{v}^2\right)-q\right\}}d\sigma \left(x\right)\end{array}$$E42

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

$$\frac{DU}{Dt}={\displaystyle \underset{{\Omega }_t}{\int }\left(-\sigma :\nabla v+r+\pi u\right)}dx+{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\left(mu-q\right)}d\sigma \left(x\right)$$E43

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

$$\rho \frac{Du}{Dt}=-\sigma :\nabla v+r+m\nabla .u-\nabla .q$$E44

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}{\displaystyle \underset{{\Omega }_t}{\int }\rho s}dx\ge {\displaystyle \underset{{\Omega }_t}{\int }\left(\pi s+{\theta }^{-1}r\right)}dx-{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\frac{q}{\theta }}d\sigma \left(x\right)\Rightarrow {\displaystyle \underset{{\Omega }_t}{\int }\left(\rho \frac{Ds}{Dt}+s\nabla .m\right)}dx\ge {\displaystyle \underset{{\Omega }_t}{\int }\left({\theta }^{-1}r\right)}dx-{\displaystyle \underset{\partial {\Omega }_t}{\int }n.\frac{q}{\theta }}d\sigma \left(x\right)$$E45

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

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

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

$$F={\nabla }_{X}x(X,t)\to J:=\det (F)$$E47

with $X,x$ the Lagrangian end Eulerian positions in the referential and actual configurations denoted $${\Omega }_R,{\Omega }_t$$ respectively, as the product of the growth deformation gradient $$F_g$$ and the growth accommodation mapping $$F_a$$

$$F=F_a.F_g$$E48

The transformation gradients $$F_a,F_g,F$$ define the mappings of the tangent spaces to the various configurations. The Jacobean of the growth mapping informs about the nature of growth:

$$J_g:=\det (F_g)$$E49

Hence $$J_g<1$$ describes growth, whereas $$J_g>1$$ 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\left(F_a;X\right)$$ with argument the reversible part of the transformation gradient (a possible explicit dependence upon the Lagrangian variable is included for heterogeneous media) as

$$P:={\partial }_{F}W\left(F_a;X\right)$$E50

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 .J_q-p_i+J_k.F_k$$;

$$\rho \dot{s}=-\nabla .J_s+{\sigma }_s$$E51

with $$J_q$$ the heat diffusion flux, $$J_s:=\frac{1}{\theta }\left(J_q-{\mu }_i{J}_i\right)$$the total entropy flux, $$J_k$$the diffusion flux of the k-specie, $$F_k\left(x,t\right)$$an external force acting on the k-specie, and $${\sigma }_s$$ the entropy production, always positive (it is dissipated). Introducing the free energy density per unit mass$$\psi :=u-\theta s$$, with $$s$$ the entropy density, we then immediately obtain the rate of variation of the free energy density

$$\rho \dot{\psi }=\rho s\dot{\theta }-\nabla .J_q+\theta \nabla .J_s+J_k.F_k-p_i-{\sigma }_s$$E52

The positivity of the entropy production $${\sigma }_s$$in previous inequality then expresses as

$$\rho \dot{\psi }\le -p_i+J_k.F_k+\nabla .\left(J_q-{\mu }_i{J}_i\right)$$E53

The principle of virtual power $$\frac{dK}{dt}=P_e+P_i$$($$K$$is the kinetic energy, $$P_e,P_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}+{\displaystyle \underset{\Omega }{\int }\rho \dot{\psi }dx}\le P_e+{\underset{¯}{J}}_k.{\underset{¯}{F}}_k+{\Phi }_m+{\Phi }_q$$E54

with $${\Phi }_q:={\displaystyle \underset{\partial \Omega }{\int }{J}_q.n}d\sigma$$ and $${\Phi }_m:=-{\displaystyle \underset{\partial \Omega }{\int }{\mu }_i{J}_i.n}d\sigma$$ respectively the flux of heat and mass through the boundary of$$\Omega$$. Previous inequality traduces the fact that the flux of mechanical work 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_a$$, the concentration of chemical specie $$n_k$$and the temperature$$\theta$$, so that Clausius-Duhem inequality (3.7) becomes in material format:

$$\begin{array}{l}\frac{D}{Dt}\left(\rho J\psi \right)\le T:\dot{F}+\nabla .\left(J_q-{\mu }_i{J}_i\right)\to \\ \rho J\psi I:{\dot{F}}_g.F_g{}^{-1}+\rho J\frac{\partial \psi }{\partial F_a}{\dot{F}}_a+\rho J\frac{\partial \psi }{\partial n_k}{\dot{n}}_k\\ +\rho J\frac{\partial \psi }{\partial \theta }\dot{\theta }\le T:\left({\dot{F}}_a.F_g+F_a.{\dot{F}}_g\right)+Div{J}_q-{\mu }_{i}Div{J}_i-J_{i}Grad{\mu }_i\end{array}$$E55

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 $${\displaystyle \underset{\partial {\Omega }_R}{\int }{J}_{\mu }.NdA}$$ and to a source term due to growth$${\displaystyle \underset{{\Omega }_R}{\int }\Gamma \rho JdX}={\displaystyle \underset{{\Omega }_R}{\int }Tr\left({\dot{F}}_g.F_g{}^{-1}\right)\rho JdX}$$, hence

$$\begin{array}{l}\frac{D}{Dt}{\displaystyle \underset{{\Omega }_t}{\int }\rho dx}={\displaystyle \underset{{\Omega }_R}{\int }Tr\left({\dot{F}}_g.F_g{}^{-1}\right)\rho JdX}={\displaystyle \underset{{\Omega }_R}{\int }\rho J{\dot{n}}_{k}dX}\\ +{\displaystyle \underset{\partial {\Omega }_R}{\int }{J}_k.NdA}\to \rho J{\dot{n}}_k+Div{J}_k=\rho J\left({\dot{F}}_g.F_g{}^{-1}:I\right)\rho J\end{array}$$E56

The last equality is nothing else than $$\Gamma =\left(\pi +\nabla .m\right)/\rho$$ expressed in material format, identifying$$\Gamma :=Tr\left(D_g\right)$$,$$M_i=J_i$$. Accordingly, (3.9) becomes

$$\begin{array}{l}0\le \left(T.F_g{}^t-\rho J\frac{\partial \psi }{\partial F_a}\right){\dot{F}}_a-\left(\frac{\partial \psi }{\partial n_k}-{\mu }_k\right)\rho J{\dot{n}}_k-J_{i}Grad{\mu }_i\\ +\left(F_a{}^t.T.F_g{}^t-\rho J\psi I-\rho J\frac{\partial \psi }{\partial n_k}\right):{\dot{F}}_g.F_g{}^{-1}+\rho Js\dot{\theta }\end{array}$$E57

Since previous equality must hold true for arbitrary variations$${\dot{F}}_a,{\dot{n}}_k$$, the following constitutive equations for the first Piola-Kirchhoff stress and the chemical potential are obtained

$$T=\rho J\frac{\partial \psi }{\partial F_a}.F_g{}^{-t};{\mu }_k=\frac{\partial \psi }{\partial n_k}$$E58

Especially, (3.12)₁ is an alternative to (3.4) using the specific free energy instead of a strain energy potential; observe that $$\psi$$ is expressed per unit mass, in contrast to$$W\left(F_a;X\right)$$ in (3.4), expressed per unit referential volume.

The residual dissipation then writes from (3.11)

$$0\le \left\{\rho Js\dot{\theta }-J_{i}Grad{\mu }_i\right\}+\rho J\left(F_a{}^t.\frac{\partial \psi }{\partial F_a}-\left(\psi +{\mu }_k\right)I\right):L_g$$E59

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

$$\rho Js\dot{\theta }-J_{i}Grad{\mu }_i\ge 0$$;

$$\left(F_a{}^t.\frac{\partial \psi }{\partial F_a}-\left(\psi +{\mu }_k\right)I\right):L_g\ge 0$$E60

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

$$L_g=f\left({\tilde{\Sigma }}_a\right)$$E61

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

$${\tilde{\Sigma }}_a:=\rho F_a{}^t.\frac{\partial \psi }{\partial F_a}-\rho \left(\psi +{\mu }_k\right)I$$E62

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 $$\alpha ,\beta$$ in the sequel), as pictured in figure 2. The transition zone between the bulk phases will be denoted by the Greek letter $$\sigma$$ 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 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.

Figure 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^{\sigma }$$ 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^{\sigma }$$ 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

$$d{F}^{\sigma }=-S^{\sigma }dT+gdA+{\mu }_{k}d{N}^{\sigma }{}_k$$E63

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 $\alpha ,\beta$ with a separating interface $\sigma$ in-between, one can write the differential of the total number of particles as

$$dN=d{N}^{\alpha }+d{N}^{\beta }+d{N}^{\sigma }$$E64

The superficial excess or molar superficial concentrations are then defined as$n^{\sigma }:=N^{\sigma }/A$, with $A$ the area of the interface. Any extensive quantity $Z$ can be decomposed as

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

with $z^{\sigma }:=Z^{\sigma }/A$ the superficial excess quantity. Regarding surface quantities, one makes a distinction between:

• The superficial energy $\gamma$ $\left(J/m^2\right)$ - a scalar - accounting for the creation of a new surface (irreversible phenomena), with a constant number of particles.

• The purely elastic variation of the surface area, expressed by a superficial stress$$\tilde{\sigma }$$, dual to an elastic surface strain$$\tilde{\epsilon }$$.

The excess total internal energy writes

$$d{U}^{\sigma }=\theta d{S}^{\sigma }+gdA+{\mu }_{k}d{N}^{\sigma }{}_k$$E66

For the whole system, using the previous decomposition

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

one has

$$dU=\theta dS-pdV+gdA+{\mu }_{k}d{N}_k$$E68

The variation of the free energy is

$$dF=-Sd\theta -pdV+gdA+{\mu }_{k}d{N}_k$$E69

Hence, $g$is defined as the partial derivative$g={\left(\frac{\partial F}{\partial A}\right)}_{T,V,N_i}$. Combining both relations

$d{F}^{\sigma }=-S^{\sigma }d\theta +gdA+{\mu }_{k}d{N}^{\sigma }{}_k$;

$$F^{\sigma }=\gamma A+{\mu }_k{N}^{\sigma }{}_k$$E70

gives

$$S^{\sigma }dT+N^{\sigma }{}_{k}d{\mu }_k-\left(g-\gamma \right)dA+Ad\gamma =0$$E71

This leads to the differential

$$d\gamma =-s^{\sigma }d\theta -n^{\sigma }{}_{k}d{\mu }_k+\left(g-\gamma \right)dA/A$$E72

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

$$g=\gamma +A{\left(\frac{\partial \gamma }{\partial A}\right)}_{T,{\mu }_i}$$E73

In the case of a purely elastic stretch, previous formula specializes to the relation$g=\gamma +A{\left(\frac{\partial \gamma }{\partial A}\right)}^{el}{}_{T,{\mu }_i}$.

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\left(P,T,n_i,A\right)$$, with respect to the formed area$$A$$, at constant temperature$$T$$, pressure$$P$$, and number of moles of each component$$n_i$$, viz

$$dG=-VdP-SdT+{\mu }_{i}d{N}_i+\gamma dA$$E74

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};\text{S}=-\frac{\partial G}{\partial T};{\mu }_i=\frac{\partial G}{\partial N_i};\gamma =\frac{\partial G}{\partial A}$$E75

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

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

with$$n_i^{\sigma }:=N_i/A$$, the surface excess of the i^th 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^{\sigma }$$, viz

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

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^{\sigma }=g^{\sigma }$$. 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 $$\sigma$$ 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$${\mu }_i$$. Imagine a modification of the temperature by$$dT$$, and of the i^th chemical potential by$$d{\mu }_i$$, at fixed surface area; $$d{n}_i$$particles from the bulk will enter the solid phase $$\sigma$$ from the bulk, and ithechange of Helmholtz free energy $$F^{\sigma }$$ will be

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

with $$S^{\sigma }$$ the entropy of the phase$$\sigma$$. Consider next a new system for which $$T$$ and $${\mu }_i$$ are returned to their initial values, but with the surface area increased by$$d{A}^{\sigma }$$, and modify thereafter the temperature by$$dT$$, and the i^th chemical potential by$$d{\mu }_i$$; the variation of free energy of this system of larger area is

$$dF{\text{'}}^{\sigma }\approx -S{\text{'}}^{\sigma }dT+{\mu }_{i}dN{\text{'}}_{i,\sigma }$$E79

Subtracting both variations of Helmholtz free energy by unit surface gives

$$\frac{dF{\text{'}}^{\sigma }-d{F}^{\sigma }}{d{A}^{\sigma }}\approx -\frac{\left(S{\text{'}}^{\sigma }-S^{\sigma }\right)}{d{A}^{\sigma }}dT+{\mu }_i\frac{d\left(N{\text{'}}_{i,\sigma }-N_{i,\sigma }\right)}{d{A}^{\sigma }}$$E80

Introducing therein the definitions of the specific surface Helmholtz energy$$f^{\sigma }:=\left(dF{\text{'}}^{\sigma }-d{F}^{\sigma }\right)/d{A}^{\sigma }$$, the specific surface entropy$$s^{\sigma }:=\left(S{\text{'}}^{\sigma }-S^{\sigma }\right)/d{A}^{\sigma }$$, and the surface excess$$n_i^{\sigma }:=d\left(N{\text{'}}_{i,\sigma }-N_{i,\sigma }\right)/d{A}^{\sigma }$$leads to

$$d{f}^{\sigma }=-s^{\sigma }dT+{\mu }_{i}d{\Gamma }_i$$E81

But one can also express the specific surface Helmholtz energy as$$f^{\sigma }=\gamma +{\mu }_i{\Gamma }_i$$, hence

$$d{f}^{\sigma }=d\gamma +{\mu }_{i}d{\Gamma }_i+{\Gamma }_{i}d{\mu }_i=-s^{\sigma }dT+{\mu }_{i}d{\Gamma }_i$$E82

and thus finally

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

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$${\tilde{\sigma }}_{ij}$$. It is related to $$\gamma$$ by

$${\tilde{\sigma }}_{ij}=\gamma {\delta }_{ij}+\frac{\partial \gamma }{\partial {\tilde{\epsilon }}_{ij}}$$E84

The second order tensor $${\tilde{\epsilon }}_{ij}$$ is the strain (a small perturbation scheme is presently adopted) induced by the component $${\tilde{\sigma }}_{ij}$$ acting in the j^th direction per unit length of the edge normal to the i^th 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:

1. 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 stretched cube along the xy plane, requiring the work$$W_2=2\left(\gamma +d\gamma \right)(1+dx)$$, with $$d\gamma$$ the variation of the specific surface work $$\gamma$$ due to the stretch $$dx$$ (factor 2 arises since two surfaces are created, and the factor $$(1+dx)$$ since the specific surface work applies per unit surface area).

2. Separate the original unit cube into two parts along an xy plane, requiring the work$$W_3=2\gamma$$, and stretch each half by $$dx$$ in the x direction, at fixed y edge, but varying z 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, hence$$W_0+W_2=W_1+W_3$$, viz

$$W_0+2\left(\gamma +d\gamma \right)(1+dx)=2\gamma +W_1$$E85

The difference $$W_1-W_0$$ is the work due to the stretching operations of both processes, and can be equalized to the x-component $${\tilde{\sigma }}_{xx}$$of a force in the newly formed surface times the distance $$2dx$$through which this force acts, hence

$$2{\tilde{\sigma }}_{xx}dx=W_1-W_0$$E86

The strain $$\Delta {\epsilon }_{xx}=dx$$ (since the other side has unit length), hence

$$2\gamma +2\gamma \Delta {\epsilon }_{xx}+2\Delta \gamma +2\Delta \gamma \Delta {\epsilon }_{xx}=2\gamma +2{\tilde{\sigma }}_{xx}\Delta {\epsilon }_{xx}$$E87

Due to the equalities

$$\Delta \gamma \Delta {\epsilon }_{xx}\approx 0$$;

$$\Delta \gamma /\Delta {\epsilon }_{xx}\approx d\gamma /d{\epsilon }_{xx}$$E88

it finally results

$${\tilde{\sigma }}_{xx}=\gamma +\frac{d\gamma }{d{\epsilon }_{xx}}$$E89

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

$${\tilde{\sigma }}_{xy}=\frac{d\gamma }{d{\epsilon }_{xy}}$$E90

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

$${\tilde{\sigma }}_{ij}=\gamma {\delta }_{ij}+\frac{d\gamma }{d{\epsilon }_{ij}}$$E91

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}=g\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$E92

Lastly, consider a square section in the xy plane of side $${\left(A^{\sigma }\right)}^{1/2}$$ and imagine an extension of the x edge by$${\epsilon }_{xx}{\left(A^{\sigma }\right)}^{1/2}$$; the required work is$$W_1={\tilde{\sigma }}_{xx}{\left(A^{\sigma }\right)}^{1/2}{\epsilon }_{xx}{\left(A^{\sigma }\right)}^{1/2}$$. Extend next the y edge by$${\epsilon }_{yy}{\left(A^{\sigma }\right)}^{1/2}$$, with an expanded work given by

$$W_2={\tilde{\sigma }}_{yy}{\left(A^{\sigma }\right)}^{1/2}\left(1+{\epsilon }_{xx}\right){\epsilon }_{yy}{\left(A^{\sigma }\right)}^{1/2}$$E93

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 }\gamma \right)=W_1+W_2={\tilde{\sigma }}_{xx}{\epsilon }_{xx}{A}^{\sigma }+{\tilde{\sigma }}_{yy}\left(1+{\epsilon }_{xx}\right){\epsilon }_{yy}{A}^{\sigma }\approx gd{A}^{\sigma }$$E94

due to the equality$$A^{\sigma }\left({\epsilon }_{xx}+{\epsilon }_{yy}\right)=d{A}^{\sigma }$$. Therefore, one has

$$gd{A}^{\sigma }=\gamma d{A}^{\sigma }+A^{\sigma }d\gamma \Rightarrow g=\gamma +\frac{d\gamma }{d{A}^{\sigma }}$$E95

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 $$\alpha ,\beta$$ are initially separated and interface-free, and are in a thought experiment imagined to be joined along a plane to generate the $$\alpha /\beta$$ interface. Since mass is conserved, any change in the thermodynamic quantities of the whole system are due to the new $$\alpha /\beta$$ interface, coined excess values of the corresponding quantities, denoted with a subscript $$\gamma$$ 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)

$$d{G}_1=\left(V_{\alpha }+V_{\beta }\right)dP-\left(S_{\alpha }+S_{\beta }\right)dT$$;

$$d{G}_2=\left(V_{\alpha }+V_{\beta }+V_{\gamma }\right)dP-\left(S_{\alpha }+S_{\beta }+S_{\gamma }\right)dT+{\gamma }^{*}dA$$E96

with $${\gamma }^{*}$$ 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_{\alpha },m_{\beta }$$ is

$$d{G}_{\gamma }=d{G}_2-d{G}_1=V_{\gamma }dP-S_{\gamma }dT+{\gamma }^{*}dA$$E97

which may be interpreted from an energetic point of view as follows: the term $$V_{\gamma }dP$$ is the mechanical work done against the external force field, the contribution $$S_{\gamma }dT$$ represents the heat of formation of the interface, and $${\gamma }^{*}dA$$ is the mechanical work done against the internal force field of both phases $$\alpha ,\beta$$ by motion of the molecules from the bulk to generate a new interface. The excess free energy of formation of the interface, potential$$G_{\gamma }$$, is the additional free energy required to form the interface from fixed masses of the pre-existing bulk phases$$\alpha ,\beta$$. The above equations implicitly use the conservation of mass, equation

$$n_{total}=n_{\alpha }+n_{\beta }$$E98

and the definition of the excess interfacial volume $$V_{\gamma }$$ from the contributions to the total volume after interface formation (balance law for the volume)

$$V_{total}=V_{\alpha }+V_{\beta }+V_{\gamma }$$E99

In contrast to this treatment, Gibbs assumes a conservation of the total volume as$$V_{total}=V_{\alpha }+V_{\beta }$$, but with addition of the new mass $$n^{\sigma }$$ such that

$$n_{total}=n_{\alpha }+n_{\beta }+n^{\sigma }$$E100

As a compensation for the volume change accompanying the formation of the interface; hence, $$n^{\sigma }$$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_{\gamma }$$ allows writing relations between partial derivatives as the analogues of the bulk phase Maxwell relations

$${\left(\frac{\partial {\gamma }^{*}}{\partial T}\right)}_{P,A,n_{\alpha ,\beta }}=-{\left(\frac{\partial S}{\partial A}\right)}_{P,T,n_{\alpha ,\beta }}=-S_{\gamma }^{*}$$;$${\left(\frac{\partial {\gamma }^{*}}{\partial P}\right)}_{T,A,n_{\alpha ,\beta }}={\left(\frac{\partial V^{*}}{\partial A}\right)}_{P,T,n_{\alpha ,\beta }}=V_{\gamma }^{}$$.

$${\left(\frac{\partial V_{\gamma }}{\partial T}\right)}_{P,A,n_{\alpha ,\beta }}=-{\left(\frac{\partial S}{\partial P}\right)}_{A,T,n_{\alpha ,\beta }}$$E101

The introduced quantities $$S_{\gamma }^{*},V_{\gamma }^{}$$ are respectively the interfacial excess entropy and the specific interfacial excess volume; the compact notation $$n_{\alpha ,\beta }$$ stands for the two quantities$$\left\{n_{\alpha ,}{n}_{\beta }\right\}$$.

The specific interfacial excess energy is obtained by simply integrating the differential $$d{G}_{\gamma }=V_{\gamma }dP-S_{\gamma }dT+{\gamma }^{*}dA\Rightarrow {\left(G_{\gamma }^{*}\right)}_{n_{\alpha ,\beta }}={\gamma }^{*}$$ at constant pressure and temperature, introducing the specific interfacial excess energy$$G_{\gamma }^{*}:=G_{\gamma }^{}/A$$. Last relation implies that the temperature and pressure dependence of $${\gamma }^{*}$$ can be determined from those of$$G_{\gamma }^{*}$$. The specific interfacial excess energy is obtained from a Legendre transform to $$d{G}_{\gamma }=V_{\gamma }dP-S_{\gamma }dT+{\gamma }^{*}dA$$ and substitution of the previous interfacial Maxwell relations, thus

$${\left(U_{\gamma }^{*}\right)}_{n_{\alpha ,\beta }}={\gamma }^{*}+T{S}_{\gamma }^{*}-P{V}_{\gamma }^{*}$$E102

It immediately results the specific interfacial excess enthalpy

$${\left(H_{\gamma }^{*}\right)}_{n_{\alpha ,\beta }}={\gamma }^{*}+T{S}_{\gamma }^{*}\equiv {\gamma }^{*}-T{\left(\frac{\partial {\gamma }^{*}}{\partial T}\right)}_{P,A,}{}_{n_{\alpha ,\beta }}$$E103

with $$H_{\gamma }^{*}$$ 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.

Figure 3.

The broken bond model for surface energy

The reversible work to form new surface area, parameter$$\gamma$$, 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 is direction dependent, and is given by the expressions $$\cos \theta /a$$and $$\sin \theta /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 $$\theta$$ 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){\epsilon }_b/\left(2a^2\right)$$E104

with $${\epsilon }_b$$ 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 $${\gamma }_i$$ over all crystal faces, a concept introduced in 1878 by J. W. Gibbs, considering constant pressure, volume, temperature and molar mass:

$$Min{\displaystyle \sum _i{A}_i{\gamma }_i}$$E105

at constant energy, hence adding the constraint$$dE=0={\displaystyle \sum _i{\gamma }_{i}d{A}_i}$$. The dependence of $$\gamma$$ 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$$\gamma \left(n\right)$$, 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).

Figure 4.

Wulff’s construction to calculate the minimizing surface for a fixed volume

with anisotropic surface tension $$\gamma \left(n\right)$$

5. Model of surface growth with application to bone remodeling

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).

5.1. Material driving forces for surface growth

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 $f_S$ (surface density) and to line densities $p_{\tau },p_{\nu }$ defined as the projections onto the unit vectors ${\tau }_g,{\nu }_g$resp. along the contour of the open growing surface $S_g$ (Figure 5); hence, those line densities are respectively tangential and normal to the surface $S_g$ (forces acting in the tangent plane).

Figure 5.

Tissue element under growth: elements of differential geometry.

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

$$\begin{array}{l}V={\displaystyle \underset{{\Omega }_g}{\int }{W}_0\left(F\right)d{x}_g}+{\displaystyle \underset{S_g}{\int }{\psi }^S\left(\tilde{F},N;X_S\right)d{\sigma }_g}+{\displaystyle \underset{S_g}{\int }{\mu }_k{n}_k^{\sigma }d{\sigma }_g}\\ -{\displaystyle \underset{S_g}{\int }{f}_S.\tilde{x}d{\sigma }_g}-{\displaystyle \underset{\partial S_g}{\int }{p}_{\tau }\tilde{\tilde{x}}.{\tau }_{g}d{l}_g}-{\displaystyle \underset{\partial S_g}{\int }{p}_{\nu }\tilde{\tilde{x}}.{\nu }_{g}d{l}_g}\end{array}$$E106

Thereby, the growing solid surface is supposed to be endowed with a volumetric density $W_0\left(F\right)$ depending upon the transformation gradient$F:={\nabla }_{X}x$, a surface energy with density ${\psi }^S\left(\tilde{F},N;X_S\right)$ per unit reference surface, depending upon the surface gradient$\tilde{F}$, the unit normal vector $N$ to$S_g$, and possibly explicitly upon the surface position vector $X_S$ on $S_g$ (no tilda notation is adopted here since the support of $X_S$ is strictly restricted to the surface$S_g$), and chemical energy${\mu }_k{n}_k^{\sigma }$, with ${\mu }_k$ the chemical potential of the surface concentration of species$n_k^{\sigma }$. The surface gradient $\tilde{F}$ maps material lengths (or material tangent vectors) onto the deformed surface; it is elaborated as the surface projection of $F$ (onto the tangent plane to${\Omega }_a$), viz

$$\tilde{F}:=F.P$$E107

The tissue element under grow is submitted to a surface force field $f_S$ (surface density) and to line densities $p_{\tau },p_{\nu }$ defined as the projections onto the unit vectors ${\tau }_g,{\nu }_g$resp. along the contour of the open growing surface $S_g$ (Figure 5); Hence, those line densities are respectively tangential and normal to the surface $S_g$ (forces acting in the tangent plane).

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

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

$$\delta \left({\displaystyle \underset{{\Omega }_g}{\int }{W}_0\left(F,X_g\right)d{x}_g}\right)={\displaystyle \underset{\partial {\Omega }_g}{\int }\left(\Sigma .\delta X_g+p.\delta x\right).Nd\left(\partial {\Omega }_g\right)}+v.t.$$E108

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

$$\Sigma :=W_{0}I-F^t.p$$

$$p:=\frac{\partial W_0}{\partial \nabla x}$$E109

As we perform material variations over an assumed fixed actual configuration, the contribution of the canonical momentum vanishes ($$\delta x=0$$). Observe that the volumetric Eshelby stress $$\Sigma$$ triggers surface growth in the sense of the boundary values taken by the normal Eshelby-like traction$$\Sigma .N$$. The variation of the surface energy contribution $${\psi }^S$$ can be expanded using the surface divergence theorem (equality (3.15) in Ganghoffer, 2010a) as

$$\delta \left({\displaystyle \underset{S_g}{\int }{\psi }^S\left(\tilde{F},N;X_S\right)d{\sigma }_g}\right)={\displaystyle \underset{S_g}{\int }\left[{\nabla }_S.\tilde{\Sigma }-\Pi .K^t.{\partial }_N{\psi }^S+{\left({\partial }_{X_S}{\psi }^s\right)}_{\exp l}+{\tilde{F}}^T.f_S\right].\delta X_{S}d{\sigma }_g}$$E110

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

$$\tilde{T}:={\partial }_{\tilde{F}}{\psi }^S\to \tilde{\Sigma }:={\tilde{F}}^T.\tilde{T}-{\psi }^s{I}_s$$E111

basing on the surface stress$$\tilde{T}$$. The Lagrangian curvature tensor is defined as$K:=-{\nabla }_{R}N$. The chemical potential as the partial derivative of the surface energy density with respect to the superficial concentration

$${\mu }_k:={\frac{\partial {\psi }^S}{\partial n_k^{\sigma }}}_{|X,F,N}\equiv {\mu }_k\left(n_k^{\sigma }\right)$$E112

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

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

$${\tilde{\Upsilon }}_g:=\Sigma .N+{\nabla }_S.\tilde{\Sigma }-P.K^t.{\partial }_N{\psi }^S+{\mu }_k{\nabla }_S{n}_k^{\sigma }-f_S$$E113

itself built from the surface stress$$\tilde{T}:={\partial }_{\tilde{F}}{\psi }^S$$, and on the curvature tensor $K:=-{\nabla }_{R}N$in the referential configuration.

5.2. Bone remodeling

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

$$\frac{d}{dt}{\displaystyle \underset{{\Omega }_g}{\int }{\rho }_{g}d{x}_g}={\displaystyle \underset{S_g}{\int }{\rho }_g{V}_S.Nd{\sigma }_g}$$E114

The quantity ${\rho }_g{V}_S.Nd{\sigma }_g$therein represents the molar flux of bone material being dissolved, hence

$${\rho }_g{V}_{N}d{\sigma }_g=MJd{\sigma }_g$$E115

with $V_N$ the normal surface velocity, $M$the bone mineral molar mass, and $J\equiv {\rho }_g{V}_N/M$ the molar influx of minerals (positive in case of bone apposition, and negative when resorption occurs). Clearly, the previous expression shows that the knowledge of the normal surface growth velocity determines the molar influx of minerals. Estimates of the order of magnitude of the dissolution rate given in (Christoffersen at al., 1997), for a pH of 7.2 (although much higher compared to the pH for which bone resorption takes place) and at a temperature of 310K, are indicative of values of the molar influx in the interval$J\in \left[{10}^{-9}{\mathrm{,1.8.10}}^{-8}\right]mol.s^{-1}.m^{-2}$. The osteoclasts responsible for bone resorption attach to the bone surface, remove the collageneous fraction of the material by transport phenomena, and diffuse within the material. This osteoclasts activity occurs at a typical scale of about$50\mu m$, which is much larger compared to the characteristic size of the ultrastructure; the resorption phase takes typically 21 days (the complete remodeling cycle lasts 3 months). The osteoclasts, generate an acid environment causing simultaneously the dissolution of the mineral - hydroxyapatite, a strong basic mineral${\left[C{a}_3{\left(P{O}_4\right)}_2\right]}_{3}Ca{\left(OH\right)}_2$, abbreviated HA in the sequel - and the degradation of the collageneous fraction of the material. The metabolic processes behind bone remodeling are very complicated, with kinetics of various chemical substances, see (Petrtyl and Danesova, 1999).

The pure chemical driving force represents the difference of the chemical potential externally supplied ${\mu }_e$ (biochemical activity generated by the osteoclasts) with the chemical potential of the mineral of the solid phase, denoted${\mu }_{\min }$; it can be estimated from the change of activity of the $H^{+}$ cation (Silva and Ulm, 2002):

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

This chemical driving force is the affinity conjugated to the superficial concentration of minerals, denoted $$n^{\sigma }(t)$$ in the sequel. The conversion to mechanical units $\Delta \mu$ is done, considering a density of HA $\rho =3000kg/m^3$ (5.1), hence$\left(\rho /M\right)\Delta \mu =-20MPa$, according to (Silva and Ulm, 2002); the negative value means that the dissolution of HA is chemically more favorable (bone resorption occurs).

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

According to experiments performed by (Currey, 1988), the elastic modulus is assumed to scale uniformly versus the bone density according to

$$E=E_{\max }{\rho }_S{\left(t\right)}^p$$E117

with $${\rho }_S\left(t\right)$$ the surface density of mineral, $E_{\max }=15GPa$(Reilly and Burstein, 1975) the maximum value of the tensile modulus, and $p$ a constant exponent, here taken equal to 3 (Currey, 1988; Ruinerman et al., 2005).

Following the representation theorems for isotropic scalar valued functions of tensorial arguments, the surface strain energy density ${\psi }_{mech}^S\left(\tilde{F},N;X_S\right)$of mechanical origin is selected as a function of the curvature tensor invariants, viz the mean and Gaussian curvatures, the invariants of the surface Cauchy-Green tensor $\tilde{C}:={\tilde{F}}^t.\tilde{F}$ and of its square. The following simple form depending on the second invariant of the linearized part of $\tilde{C}\cong I+2\tilde{\epsilon }$ is selected, adopting the small strain framework, viz, hence

$${\psi }_{mech}^S\left(\tilde{\epsilon }\right)=\frac{A}{2}Tr{(\tilde{\epsilon })}^2+B\left(\tilde{\epsilon }:\tilde{\epsilon }\right)$$E118

Figure 6.

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

with $\tilde{\epsilon }=P.\epsilon =I_S.\epsilon -\epsilon .\left(e_r\otimes e_r\right)$ the surface strain (induced by the existing volumetric strain), and $A,B$ mechanical properties of the surface, expressing versus the surface density of minerals and the maximum value of the traction modulus as (the Poisson ratio is selected as$\nu =0.3$)

$$A=\frac{E_{\max }{\rho }_S{\left(t\right)}^3\nu }{\left(1-2\nu \right)\left(1+\nu \right)};B=\frac{E_{\max }{\rho }_S{\left(t\right)}^3}{2\left(1+\nu \right)}$$E119

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

$$\tilde{T}\equiv \tilde{\sigma }:=\frac{\partial {\psi }_{mech}^S}{\partial \tilde{\epsilon }}=A\tilde{\epsilon }+2Btr\left(\tilde{\epsilon }\right)I_S$$E120

The unknowns of the remodeling problem are the normal velocity of the bone surface$$V_N\left(t\right)$$, the surface density of minerals $${\rho }_S\left(t\right)$$ and its superficial concentration. We shall herewith simulate the resorption of a hollow bone submitted to a composite applied stress, consisting of the superposition of an axial and a radial component, as

$$\sigma ={\sigma }_{rr}{e}_r\otimes e_r+{\sigma }_{zz}{e}_z\otimes e_z$$E121

in the cylindrical basis$\left(e_r,e_{\theta },e_z\right)$; this applied stress generates a preexisting homogeneous stress state within the bulk material, inducing a surface stress given by

$$\tilde{\sigma }=P.\sigma ={\sigma }_{zz}{e}_z\otimes e_z$$E122

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

$${\tilde{\Upsilon }}_g{}_{gN}=\frac{1}{r_i\left(t\right)}\left\{\frac{1}{8}\Delta \mu n_{}^{\sigma }\left(t\right)+\frac{A+2B}{8\left(A+B\right)B}{\sigma }_{zz}{}^2\right\}$$E123

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

$$V_N\left(t\right)=C{\tilde{\Upsilon }}_g{}_{gN}\left(t\right)$$E124

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

$$J={\rho }_g{V}_N/M={10}^{-8}mol.s^{-1}.m^{-2}\Rightarrow V_N=JM/{\rho }_g\approx {\mathrm{3.3.10}}^{-12}m/s\approx 0.286\mu m/day$$E125

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

The mass balance equation for the surface density of minerals $${\rho }_S$$writes

$${\dot{\rho }}_S+{\rho }_S{\nabla }_S.\tilde{V}={\Gamma }^S{\rho }_S$$E126

expressing as the following conservation law

$$\frac{{\dot{\rho }}_S}{{\rho }_S}-\frac{V_N}{r_i(t)}={\Gamma }_0^S\iff {\rho }_S\left(t\right)={\rho }_s^0\frac{r_0}{r_i\left(t\right)}\exp \left({\Gamma }_0^{S}t\right)$$E127

The initial surface density of minerals$${\rho }_s^0={\rho }_S\left(t_0\right)$$, is evaluated from the bulk density of HA, viz$3000kg/m^3$, and the estimated thickness of the attachment region of osteoclasts, about $7\mu m$ (Blair, 1998), hence$${\rho }_s^0\approx {\mathrm{2.1.10}}^{-2}kg/m^2$$.

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

$${\left[C{a}_3{\left(P{O}_4\right)}_2\right]}_{3}Ca{\left(OH\right)}_2+8H^{+}\to 10C{a}^{2+}+6HP{O}_4^{2-}+2H_{2}O$$E128

The kinetic equation is chosen as:

$$\frac{\partial n^{\sigma }\left(t\right)}{\partial t}=-\tilde{\gamma }{\rho }_s\left(t\right)n^{\sigma }\left(t\right)=-\tilde{\gamma }{r}_i\left(t\right)\frac{{\rho }_s^0}{r_0}\exp \left({\Gamma }_0^{S}t\right)n^{\sigma }\left(t\right)$$E129

incorporating the density of minerals. The rate coefficient of dissolution of HA, namely the parameter$$\tilde{\gamma }$$, is taken at room temperature from literature values available for CHA (carbonated HA, similar to bone), viz $$\tilde{\gamma }\approx {\mathrm{2.2.10}}^{-4}{s}^{-1}$$ (Hankermeyer et al., 2002).

5.3 Simulation results

The present model involves a dependency of the triplet of variables $$\left\{r_i\left(t\right),{\rho }_S\left(t\right),n^{\sigma }\left(t\right)\right\}$$ solution of the set of equations (5.15) through (5.19) on a set of parameters, arising from initial conditions satisfied by those variables:

• The initial concentration of minerals $$n_0^{\sigma }$$ is taken as unity, viz$$n_0^{\sigma }=1\text{mol}{\text{.m}}^{-3}$$.

• The initial radius $$r_0:=r_i\left(0\right)$$ is estimated as $$r_0=1.6cm$$ for the diaphysis of the human femur (Huiskes and Sloof, 1981). The evolution versus time of the internal radius obtained by time integration of the normal velocity expressed in (5.16).

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

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

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

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

$${\sigma }_{zz}^{crit}(t)\approx {\mathrm{9.4.10}}^{10}{\rho }_S{\left(t\right)}^{3/2}{n}_{}^{\sigma }{\left(t\right)}^{1/2}$$E130

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

Figure 7.

Evolution vs. time of the surface growth velocity for three stress levels: ${\sigma }_{zz}=1MPa$(thick line), ${\sigma }_{zz}=2MPa$(dashed line), ${\sigma }_{zz}=5MPa$(dash-dotted line).

Figure 8.

Evolution of the superficial density of HA versus time for three stress levels. ${\sigma }_{zz}=1MPa$(thick line), ${\sigma }_{zz}=2MPa$(dashed line), ${\sigma }_{zz}=5MPa$(dash-dotted line).

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

Figure 9.

Evolution of the internal radius of the diaphysis of the human femur (in microns) versus time. Applied stress above the critical stress level:$${\sigma }_{zz}=0.02MPa$$.

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

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

6. Concluding remarks

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

A macroscopic model of bone external remodeling has been developed, basing on the thermodynamics of surfaces and with the identified configurational driving forces promoting surface evolution. The interactions between the surface diffusion of minerals and the mechanical driving factors have been quantified, resulting in a relatively rich model in terms of physical and mechanical parameters. Applications of the developed formalism to real geometries

Works accounting for the multiscale aspect of bone remodeling have emerged in the literature since the late nineteen’s considering cell-scale (a few microns) up to bone-scale (a few centimeters) remodeling, showing adaptation of the 3D trabeculae architecture in response to mechanical stimulation, see the recent contributions (Tsubota et al., 2009; Coelho et al., 2009) and the references therein. It is likely that one has in the future to combine models at both micro and macro scales in a hierarchical approach to get deeper insight into the mechanisms of Wolff’s law.

The present modeling framework shall serve as a convenient platform for the simulation of bone remodeling with the consideration of real geometries extracted from CT scans. The predictive aspect of those simulations is interesting in a medical context, since it will help doctors in adapting the medical treatment according to short and long term predictions of the simulations.