InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Chemical Engineering » "Thermodynamics - Interaction Studies - Solids, Liquids and Gases", book edited by Juan Carlos Moreno-Pirajan, ISBN 978-953-307-563-1, Published: November 2, 2011 under CC BY 3.0 license. © The Author(s).

Chapter 14

Thermodynamics of Surface Growth with Application to Bone Remodeling

By Jean-François Ganghoffer
DOI: 10.5772/19978

Article top


Domain variation due to the virtual velocity field w
Figure 1. Domain variation due to the virtual velocity field w
Formation of an interface from a fixed number of moles of α and β.
Figure 2. Formation of an interface from a fixed number of moles of α and β.
The broken bond model for surface energy
Figure 3. The broken bond model for surface energy
Wulff’s construction to calculate the minimizing surface for a fixed volume
Figure 4. Wulff’s construction to calculate the minimizing surface for a fixed volume
Tissue element under growth: elements of differential geometry.
Figure 5. Tissue element under growth: elements of differential geometry.
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].
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].
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 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).
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).
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).
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.
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.

Thermodynamics of Surface Growth with Application to Bone Remodeling

Jean-François Ganghoffer1

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


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)


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


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


Observe that the total flux of mass is the sum of a convective flux ρku and 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ρk writes (Vidal et al., 1994) accounting for the relationj=1nJk=j=1nρkuk=0, as


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 mk the 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 nk satisfies the balance equation


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


The two previous equalities enter into Gibbs relation as


with θ the temperature and μ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 ξ˙j


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


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


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


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


Hence, the rate of the entropy density decomposes into


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


and of the internal entropy production


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


One has assumed in this alternative that the mechanical power ρw˙=σeq:u does not include a flux contribution, hence only the heat diffusion contributes to the flux of internal energy. The contribution σ:ε˙+k(μk/V)n˙k is 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


Hence, the internal entropy production is identified as


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


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 Ωt expresses then as


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


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


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


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


with w the 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


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


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


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


The strong form of the balance law of mass writes finally


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


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


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


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)


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

Expressing the total mass of the domain Ωt asm(Ω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


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


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ρJdX also 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)


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


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


with r the volumetric heat supply (generated by growth), and q the 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


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)


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


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


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


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


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


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


with X,x the Lagrangian end Eulerian positions in the referential and actual configurations denoted ΩR,Ωt respectively, as the product of the growth deformation gradient Fg and the growth accommodation mapping Fa

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

Hence Jg<1 describes growth, whereas Jg>1 represents 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

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 s the density of internal energy and entropy per unit mass respectively; the first and second principles of thermodynamics write (Munster, 1970)


with Jq the 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 σs the entropy production, always positive (it is dissipated). Introducing the free energy density per unit massψ:=uθs, with s the entropy density, we then immediately obtain the rate of variation of the free energy density


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


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:


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:


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


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


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


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)


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


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

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


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


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


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


with zσ:=Zσ/A the 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


For the whole system, using the previous decomposition

one has


The variation of the free energy is


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




This leads to the differential


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


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


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


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,

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

with e the thickness of the surface; in most cases, the parameter e is 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 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μ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


with Sσ the entropy of the phaseσ. Consider next a new system for which T and μi are 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


Subtracting both variations of Helmholtz free energy by unit surface gives


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


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


and thus finally


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


The second order tensor ε˜ij is the strain (a small perturbation scheme is presently adopted) induced by the component σ˜ij acting in the jth direction per unit length of the edge normal to the ith direction, with both indices i, j lying in the plane of the surface. Previous equation is valid for an anisotropic solid, and reduces in the case of an isotropic surface to the previously written Couchmann-Everett formula, with g half the trace of the surface stress tensor. The proof of previous formula follows Mullins derivation: let imagine a unit cube with edges parallel to the 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 W0 the 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 dx in the x direction, at fixed y edge, but varying z edge. Let W1 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, henceW0+W2=W1+W3, viz


The difference W1W0 is 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


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


Due to the equalities


it finally results


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


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


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

