Open access peer-reviewed chapter

Thermodynamics of Surface Growth with Application to Bone Remodeling

By Jean-François Ganghoffer

Submitted: November 5th 2010Reviewed: March 19th 2011Published: November 2nd 2011

DOI: 10.5772/19978

Downloaded: 1596

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(A.B)ij=AikBkj. 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 Awith volumetric density a=a(x,t)satisfies a prototype balance law of the form

a(x,t)t=.Ja(x,t)+σa(x,t)E1

with Ja(x,t)the flux density of a(x,t)and σa(x,t)the local production (or destruction) ofa(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ρk, obey the local balance law (Vidal et al., 1994)

ρkt=.(ρku+Jk)+Mkα=1..rναkJαE2

with u:=1ρk=1nρkukthe local barycentric velocity, Mkthe molar mass, and ναkthe stoechiometric coefficients in the reactionα, such that the variation of the mass dmkof the species k due to chemical reactions expresses as

dmk=Mkα=1..rναkξα,k=1..nE3

wherein ξαdenotes the degree of advancement of reactionα. The molar masses Mksatisfy the global conservation law (due to Lavoisier)

k=1nναkMk=0,α=1..rE4

Observe that the total flux of mass is the sum of a convective flux ρkuand a diffusive fluxJk; the mass production is identified as the contributionMkα=1..rναkJα. In this viewpoint, the system is in fact closed, since the balance law satisfied by the global density ρ=k=1nρkwrites (Vidal et al., 1994) accounting for the relationj=1nJk=j=1nρkuk=0, as

ρt=.(ρu)+k=1nα=1..rMkναkJα.(ρu)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, nk=mk/Mk, with mkthe mass of the same constituent. The molar concentration is defined asck=nk/V, its inverse being called the partial molar volume. The partial mole number nksatisfies the balance equation

nkt=divJk+inktE6

with Jkthe flux of species k and inktits production term, given by De Donder definition of the rate of progress of the jth chemical reaction

inkt=j=1rνkjξ˙jE7

The two previous equalities enter into Gibbs relation as

ρu˙=θρs˙e+θρs˙i+σ:ε˙kμkρMdivJk+kμkρMjνkjξ˙jE8

with θthe temperature and μkthe chemical potential of constituent k. The chemical affinity in the sense of De Donder is defined as the force conjugated to the rate ξ˙j

Aj=k(μkρM)νkj=k(μk/V)νkjE9

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

ρs˙=1θρu˙1θσ:ε˙+k1θ(μk1V)divJk+j1θAjξ˙jE10

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

ρu˙=.Jq+ρw˙E11

with Jqthe heat flux, and the term ρwis relative to all forms of work. One shall isolate the flux-like contributions in the entropy variation, which after a few transformations writes

ρs˙=ρs˙e+ρs˙i=1θ.Jq+k.(Jkμkθ1V)kJk.(μkθ1V)+1θ(ρw˙σ:ε˙)+j1θAjξ˙jE12

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

1θσ:ε˙=1θσ:u=.(1θσ:u)u..(1θσ)E13

Hence, the rate of the entropy density decomposes into

ρs˙=ρs˙e+ρs˙i=.(Jq+kJkμkθ1V)+Jq.(1θ)kJk.(μkθ1V)+1θ(ρw˙σ:ε˙)+j1θAjξ˙jE14

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

Js=Jq+kJkμkθ1VE15

and of the internal entropy production

ρs˙i=Jq.(1θ)kJk.(μkθ1V)1θ(σ:ε˙ρw˙)+j1θAjξ˙jE16

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

ρu˙=.Jq+σ:ε˙+k(μk/V)n˙kE17

One has assumed in this alternative that the mechanical power ρw˙=σeq:udoes not include a flux contribution, hence only the heat diffusion contributes to the flux of internal energy. The contribution σ:ε˙+k(μk/V)n˙kis identified to the termρw˙. Previous equality combined with the second principle, equality ρs˙=.(Jqθ)+ρs˙i(the entropy flux resumes to the sole heat flux), delivers after a few manipulations the variation of the internal energy as

ρu˙=(θρs˙JqθθTρs˙i)+σ:ε˙+k(μk/V)n˙kE18

Hence, the internal entropy production is identified as

θρs˙i=ρ(u˙Ts˙)Jqθθ+σ:ε˙+k(μk/V)n˙kE19

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

θρs˙i=ρ(ψ˙+sθ˙)Jq.θθ+σ:ε˙+k(μk/V)n˙kE20

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Ω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 asF(a)=av; the overall convective flux of a across the closed surface Ωtexpresses then as

ϕconv(a)=Ωta(vw).ndAE21

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

ϕdiff(a)=ΩtJ(a).ndAE22

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Ωt.

The convective derivative along the vector field wof the field a=a(x,t)writes

δwaδt=(at)x+a.wE23

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

dadtδvaδt=δwaδt+a.(vw)E24

The derivative of the volume integral A:=Ωtadxis next calculated, according to Leibniz rule:

DDtΩtadx=Ωtatdx+Ωta(w.n)dAE25

with wthe velocity field of the points onΩ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

DDtΩtadx=ΩtΠdxΩt{a(vw)+J(a)}.ndAE26

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

Ωtatdx+Ωta(w.n)dA=ΩtΠdxΩt{a(vw)+J(a)}.ndAE27

The corresponding local balance law is obtained after elimination of the velocityw, hence

δvaδt+adivv=ΠdivJ(a)at+div(av)=ΠdivJ(a)E28

Mass balance: the mass balance equation is deduced from the identificationa=ρ, the actual density. Hence, (23) gives

DDtΩtρdx=ΩtπdxΩtJ(ρ).ndAE29

The strong form of the balance law of mass writes finally

δvρδt=ΠdivJ(ρ)ρ.vE30

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

DDtΩtρdx=Ωt(DρDt+ρ.v)dx=Ωtπdx+ΩtmdsΩtΓρdxE31

with ρ(x,t)the actual density, πthe physical source of mass, and m:=m.nthe 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

DρDt=π+.mρ.vE32

with v(x,t):=(xt)Xthe 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(x,t), obtained using the mass balance (28)

DDtΩtρadx=Ωt{ρDaDt+a(π+.m)}dxE33

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

ρ˙+ρdiv(v)=Φρ+σρE34

with Φρ.mthe 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 Ωtasm(Ωt)=Ωtρ(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

(dmdt)source:=Ωtπdx=Ωtdxπ=ΓρE35

The time variation of chemical concentration of nutrients is due to exchange through the boundary accounted for by a fluxΩtjk.nds, and to a source term due to growthΩtΓρdx=ΩtTr(F˙g.Fg1)ρdx, hence (see (28))

DDtΩtρdx=ΩtTr(F˙g.Fg1)ρdx=Ωtρn˙kdx+Ωtjk.Ndsρn˙k+divjk=ρ(F˙g.Fg1:I)=ρΓE36

The last equality is nothing else than Γ=(π+.m)/ρ- a consequence of (28) - expressed in material format, with the identificationsΓ:=Tr(Dg),mi=ji;π=ρn˙k. The global form DDtΩRρdx=DDtΩRρJdXalso fits within the general balance law for an open system, relation (30) witha1.

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

DDtΩtρvdx=Ωtfdx+Ωtn.σdσt+Ωtπvdx+Ωtn.(mv)dσtE37

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

ρDvDt=f+divσ+(m.)vE38

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 uthe specific internal energy density, one may write the energy balance in the actual configuration as

DDtΩtρ(u+12v2)dx=Ωt{(f.v+r)+π(u+12v2)}dx+Ωtn.{σ.v+m(u+12v2)q}dσ(x)E39

with rthe volumetric heat supply (generated by growth), and qthe heat flux acrossΩt. This writing of the energy balance can be simplified using the balance of kinetic energy with volumetric densityk, obtained by multiplying (33) by the velocity and integrating overΩt, hence

k:=12ρDv2Dt=f.v+v.divσ+v.(m.)vf.v+v.divσ+m.v22Ωt12ρDv2Dtdx=Ωt(f.v+v.divσ+m.(v22))dxE40

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) witha=12v2, hence (35)

DKDt=Ωt{ρD(v22)Dt+v22(π+.m)}dx=Ωt(f.v+v.divσ+.(mv22)+πv22)dx=Ωt(f.vσ:v+πv22)dx+Ωtn.{σ.v+mv22}dσE41

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

DDt(U+K)=ΩtρDDt(u+12v2)dx+Ωt(u+12v2)(π+.m)dx=Ωt{(f.v+r)+π(u+12v2)}dx+Ωtn.{σ.v+m(u+12v2)q}dσ(x)E42

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

DUDt=Ωt(σ:v+r+πu)dx+Ωtn.(muq)dσ(x)E43

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

ρDuDt=σ:v+r+m.u.qE44

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

DDtΩtρsdxΩt(πs+θ1r)dxΩtn.qθdσ(x)Ωt(ρDsDt+s.m)dxΩt(θ1r)dxΩtn.qθdσ(x)E45

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

ρDsDtθ1rdiv(qθ)s.mE46

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=Xx(X,t)J:=det(F)E47

with X,xthe Lagrangian end Eulerian positions in the referential and actual configurations denoted ΩR,Ωtrespectively, as the product of the growth deformation gradient Fgand the growth accommodation mapping Fa

F=Fa.FgE48

The transformation gradients Fa,Fg,Fdefine the mappings of the tangent spaces to the various configurations. The Jacobean of the growth mapping informs about the nature of growth:

Jg:=det(Fg)E49

Hence Jg<1describes growth, whereas Jg>1represents resorption. Growth essentially occurs between the referential and the actual configurations.

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

P:=FW(Fa;X)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 uand sthe density of internal energy and entropy per unit mass respectively; the first and second principles of thermodynamics write (Munster, 1970)

ρu˙=.Jqpi+Jk.Fk;
ρs˙=.Js+σsE51

with Jqthe heat diffusion flux, Js:=1θ(JqμiJi)the total entropy flux, Jkthe diffusion flux of the k-specie, Fk(x,t)an external force acting on the k-specie, and σsthe entropy production, always positive (it is dissipated). Introducing the free energy density per unit massψ:=uθs, with sthe entropy density, we then immediately obtain the rate of variation of the free energy density

ρψ˙=ρsθ˙.Jq+θ.Js+Jk.FkpiσsE52

The positivity of the entropy production σsin previous inequality then expresses as

ρψ˙pi+Jk.Fk+.(JqμiJi)E53

The principle of virtual power dKdt=Pe+Pi(Kis the kinetic energy, Pe,Pibeing the virtual power of external and internal forces respectively), leads to the global form of previous inequality in Eulerian format:

dKdt+Ωρψ˙dxPe+J¯k.F¯k+Φm+ΦqE54

with Φq:=ΩJq.ndσand Φm:=ΩμiJi.ndσrespectively the flux of heat and mass through the boundary ofΩ. 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 gradientFa, the concentration of chemical specie nkand the temperatureθ, so that Clausius-Duhem inequality (3.7) becomes in material format:

DDt(ρJψ)T:F˙+.(JqμiJi)ρJψI:F˙g.Fg1+ρJψFaF˙a+ρJψnkn˙k+ρJψθθ˙T:(F˙a.Fg+Fa.F˙g)+DivJqμiDivJiJiGradμiE55

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 ΩRJμ.NdAand to a source term due to growthΩRΓρJdX=ΩRTr(F˙g.Fg1)ρJdX, hence

DDtΩtρdx=ΩRTr(F˙g.Fg1)ρJdX=ΩRρJn˙kdX+ΩRJk.NdAρJn˙k+DivJk=ρJ(F˙g.Fg1:I)ρJE56

The last equality is nothing else than Γ=(π+.m)/ρexpressed in material format, identifyingΓ:=Tr(Dg),Mi=Ji. Accordingly, (3.9) becomes

0(T.FgtρJψFa)F˙a(ψnkμk)ρJn˙kJiGradμi+(Fat.T.FgtρJψIρJψnk):F˙g.Fg1+ρJsθ˙E57

Since previous equality must hold true for arbitrary variationsF˙a,n˙k, the following constitutive equations for the first Piola-Kirchhoff stress and the chemical potential are obtained

T=ρJψFa.Fgt;μk=ψnkE58

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 toW(Fa;X)in (3.4), expressed per unit referential volume.

The residual dissipation then writes from (3.11)

0{ρJsθ˙JiGradμi}+ρJ(Fat.ψFa(ψ+μk)I):LgE59

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

ρJsθ˙JiGradμi0;
(Fat.ψFa(ψ+μk)I):Lg0E60

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

Lg=f(Σ˜a)E61

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

Σ˜a:=ρFat.ψFaρ(ψ+μk)IE62

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 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σ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σ=SσdT+gdA+μkdNσkE63

The quantity gaccounts 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

dN=dNα+dNβ+dNσE64

The superficial excess or molar superficial concentrations are then defined asnσ:=Nσ/A, with Athe area of the interface. Any extensive quantity Zcan be decomposed as

Z=Zα+Zβ+Zσ=zαVα+zβVβ+zσAE65

with zσ:=Zσ/Athe superficial excess quantity. Regarding surface quantities, one makes a distinction between:

  • The superficial energy γ(J/m2)- 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σ˜, dual to an elastic surface strainε˜.

The excess total internal energy writes

dUσ=θdSσ+gdA+μkdNσkE66

For the whole system, using the previous decomposition

Z=Zα+Zβ+ZσE67

one has

dU=θdSpdV+gdA+μkdNkE68

The variation of the free energy is

dF=SdθpdV+gdA+μkdNkE69

Hence, gis defined as the partial derivativeg=(FA)T,V,Ni. Combining both relations

dFσ=Sσdθ+gdA+μkdNσk;
Fσ=γA+μkNσkE70

gives

SσdT+Nσkdμk(gγ)dA+Adγ=0E71

This leads to the differential

dγ=sσdθnσkdμk+(gγ)dA/AE72

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=γ+A(γA)T,μiE73

In the case of a purely elastic stretch, previous formula specializes to the relationg=γ+A(γA)elT,μ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), quantityG(P,T,ni,A), with respect to the formed areaA, at constant temperatureT, pressureP, and number of moles of each componentni, viz

