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. 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 with volumetric density satisfies a prototype balance law of the form
with the flux density of and the local production (or destruction) of. 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, obey the local balance law (Vidal et al., 1994)
with the local barycentric velocity, the molar mass, and the stoechiometric coefficients in the reaction, such that the variation of the mass of the species k due to chemical reactions expresses as
wherein denotes the degree of advancement of reaction. The molar masses satisfy the global conservation law (due to Lavoisier)
Observe that the total flux of mass is the sum of a convective flux and a diffusive flux; the mass production is identified as the contribution. In this viewpoint, the system is in fact closed, since the balance law satisfied by the global density writes (Vidal et al., 1994) accounting for the relation, 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, , with the mass of the same constituent. The molar concentration is defined as, its inverse being called the partial molar volume. The partial mole number satisfies the balance equation
with the flux of species k and 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 the chemical potential of constituent k. The chemical affinity in the sense of De Donder is defined as the force conjugated to the rate
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 the heat flux, and the term 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 does not include a flux contribution, hence only the heat diffusion contributes to the flux of internal energy. The contribution is identified to the term. Previous equality combined with the second principle, equality (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 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; a convective term is also added, corresponding to the transport of nutrients by the velocity field of the underlying continuum. For any quantity a, the convective flux is locally defined in terms of its surface density as; the overall convective flux of a across the closed surface expresses then as
(the minus accounts for the unit exterior normal). This diffusive flux corresponds to a macroscopic flux
The density of microscopic flux 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.
The convective derivative along the vector field of the field writes
In the case coincides with the velocity of the material particles, previous relation delivers the definition of the material (or particular) derivative
The derivative of the volume integral is next calculated, according to Leibniz rule:
with the velocity field of the points on, which is associated to a variation of the domain occupied by the material points of the growing solid body (Figure 1).
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, and the third contribution to diffusion through the boundary of the moving volume. 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 velocity, hence
Mass balance: the mass balance equation is deduced from the identification, 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 the actual density, the physical source of mass, and the scalar physical mass flux across the boundary, projection of the flux (vector). The previous balance law is quite general, as we account for both the variation of the integration volume through the term, and for the source and flux of mass reflected by the right hand side of (28). Localization of previous integral equation gives
with 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), obtained using the mass balance (28)
with 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 as, 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, and to a source term due to growth, hence (see (28))
The last equality is nothing else than - a consequence of (28) - expressed in material format, with the identifications,. The global form also fits within the general balance law for an open system, relation (30) with.
Balance of momentum: the Eulerian version of the balance of momentum writes (Epstein and Maugin, 2000)
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 the specific internal energy density, one may write the energy balance in the actual configuration as
with the volumetric heat supply (generated by growth), and the heat flux across. This writing of the energy balance can be simplified using the balance of kinetic energy with volumetric density, obtained by multiplying (33) by the velocity and integrating over, hence
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
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 the Lagrangian end Eulerian positions in the referential and actual configurations denoted respectively, as the product of the growth deformation gradient and the growth accommodation mapping
The transformation gradients define the mappings of the tangent spaces to the various configurations. The Jacobean of the growth mapping informs about the nature of growth:
Hence describes growth, whereas represents resorption. Growth essentially occurs between the referential and the actual configurations.
Adopting the framework of hyperelasticity, the first Piola-Kirchhoff stress expresses from the strain energy density per unit volume in the reference configuration 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 and the density of internal energy and entropy per unit mass respectively; the first and second principles of thermodynamics write (Munster, 1970);
with the heat diffusion flux, the total entropy flux, the diffusion flux of the k-specie, an external force acting on the k-specie, and the entropy production, always positive (it is dissipated). Introducing the free energy density per unit mass, with the entropy density, we then immediately obtain the rate of variation of the free energy density
The positivity of the entropy production in previous inequality then expresses as
The principle of virtual power (is the kinetic energy, being the virtual power of external and internal forces respectively), leads to the global form of previous inequality in Eulerian format:
with and 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 gradient, the concentration of chemical specie and 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 and to a source term due to growth, hence
The last equality is nothing else than expressed in material format, identifying,. Accordingly, (3.9) becomes
Since previous equality must hold true for arbitrary variations, 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 toin (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.
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 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 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 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 as, with the area of the interface. Any extensive quantity can be decomposed as
with the superficial excess quantity. Regarding surface quantities, one makes a distinction between:
The superficial energy - 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
The variation of the free energy is
Hence, is defined as the partial derivative. 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 relation.
The reversible work needed to form a unit area of new surface is defined at constant temperature and pressure as the partial derivative of the Gibbs free energy of the entire system (bulk phases and surface), quantity, with respect to the formed area, at constant temperature, pressure, and number of moles of each component, viz
whereby the multiplicative factors of the differential elements on the right-hand side of 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, denoted, and the change per unit area of surface created from the surrounding bulk phases, evaluated as the sum over all components,
with the thickness of the surface; in most cases, the parameter is small, and one may neglect the contribution, hence one has the identification. 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 and as a chemical reservoir for each component; accordingly, the components concentrations can be adjusted to maintain the chemical potentials at fixed specified values. Imagine a modification of the temperature by, and of the ith chemical potential by, at fixed surface area; particles from the bulk will enter the solid phase from the bulk, and ithechange of Helmholtz free energy will be
with the entropy of the phase. Consider next a new system for which and are returned to their initial values, but with the surface area increased by, and modify thereafter the temperature by, and the ith chemical potential by; 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 energy, the specific surface entropy, and the surface excessleads to
But one can also express the specific surface Helmholtz energy as, 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. It is related to by
The second order tensor is the strain (a small perturbation scheme is presently adopted) induced by the component 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 half the trace of the surface stress tensor. The proof of previous formula follows Mullins derivation: let imagine a unit cube with edges parallel to the axes, and perform two distinct operations on it:
Stretch the cube reversibly along the x axis by an amount, 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 the work expanded in this transformation. Let next separate the stretched cube along the xy plane, requiring the work, with the variation of the specific surface work due to the stretch (factor 2 arises since two surfaces are created, and the factor since the specific surface work applies per unit surface area).
Separate the original unit cube into two parts along an xy plane, requiring the work, and stretch each half by in the x direction, at fixed y edge, but varying z edge. Let be the work expanded in this operation. The final configuration is the same as that obtained in the first process, hence the same total work has been expanded, hence, viz
The difference is the work due to the stretching operations of both processes, and can be equalized to the x-component of a force in the newly formed surface times the distance through which this force acts, hence
The strain (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 and imagine an extension of the x edge by; the required work is. Extend next the y edge by, 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 equality. 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 is
which may be interpreted from an energetic point of view as follows: the term is the mechanical work done against the external force field, the contribution represents the heat of formation of the interface, and is the mechanical work done against the internal force field of both phases by motion of the molecules from the bulk to generate a new interface. The excess free energy of formation of the interface, potential, 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 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 as, but with addition of the new mass such that
As a compensation for the volume change accompanying the formation of the interface; hence, 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 allows writing relations between partial derivatives as the analogues of the bulk phase Maxwell relations;.
The introduced quantities are respectively the interfacial excess entropy and the specific interfacial excess volume; the compact notation stands for the two quantities.
The specific interfacial excess energy is obtained by simply integrating the differential at constant pressure and temperature, introducing the specific interfacial excess energy. Last relation implies that the temperature and pressure dependence of can be determined from those of. The specific interfacial excess energy is obtained from a Legendre transform to and substitution of the previous interfacial Maxwell relations, thus
It immediately results the specific interfacial excess enthalpy
with 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.
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 and in the x and y directions respectively (figure 3), with 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 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 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 constraint. 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, with 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).
with anisotropic surface tension
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 (surface density) and to line densities defined as the projections onto the unit vectors resp. along the contour of the open growing surface (Figure 5); hence, those line densities are respectively tangential and normal to the surface (forces acting in the tangent plane).
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 depending upon the transformation gradient, a surface energy with density per unit reference surface, depending upon the surface gradient, the unit normal vector to, and possibly explicitly upon the surface position vector on (no tilda notation is adopted here since the support of is strictly restricted to the surface), and chemical energy, with the chemical potential of the surface concentration of species. The surface gradient maps material lengths (or material tangent vectors) onto the deformed surface; it is elaborated as the surface projection of (onto the tangent plane to), viz
The tissue element under grow is submitted to a surface force field (surface density) and to line densities defined as the projections onto the unit vectors resp. along the contour of the open growing surface (Figure 5); Hence, those line densities are respectively tangential and normal to the surface (forces acting in the tangent plane).
The variation of the previously built potential energy of the growing tissue elementis next evaluated, assuming applied forces act as dead loads, using the fact that the variation is performed over a changing domain (Petryk and Mroz, 1986), here the growing surface. 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) 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 (). Observe that the volumetric Eshelby stress triggers surface growth in the sense of the boundary values taken by the normal Eshelby-like traction. The variation of the surface energy contribution 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 stress. The Lagrangian curvature tensor is defined as. 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 as the vector acting on the variation of the surface position
itself built from the surface stress, and on the curvature tensor in the referential configuration.
5.2. Bone remodeling
Bone is considered as a homogeneous single phase continuum material; from a microstructural viewpoint, bone consists mainly of hydroxyapatite, a type-I collagen, providing the structural rigidity. The collageneous fraction will be discarded, as the mineral carries most of the strain energy (Silva and Ulm, 2002). The ultrastructure may be considered as a continuum, subjected to a portion of its boundary to the chemical activity generated by osteoclasts, generating an overall change of mass of the solid (the mineral fraction) given by
The quantity therein represents the molar flux of bone material being dissolved, hence
with the normal surface velocity, the bone mineral molar mass, and the molar influx of minerals (positive in case of bone apposition, and negative when resorption occurs). Clearly, the previous expression shows that the knowledge of the normal surface growth velocity determines the molar influx of minerals. Estimates of the order of magnitude of the dissolution rate given in (Christoffersen at al., 1997), for a pH of 7.2 (although much higher compared to the pH for which bone resorption takes place) and at a temperature of 310K, are indicative of values of the molar influx in the interval. The osteoclasts responsible for bone resorption attach to the bone surface, remove the collageneous fraction of the material by transport phenomena, and diffuse within the material. This osteoclasts activity occurs at a typical scale of about, 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, 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 (biochemical activity generated by the osteoclasts) with the chemical potential of the mineral of the solid phase, denoted; it can be estimated from the change of activity of the cation (Silva and Ulm, 2002):
This chemical driving force is the affinity conjugated to the superficial concentration of minerals, denoted in the sequel. The conversion to mechanical units is done, considering a density of HA (5.1), hence, 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
Following the representation theorems for isotropic scalar valued functions of tensorial arguments, the surface strain energy density 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 and of its square. The following simple form depending on the second invariant of the linearized part of is selected, adopting the small strain framework, viz, hence
with the surface strain (induced by the existing volumetric strain), and 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)
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 surface, the surface density of minerals 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; 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 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 given in (5.12), and the axial stress possibly function of time. A simple linear relation of the velocity of the growing surface to the driving force is selected, viz
with 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 mass, following (Silva and Ulm, 2002). This value is an initial condition for the radius evolution (its rate is prescribed), leading to; it is however much lower compared to typical values of the bulk growth velocity, about.
The mass balance equation for the surface density of minerals writes
expressing as the following conservation law
The initial surface density of minerals, is evaluated from the bulk density of HA, viz, and the estimated thickness of the attachment region of osteoclasts, about (Blair, 1998), hence.
The surface growth rate of mass 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): is estimated by considering that 80% of the superficial minerals have been dissolved in a 2 months period, hence. 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 (Hankermeyer et al., 2002).
5.3 Simulation results
The present model involves a dependency of the triplet of variables 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 is taken as unity, viz.
The initial radius is estimated as 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,. 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 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, 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.
For an applied stress 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).
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.