Lastly, consider a square section in the xy plane of side (Aσ)1/2 and 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


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


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


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)


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


which may be interpreted from an energetic point of view as follows: the term VγdP is the mechanical work done against the external force field, the contribution SγdT represents the heat of formation of the interface, and γ*dA is the mechanical work done against the internal force field of both phases α,β by motion of the molecules from the bulk to generate a new interface. The excess free energy of formation of the interface, 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

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


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


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


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+γ*dA and substitution of the previous interfacial Maxwell relations, thus


It immediately results the specific interfacial excess enthalpy


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θ/a in the x and y directions respectively (figure 3), with a the distance to nearest neighbour (function of the type of atomic packing) and θ the inclination of the overall crystal shape resulting from the total number of steps being created.

The surface energy is given by the expression


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

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 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 γ(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


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 N toSg, and possibly explicitly upon the surface position vector XS on Sg (no tilda notation is adopted here since the support of XS is strictly restricted to the surfaceSg), and chemical energyμknkσ, with μk the 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

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


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


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 ψS can be expanded using the surface divergence theorem (equality (3.15) in Ganghoffer, 2010a) as


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


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


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


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


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


with VN the normal surface velocity, Mthe bone mineral molar mass, and JρgVN/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 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):


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

with ρS(t) the surface density of mineral, Emax=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 ψ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


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,B mechanical properties of the surface, expressing versus the surface density of minerals and the maximum value of the traction modulus as (the Poisson ratio is selected asν=0.3)


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


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


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


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


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


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


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


expressing as the following conservation law


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 Γ0S 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): Γ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