dG=VdPSdT+μidNi+γdAE74

whereby the multiplicative factors of the differential elements on the right-hand side of dGare the partial derivatives

V=GP;S=GT;μi=GNi;γ=GAE75

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, denotedgσ, and the change per unit area of surface created from the surrounding bulk phases, evaluated as the sum μiniσover all components,

γ=gσμiniσE76
withniσ:=Ni/A, the surface excess of the ith species. In terms of the Helmholtz energy of the whole systemF, one has the similar relation involving the Helmholtz free energy per unit areafσ, viz
γ=fσ+PeμiniσE77

with ethe thickness of the surface; in most cases, the parameter eis small, and one may neglect the contributionPe, hence one has the identificationfσ=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 Tand as a chemical reservoir for each component; accordingly, the components concentrations can be adjusted to maintain the chemical potentials at fixed specified valuesμi. Imagine a modification of the temperature bydT, and of the ith chemical potential bydμi, at fixed surface area; dniparticles from the bulk will enter the solid phase σfrom the bulk, and ithechange of Helmholtz free energy Fσwill be

dFσ=SσdT+(μi+dμi)dNi,σSσdT+μidNi,σE78

with Sσthe entropy of the phaseσ. Consider next a new system for which Tand μiare returned to their initial values, but with the surface area increased bydAσ, and modify thereafter the temperature bydT, and the ith chemical potential bydμi; the variation of free energy of this system of larger area is

dF'σS'σdT+μidN'i,σE79

Subtracting both variations of Helmholtz free energy by unit surface gives

dF'σdFσdAσ(S'σSσ)dAσdT+μid(N'i,σNi,σ)dAσE80

Introducing therein the definitions of the specific surface Helmholtz energyfσ:=(dF'σdFσ)/dAσ, the specific surface entropysσ:=(S'σSσ)/dAσ, and the surface excessniσ:=d(N'i,σNi,σ)/dAσleads to

dfσ=sσdT+μidΓiE81

But one can also express the specific surface Helmholtz energy asfσ=γ+μiΓi, hence

dfσ=dγ+μidΓi+Γidμi=sσdT+μidΓiE82

and thus finally

dγ=sσdTΓidμiE83

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

σ˜ij=γδij+γε˜ijE84

The second order tensor ε˜ijis the strain (a small perturbation scheme is presently adopted) induced by the component σ˜ijacting in the jth direction per unit length of the edge normal 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 ghalf 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 axesx,y,z, and perform two distinct operations on it:

  1. Stretch the cube reversibly along the x axis by an amountdx, 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 W0the work expanded in this transformation. Let next separate the stretched cube along the xy plane, requiring the workW2=2(γ+dγ)(1+dx), with dγthe variation of the specific surface work γ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 workW3=2γ, and stretch each half by dxin the x direction, at fixed y edge, but varying z edge. Let W1be 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, henceW0+W2=W1+W3, viz

W0+2(γ+dγ)(1+dx)=2γ+W1E85

The difference W1W0is the work due to the stretching operations of both processes, and can be equalized to the x-component σ˜xxof a force in the newly formed surface times the distance 2dxthrough which this force acts, hence

2σ˜xxdx=W1W0E86

The strain Δεxx=dx(since the other side has unit length), hence

2γ+2γΔεxx+2Δγ+2ΔγΔεxx=2γ+2σ˜xxΔεxxE87

Due to the equalities

ΔγΔεxx0;
Δγ/Δεxxdγ/dεxxE88

it finally results

σ˜xx=γ+dγdεxxE89

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

σ˜xy=dγdεxyE90

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

σ˜ij=γδij+dγdεijE91

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

σ˜ij=g(1001)E92

Lastly, consider a square section in the xy plane of side (Aσ)1/2and imagine an extension of the x edge byεxx(Aσ)1/2; the required work isW1=σ˜xx(Aσ)1/2εxx(Aσ)1/2. Extend next the y edge byεyy(Aσ)1/2, with an expanded work given by

W2=σ˜yy(Aσ)1/2(1+εxx)εyy(Aσ)1/2E93

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(Aσγ)=W1+W2=σ˜xxεxxAσ+σ˜yy(1+εxx)εyyAσgdAσE94

due to the equalityAσ(εxx+εyy)=dAσ. Therefore, one has

gdAσ=γdAσ+Aσdγg=γ+dγdAσ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 α,β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)

dG1=(Vα+Vβ)dP(Sα+Sβ)dT;
dG2=(Vα+Vβ+Vγ)dP(Sα+Sβ+Sγ)dT+γ*dAE96

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

dGγ=dG2dG1=VγdPSγdT+γ*dAE97

which may be interpreted from an energetic point of view as follows: the term VγdPis the mechanical work done against the external force field, the contribution SγdTrepresents the heat of formation of the interface, and γ*dAis 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, potentialGγ, is the additional free energy required to form the interface from fixed masses of the pre-existing bulk phasesα,β. The above equations implicitly use the conservation of mass, equation

ntotal=nα+nβE98

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

Vtotal=Vα+Vβ+VγE99

In contrast to this treatment, Gibbs assumes a conservation of the total volume asVtotal=Vα+Vβ, but with addition of the new mass nσsuch that

ntotal=nα+nβ+nσE100

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

(γ*T)P,A,nα,β=(SA)P,T,nα,β=Sγ*;(γ*P)T,A,nα,β=(V*A)P,T,nα,β=Vγ.
(VγT)P,A,nα,β=(SP)A,T,nα,βE101

The introduced quantities Sγ*,Vγare respectively the interfacial excess entropy and the specific interfacial excess volume; the compact notation nα,βstands for the two quantities{nα,nβ}.

The specific interfacial excess energy is obtained by simply integrating the differential dGγ=VγdPSγdT+γ*dA(Gγ*)nα,β=γ*at constant pressure and temperature, introducing the specific interfacial excess energyGγ*:=Gγ/A. Last relation implies that the temperature and pressure dependence of γ*can be determined from those ofGγ*. The specific interfacial excess energy is obtained from a Legendre transform to dGγ=VγdPSγdT+γ*dAand substitution of the previous interfacial Maxwell relations, thus

(Uγ*)nα,β=γ*+TSγ*PVγ*E102

It immediately results the specific interfacial excess enthalpy

(Hγ*)nα,β=γ*+TSγ*γ*T(γ*T)P,A,nα,βE103

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.

Figure 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 is direction dependent, and is given by the expressions cosθ/aand sinθ/ain the x and y directions respectively (figure 3), with athe 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

Esurf=(cosθ+sinθ)εb/(2a2)E104

with εbthe 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 γiover all crystal faces, a concept introduced in 1878 by J. W. Gibbs, considering constant pressure, volume, temperature and molar mass:

MiniAiγiE105

at constant energy, hence adding the constraintdE=0=iγidAi. 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 nthe 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 γ(n)

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 fS(surface density) and to line densities pτ,pνdefined as the projections onto the unit vectors τg,νgresp. along the contour of the open growing surface Sg(Figure 5); hence, those line densities are respectively tangential and normal to the surface Sg(forces acting in the tangent plane).

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

V=ΩgW0(F)dxg+SgψS(F˜,N;XS)dσg+SgμknkσdσgSgfS.x˜dσgSgpτx˜˜.τgdlgSgpνx˜˜.νgdlgE106

Thereby, the growing solid surface is supposed to be endowed with a volumetric density W0(F)depending upon the transformation gradientF:=Xx, a surface energy with density ψS(F˜,N;XS)per unit reference surface, depending upon the surface gradientF˜, the unit normal vector NtoSg, and possibly explicitly upon the surface position vector XSon Sg(no tilda notation is adopted here since the support of XSis strictly restricted to the surfaceSg), and chemical energyμknkσ, with μkthe chemical potential of the surface concentration of speciesnkσ. The surface gradient 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Ωa), viz