The kinetic equation is chosen as:


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.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,σ{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/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σzzcrit(t), depending upon the density and concentration, given from (5.18), (5.19) as


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.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:σ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.


1 - H. C. Blair, 1998 How the osteoclast degrades bone. Bioessays, 20 837846
2 - D. R. . Carter, T. E. Orr, D. P. Fyrhie, 1989 Relationships between loading history and femoral cancellous bone architecture. J. Biomech., 22 3 231244
3 - P. R. Couchman, R. G. Linford, 1980 Aspects of solid surface thermodynamics: relationships of the Shuttleworth type. Jnal of Electroanalytical Chemistry. 115 143148
4 - S. C. Cowin, 2001 Bone Mechanics Handbook. CRC Press
5 - J. D. Eshelby, 1951 The force on an elastic singularity. Phil. Trans. R. Soc., A244 87112
6 - H. M. Frost, 1964 Dynamics of bone remodeling. Bone Biodynam., 2 315333
7 - J. F. Ganghoffer, B. Haussy, 2005 Mechanical modeling of growth considering domain variation. Part I: constitutive framework. Int. J. Solids Struct., 42 15 43114337
8 - J. F. Ganghoffer, 2010a Mechanical modeling of growth considering domain variation-Part II: Volumetric and surface growth involving Eshelby tensors. J. Mech. Phys. Solids, 58 9 14341459
9 - J. F. Ganghoffer, 2010b On Eshelby tensors in the context of the thermodynamics of open systems: Application to volumetric growth. Int. J. Engng Sci., doi:10.1016/j.ijengsci.2010.04.003
10 - K. Garikipati, 2009 The kinematics of biological growth. Appl. Mech. Rev. 62 (3), 030801 EOF , doi:10.1115/1.3090829
11 - F. C. Goodrich, 1969 In Surface and Colloid Science, 1 Ed. By E. Matijevic, 1
12 - M. E. Gurtin, A. I. Murdoch, 1975 A continuum theory of elastic material Surfaces. Archive for Rational Mechanics and Analysis, 57 4 291323
13 - M. E. Gurtin, A. Struthers, 1990 Multiphase thermomechanics with interfacial structure. Part 3. Arch. Rat. Mech. Anal., 112 97160
14 - M. E. Gurtin, 1995 On the nature of configurational forces. Arch. Rat. Mech. Anal., 131 67100
15 - M. E. Gurtin, 2000 In: Configurational forces as Basic Concepts of Continuum Physics. Springer, New York
16 - R. Huiskes, T. J. Sloff, 1981 Geometrical and mechanical properties of the human femur. Biomechanics VII-A: 7th international congress on biomechanics, Ed. A. Morecki, 3A 5764
17 - R. Huiskes, H. Weinans, H. J. Grootenboer, M. Dalstra, B. Fudala, T. J. Slooff, 1987 J. Biomech., 20 11351150
18 - D. A. Kessler, 1990 Coupled-map lattice model for crystal growth. Phys. Rev., Vol.A 42, 61256128
19 - D. Kondepudi, I. Prigogine, 1998 Modern Thermodynamics: From Heat Engines to Dissipative Structures. Wiley
20 - J. S. Langer, 1980 Instabilities and pattern formation in crystal growth. Rev. Mod. Phys., 52 128
21 - P. H. Leo, R. F. Sekerka, 1989 The effect of surface stress on crystal-melt and crystal-crystal equilibrium. Acta Metall., 37 12 31193138
22 - R. G. Linford, 1973 Surface thermodynamics of solids. Solid State. Surface Science, 2 1152
23 - T. W. Lu, J. J. O’Connor, S. J. G. Taylor, P. S. Walker, 1997 Influence of muscle activity on the forces in the femur: an in vivo study. J. Biomech. 30 11 11011106
24 - G. A. Maugin, C. Trimarco, 1995 The dynamics of configurational forces at phase-transition fronts. Meccanica 30 605619
25 - G. A. Maugin, 1993 Material Inhomogeneities in elasticity. Chapman et al., London
26 - R. D. Mindlin, 1965 Second gradient of strain and surface-tension in linear elasticity. Int. J. Solids Struct., 1 417438
27 - P. Muller, R. Kern, 2001 In Stress and Strain in epitaxy: theoretical concept measurements. Ed. M. Handbucken, J.P. Deville, Elsevier
28 - W. W. Mullins, 1963 in Metals surfaces, structures, energetics, and kinetics. Amer. Soc. Metals. Metals Park, Ohio, 7
29 - A. Munster, 1970 Classical Thermodynamics, John Wiley and Sons
30 - D. T. Reilly, A. H. Burstein, 1975 The elastic and ultimate properties of compact bone tissue, J. Biomech., 8 393405
31 - J. C. Rice, S. C. Cowin, J. A. Bowman, 1988 On the dependence of the elasticity and strength of cancellous bone on apparent density. J. Biomech., 21 155168
32 - R. Ruimerman, P. Hilbers, B. Van Rietbergen, R. Huiskes, 2005 A theoretical framework for strain-related trabecular bone maintenance and adaptation. J. Biomech., 38 931941
33 - E. C. C. M. Silva, F. J. Ulm, 2002 A bio-chemo-mechanics approach to bone resorption and fracture. Proc. 15th ASCE Engineering Mechanics Confererence. 0205 /06/2002, Columbia University, New York, USA
34 - R. Skalak, G. Dasgupta, M. Moss, 1982 Analytical description of growth. J. Theor. Biol., 94 555577
35 - R. Skalak, D. A. Farrow, A. Hoger, 1997 Kinematics of surface growth. J. Math. Biol., 35 869907
36 - L. Taber, 1995 Biomechanics of growth, remodeling and morphogenesis. Appl. Mech. Rev., 48 487545
37 - D. W. Thompson, 1992 On Growth and Form. Dover reprint of 1942. 2nd edition
38 - C. Vidal, P. Dewel, P. Borckmans, 1994 Au-delà de l’équilibre. Hermann
39 - H. Weinans, R. Huskes, H. J. Grootenboer, 1992 The behavior of adaptive bone remodeling simulation models. J. Biomech., 25 14251441
40 - G. Yang, J. Kaber, B. Van Rietbergen, A. Odgaard, R. Huiskes, S. C. Cowin, 1999 The anisotropic Hooke’s lax for cancellous bone and wood. J. Elast., 53 125146