F˜:=F.PE107

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

The variation of the previously built potential energy of the growing tissue elementVis 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 surfaceSg. We refer to the recent work in (Ganghoffer, 2010a) giving the detailed calculation of the material forces for surface growth, very similar to present developments.

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

δ(ΩgW0(F,Xg)dxg)=Ωg(Σ.δXg+p.δx).Nd(Ωg)+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

Σ:=W0IFt.p
p:=W0xE109

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

δ(SgψS(F˜,N;XS)dσg)=Sg[S.Σ˜Π.Kt.NψS+(XSψs)expl+F˜T.fS].δXSdσgE110

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

T˜:=F˜ψSΣ˜:=F˜T.T˜ψsIsE111

basing on the surface stressT˜. The Lagrangian curvature tensor is defined asK:=RN. The chemical potential as the partial derivative of the surface energy density with respect to the superficial concentration

μk:=ψSnkσ|X,F,Nμk(nkσ)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 Vas the vector acting on the variation of the surface position

ϒ˜g:=Σ.N+S.Σ˜P.Kt.NψS+μkSnkσfSE113

itself built from the surface stressT˜:=F˜ψS, and on the curvature tensor K:=RNin 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

ddtΩgρgdxg=SgρgVS.NdσgE114

The quantity ρgVS.Ndσgtherein represents the molar flux of bone material being dissolved, hence

ρgVNdσg=MJdσgE115

with VNthe normal surface velocity, Mthe bone mineral molar mass, and JρgVN/Mthe 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 intervalJ[109,1.8.108]mol.s1.m2. 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 about50μ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[Ca3(PO4)2]3Ca(OH)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 μe(biochemical activity generated by the osteoclasts) with the chemical potential of the mineral of the solid phase, denotedμmin; it can be estimated from the change of activity of the H+cation (Silva and Ulm, 2002):

Δμ:=μeμmin=Rθln[H+]eq2[H+]ex2E116

This chemical driving force is the affinity conjugated to the superficial concentration of minerals, denoted nσ(t)in the sequel. The conversion to mechanical units Δμis done, considering a density of HA ρ=3000kg/m3(5.1), hence(ρ/M)Δμ=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=EmaxρS(t)pE117

with ρS(t)the surface density of mineral, Emax=15GPa(Reilly and Burstein, 1975) the maximum value of the tensile modulus, and pa 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 ψmechS(F˜,N;XS)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 C˜:=F˜t.F˜and of its square. The following simple form depending on the second invariant of the linearized part of C˜I+2ε˜is selected, adopting the small strain framework, viz, hence

ψmechS(ε˜)=A2Tr(ε˜)2+B(ε˜:ε˜)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 ε˜=P.ε=IS.εε.(erer)the surface strain (induced by the existing volumetric strain), and A,Bmechanical properties of the surface, expressing versus the surface density of minerals and the maximum value of the traction modulus as (the Poisson ratio is selected asν=0.3)

A=EmaxρS(t)3ν(12ν)(1+ν);B=EmaxρS(t)32(1+ν)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

T˜σ˜:=ψmechSε˜=Aε˜+2Btr(ε˜)ISE120

The unknowns of the remodeling problem are the normal velocity of the bone surfaceVN(t), the surface density of minerals ρS(t)and its superficial concentration. We shall herewith simulate the resorption of a hollow bone submitted to a composite applied stress, consisting of the superposition of an axial and a radial component, as

σ=σrrerer+σzzezezE121

in the cylindrical basis(er,eθ,ez); this applied stress generates a preexisting homogeneous stress state within the bulk material, inducing a surface stress given by

σ˜=P.σ=σzzezezE122

The radial component of Eshelby stress Σrris 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:

ϒ˜ggN=1ri(t){18Δμnσ(t)+A+2B8(A+B)Bσzz2}E123

with the material coefficients A,Bgiven in (5.12), and the axial stress σzzpossibly function of time. A simple linear relation of the velocity of the growing surface to the driving force is selected, viz

VN(t)=Cϒ˜ggN(t)E124

with Ca 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=ρgVN/M=108mol.s1.m2VN=JM/ρg3.3.1012m/s0.286μm/dayE125

selecting a molar massM1.004kg/mol, following (Silva and Ulm, 2002). This value is an initial condition for the radius evolution (its rate is prescribed), leading toC=3.5.1023m2.kg1.s; it is however much lower compared to typical values of the bulk growth velocity, about10μm/day.

The mass balance equation for the surface density of minerals ρSwrites

ρ˙S+ρSS.V˜=ΓSρSE126

expressing as the following conservation law

ρ˙SρSVNri(t)=Γ0SρS(t)=ρs0r0ri(t)exp(Γ0St)E127

The initial surface density of mineralsρs0=ρS(t0), is evaluated from the bulk density of HA, viz3000kg/m3, and the estimated thickness of the attachment region of osteoclasts, about 7μm(Blair, 1998), henceρs02.1.102kg/m2.

The surface growth rate of mass Γ0Sis 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): Γ0Sis estimated by considering that 80% of the superficial minerals have been dissolved in a 2 months period, henceΓ0S2.2.107s1. 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

[Ca3(PO4)2]3Ca(OH)2+8H+10Ca2++6HPO42+2H2OE128

The kinetic equation is chosen as:

nσ(t)t=γ˜ρs(t)nσ(t)=γ˜ri(t)ρs0r0exp(Γ0St)nσ(t)E129

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

5.3 Simulation results

The present model involves a dependency of the triplet of variables {ri(t),ρS(t),nσ(t)}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 n0σis taken as unity, vizn0σ=1mol.m3.

  • The initial radius r0:=ri(0)is estimated as r0=1.6cmfor 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,σ{1MPa,2MPa,5MPa}. The surface velocity (Figure 7) shows an acceleration of the resorption process with time, which is enhanced by the stress level, as expected from the higher magnitude of the driving force.

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

An order of magnitude of the simulated radial surface velocity is about 10μm/dayfor 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σzzcrit(t), depending upon the density and concentration, given from (5.18), (5.19) as

σzzcrit(t)9.4.1010ρS(t)3/2nσ(t)1/2E130

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: σzz=1MPa(thick line), σzz=2MPa(dashed line), σzz=5MPa(dash-dotted line).

Figure 8.

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

For an applied stress σzz=0.2MPalying 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:σ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.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Jean-François Ganghoffer (November 2nd 2011). Thermodynamics of Surface Growth with Application to Bone Remodeling, Thermodynamics - Interaction Studies - Solids, Liquids and Gases, Juan Carlos Moreno-Pirajan, IntechOpen, DOI: 10.5772/19978. Available from:

chapter statistics

1596total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Thermodynamic Aspects of CVD Crystallization of Refractory Metals and Their Alloys

By Yu. V. Lakhotkin

Related Book

First chapter

Thermodynamics of Molecular Recognition by Calorimetry

By Luis García-Fuentes, Ramiro, Téllez-Sanz, Indalecio Quesada-Soriano and Carmen Barón

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us