Electro-mechanics of bone remodelling

International Journal of Engineering Science 46 (2008) 1173–1182

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Electro-mechanics of bone remodelling Salah Ramtani * LPMTM – CNRS/UPR9001, Université Paris 13, IUT de Saint-Denis, Institut Galilée, 99, Avenue J.B. Clément, 93430 Villetaneuse, France

a r t i c l e

i n f o

Article history: Received 15 January 2007 Received in revised form 30 May 2008 Accepted 3 June 2008 Available online 26 July 2008 Communicated by Prof. K.R. Rajagopal Keywords: Electro-mechanics Bone remodelling Continuum damage and thermodynamics of interacting continua

a b s t r a c t Clinically, it is well known that an electromagnetic field applied to bones hastens the healing process. This is the main motivation behind the present theoretical contribution concerned with the introduction of damage in the electromechanical model previously stated by Demiray and Dost [H. Demiray, S. Dost, The effect of quadrupole on bone remodelling, Int. J. Eng. Sci. 3 (1996) 257–268]. The interaction between electro-mechanical and chemical fields are known to play an important role in the regeneration, growth and maintenance of our skeleton, as well as damage that occurs naturally when the bone is overloaded during day-to-day activities or by cellular activity of resorption. As a result, the propagation of a longitudinal wave in wet bone is considered as an illustrative example. It is shown that in addition to given electric field, both damage and damage gradient affects the reactive mass density of the bone solid matrix. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The nature of the dependence between the form of bones and the load they carry was described by Wolff [3,4] who stated that ‘‘every change in the form or function of a living bone is followed by adaptive changes in its internal architecture and its external shape”. Since the above description, living bone is viewed as a dynamical system which is continuously adapting and constantly repairing itself in response to repetitive external applied loads [5–11]. This adaptation process, including modelling and remodelling mechanisms, has an enormous effect on the overall behavior and health of the entire body. For example, it is well argued that normal, daily activities generate damage in the form of tiny micro cracks throughout bone tissue [12–19] and that have been implicated in physiological phenomena such as remodelling and adaptation. It is thought that the damage process begins at the ultra structural level and evolves through each hierarchical level of bone organization manifesting as micro cracks observable by light microscopy [20–22]. Justus and Luft [23] have demonstrated that straining bone will increase the calcium concentration in the interstitial fluid, probably due to a change in the solubility of the hydroxyapatite crystals and also stated this mechano-chemical effect with inorganic crystals in vitro. Then, the presence of micro cracks may also be regarded as a break in the mineralized matrix which exposes new surfaces lined by calcium ions [8,20]. Electricity is ubiquitous in living things. In humans, it provides the basis for thoughts, senses, movement, and the rhythm of our hearts. As we have been learning over the last 50 years or so, it may also play a crucial role in the functioning of our skeletal system. Electrical properties of bone are relevant not only as an hypothesized feedback mechanism for bone adaptation and remodelling, but also in the context of external electrical stimulation of bone in order to aid it’s healing and repair [24–27].

* Tel.: +33 (0)1 49 40 39 53x61 65; fax: +33 (0)1 49 40 39 38. E-mail address: [email protected] 0020-7225/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2008.06.001

1174

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

Fukada and Yasuda [28] has first demonstrated that dry bone is piezoelectric in the classical sense, i.e. mechanical stress results in electric polarization (the indirect effect); and an applied electric field causes strain (the converse effect). Since that time, many others have confirmed the capacity of bones to produce piezoelectric potentials [29–43]. What is well known is that bone matrix may be considered as an amphoteric ion exchanger with both positive and negative fixed charges, and composed in large part by: (a) inorganic hydroxyapatite crystals which are not asymmetrical and thus cannot contribute to any piezoelectric effect, and (b) collagen that offers an asymmetric structure and believed to be capable of producing piezoelectricity [37]. The magnitude of the piezoelectric sensitivity coefficients of bone depends upon frequency, direction of load and relative humidity [35–38,41,43,44]. Theoretical analyses of bone piezoelectricity [45–49] may be relevant to the issue of bone remodelling. Recent, thorough, studies have explored electromechanical effects in wet and dry bone [38,47]. They suggest that two different mechanisms are responsible for these effects: classical piezoelectricity mainly due to the molecular asymmetry of collagen in dry bone, and streaming potentials found in moist or living bone, and generated by the flow of a liquid across charged surfaces. The second mechanism has been argued by dielectric measurements and has suggested that the electromechanical effect in wet (fluid saturated) bone is not due to a piezoelectric effect [46,49]. Thus, fluid transport plays a significant role not only in various aspects of bone metabolism such as mineralization, but also in the electromechanical properties of bone [48,49]. From the mechanical point of view, the bone can be considered as a composite material with interconnected pores containing both a solid piezoelectric electric bone matrix and a fluid with independent mass and charge densities [1,2]. In particular, the following assumptions are made [1–4,51–53]: (a) The solid component is composed of bone cells embedded in an extra cellular matrix that includes collagen for elasticity and hydroxyapatite crystals for strength. The fluid component contains blood and extra cellular fluid (perfusant), (b) the transfer of mass, such as blood salts Ca3(PO4)2, occurs as the result of the interaction between the biochemical and the electro-mechanical reaction which are mediated by bone cells. Energy and entropy may also be transferred to or from the matrix structure by the above complex process, (c) the blood plasma, which supplies the materials necessary for the synthesis of bone matrix, is always in contact with the extra cellular fluid, (d) the characteristic time of chemical reactions are several orders of magnitude greater than the characteristic time associated with a complete perfusion of the blood plasma. Hence, any excess heat generated by the chemical reactions is quickly carried out by the circulation, (e) as the remodelling process is known as very slow, magnetic effects have been neglected. Following the reference paper of Demiray and Dost [1,2], the quadrupole effects as well as polarization effect of the bone solid matrix are also taken into account. This work is mainly concerned with the introduction of damage in the electromechanical model previously stated by Demiray and Dost [1]. In the light of the newly introduced damage variable, the propagation of a longitudinal wave in wet bone, also investigated by [1,2], is taken again as an illustrative example. As a result of this study, it is shown that, in addition to applied electric field, both damage and damage gradient affects the reactive mass density of the bone solid matrix. Mechanics of interacting continua (Mixture theory) is a suitable framework to describe the coupled interaction between fluid flow, electricity, continuum micro-damage and solid deformation in porous materials, and is a common modelling approach in both engineering and biology. The reader could consult [1,2,50,53–60,69] for understanding details which are voluntarily omitted herein.

2. Kinematics and field equations 2.1. Kinematics We shall consider a mixture of two continua, a solid piezoelectric bone matrix S(s) and a charged conducting fluid S(f). It is assumed that at any time t, each place in the region of space is occupied by particles belonging to both S(s) and S(f). The motion of the solid and fluid is represented by the mappings

x ¼ xðaÞ ðXðaÞ ; tÞ;

ða ¼ f ; sÞ;

ð1Þ

where, X a denote the reference positions of typical particles of S(s) and S(f). These motions are assumed to be one-to-one, ( )

continuous, and invertible. The velocity vectors associated with these motions are

   ox  DðaÞ x  vðaÞ  otðaÞ  ðaÞ  Dt ðaÞ ;  X   o2 x  DðaÞ v   aðaÞ  ot2ðaÞ  ðaÞ  Dt ðaÞ ; X

where

DðaÞ Dt

denotes material time derivative operator following the path of a species

ð2Þ

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

 DðaÞ o ðaÞ o :   þ vk ot x oxk Dt

1175

ð3Þ

The deformation gradient and the Jacobian of the motions associated with the ath components are

  FðaÞ  oxðaÞ ;  oX   J ðaÞ  detðFðaÞ Þ–0:

ð4Þ

2.2. Mass conservation The axiom of balance of mass for a mixture consists of two parts. The first is a statement of balance of mass for each constituent, and the second a statement of balance of mass for the mixture as a whole [1,2,49,55,56]. Then, the local statement of balance of mass for the ath constituent can be written as:

 oq  ðaÞ þ ðq vðaÞ Þ ¼ c ; ðaÞ  ot ðaÞ k ;k  P  c ¼ 0; ð a Þ 

ð5Þ

a¼f;s

where q(a) is the partial mass density and c(a) is the rate of mass transfer from other species into the ath components by electro-chemical interactions. Damage is conventionally defined as the progressive deterioration of materials due to nucleation and growth of micro cracks. At the ultra structural level, damage is hypothesized to be related with debonding of the collagen-hydroxyapatite composite. At the micro structural level, damage is associated with slipping of lamellae along cement lines, cracking along cement lines or lacunae, shear cracking in cross-hatched patterns and progressive failure of the weakest trabeculae. At the macroscopic level, damage is hardly visible before a large crack and subsequent global failure occurs, but may involve a substantial alteration in mechanical function at early stages. In this case different levels of damage are related to the principle directions, and thus a simple scalar damage parameter is no longer sufficient to quantify damage in all directions. Thus, the anisotropic phenomenon of the micro cracks distribution in the material is interpreted using a symmetric second-order damage tensor Dij. ^ ðsÞ for the damaged solid matrix as follows [1] Let define the reactive mass density q

q^ ðsÞ ¼ qðsÞ 

q0ðsÞ J ðsÞ

-ðDij Þ;

ð6Þ

where q0ðsÞ is the initial mass density of the undamaged solid matrix, Dij, -ðDij Þ are respectively, the second order damage tensor and the scalar ‘‘damage function” chosen under the following restrictions:

d6 0: dDij

-ð0Þ ¼ 1; 0 < -ðDij Þ 6 1 and

ð7Þ

Introducing (6) into the local statement of the bone matrix mass conservation (5), one obtained the local form of the reactive mass conservation for the damaged bone matrix component as follows:

^ ðsÞ oq  ^ ðsÞ vðsÞ þ ðq k Þ;k ¼ cðsÞ ; ot

ð81 Þ

where the rate of mass transfer cðsÞ from the fluid into the damage-bone matrix by electro-chemical interactions is stated as

0

0 @ ðsÞ

cðsÞ ¼ cðsÞ  q

1 d- _ Dkl þ J ðsÞ dDkl J ðsÞ

1

!

ðsÞ v k A:

ð82 Þ

;k

One can observe that the rate of reactive mass transfer (81) from the fluid phase into the damage-bone matrix phase is affected by both the damage rate, damage and the damage gradient. 2.3. Momentum conservation The axiom of balance of momentum involves statements about the linear and angular momentum of the constituents as well as of the mixture. Following [1,2,69], the local form of linear momentum is expressed as: ðaÞ

ðaÞ

ðaÞ

ðaÞ

ðaÞ

tkl;k þ g l þ Rl þ cðaÞ wl  qðaÞ al ðaÞ tkl

ðaÞ gl

¼ 0;

ð9Þ ðaÞ Rl

where is the stress tensor, is the electromagnetic body force per unit volume, is the rate of linear momentum P ðaÞ ðaÞ transfer from other species into ath species submitted to the condition a¼f;s Rl ¼ 0, and wl is defined from the reference l  vlðaÞ . l as wðl aÞ ¼ v velocity v

1176

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

Here the effect of mechanical body force have been neglected and the electrical force densities are expressed as:

 ðsÞ  g ¼ E PðsÞ þ E Q ðsÞ ; l;k k l;mn mn  l  ðfÞ  g ¼ qE ; l l

ð10Þ

where E, P and Q are respectively, the electric field density vector at x, the polarization density vector of the solid matrix, the quadrupole density tensor of the solid matrix and q is the charge density of the conducting fluid phase. Here we assumed that all electrical charges are carried by the fluid phase and neglected the polar effect of the fluid [1,2]. The local form of the balance of angular momentum for mixtures leads to the following relation [1,2] ðsÞ

ðfÞ

ðsÞ

ðsÞ

t½kl þ t ½kl þ P½k El þ 2Q ½km El;m ¼ 0;

ð11Þ

where the indices enclosed in brackets are used to denote the skew-symmetric part of the associated tensor, e.g. ðaÞ ðaÞ ðaÞ t½kl ¼ 12 ðtkl  tlk Þ. Eq. (11) shows that the partial stress tensors in the mixture are not symmetric, as opposed to the case of a single continuum. 2.4. Energy conservation With the effect of heat source neglected, the local form of the energy conservation is expressed as [1,2]

1 2

ðaÞ ðaÞ ðaÞ ðaÞ ðaÞ ðaÞ qðaÞ e0ðaÞ þ cðaÞ eðaÞ ¼ tðklaÞ vl;k  hk;k þ cðaÞ wk wk þ Rk wk þ EðaÞ þ eðaÞ ;

ð12Þ

ðaÞ hk

where e(a) is the partial internal energy density, is the partial heat flux, e(a) is the rate of energy transfer submitted to the P condition a¼f;s eðaÞ ¼ 0 and E(a) is the electrical power density of the ath species given by

  EðfÞ ¼ 0;   b ðsÞ ; b ðsÞ þ Ek;m Q  EðsÞ ¼ Ek P k km

ð13Þ

ðsÞ b ðsÞ ðsÞ ðsÞ b ðsÞ are respectively defined as P _ ðsÞ b ðsÞ ¼ P_ ðsÞ þ vðsÞ b ðsÞ and Q P r;r P k , Q kl ¼ Q kl þ vr;r Q kl , and where a prime (or over dot) denotes the k kl k k material time derivative of the corresponding quantity.

2.5. Electromagnetic field equations The local form of the field equations are the Faraday’s Law in which the gradient of the electric field is symmetric for this particular quasi-static case, the Gauss’s Law and the Ampere’s Law expressed respectively as

  eijk Ek;j ¼ 0;   D ¼ q;  i;i  oD  i þ J ¼ 0;

ð14Þ

i

ot

where eijk is the permutation tensor, Di and Ji are respectively the dielectric displacement vector and the total electrical current density given by the following local relations [1]

  J ¼ qvðfÞ ;  i i   Di ¼ e0 Ei þ PðsÞ  Q ðsÞ ; i mi;m

ð15Þ

where e0 is the dielectric constant of the vacuum.

2.6. Entropy inequality Each constituent is assigned an entropy density g(a), a common absolute temperature of the mixture h. Here a single entropy inequality for the whole mixture is used as

X

2

ða Þ

4q g_ ðaÞ þ cðaÞ g þ hk ðaÞ ða Þ h a¼f;s

! 3 5 P 0:

ð16Þ

;k

Introducing the Legendre transformation as

 ðsÞ ðsÞ   WðsÞ  eðsÞ  hgðsÞ  ðEk Pk þEk;l Q kl Þ qðsÞ    WðfÞ  eðfÞ  hgðfÞ

ð17Þ

and eliminating the internal energies e(a) in favor of the Helmoltz free energy W(a), Eq. (16) gives the following inequality

1177

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

_ ðsÞ þ h_ g Þ  q ðW _ ðfÞ þ h_ g Þ þ t ðsÞ vðsÞ þ t ðfÞ vðfÞ  E_ k PðsÞ  E_ k;l Q ðsÞ  hk h;k  cðWðsÞ  WðfÞ Þ þ c wk wk þ Rk wk P 0; qðsÞ ðW ðsÞ ðfÞ ðfÞ kl l;k kl l;k k kl 2 h ð18Þ ðsÞ

ðfÞ

ðsÞ

which must be valid for all thermodynamically admissible processes and were we have set c(s)  c, hk  hk  hk , Rk  Rk ðfÞ ðsÞ and wk  wk  wk the relative velocity of fluid with respect to the solid matrix. 3. Constitutive equations

In this section, we shall state the constitutive assumptions that will serves to classify the particular type of mixture consisting of an electrically polar solid (bony matrix) and an inviscid fluid (blood). For such a mixture, constitutive independent variables and free energy for both the fluid phase W(f) and the damaged solid phase bone matrix W(s) are motivated following previous works [1,2,17,18,35,43,69]

^ ðsÞ ; q; Ek ; h; wk ; Dkl ; F kK ; q

WðfÞ ¼ WðfÞ ðqðfÞ ; Ek ; q; hÞ and hk ¼ 0;

ð19Þ

WðsÞ ¼ WðsÞ ðqðsÞ ; F kK ; Ek ; Dkl ; hÞ:

ð20Þ

One can note that heat flux hk for the whole of mixture is zero, but it may not vanish for the individual species material derivative of (19) and (20) and introducing the result into the entropy inequality (18), we have

ðaÞ hk .

Taking the

" #     oWðfÞ _ oWðfÞ oWðfÞ oWðsÞ _ ðfÞ ðfÞ  qðsÞ gðsÞ þ dkl þ qðfÞ q h  qðfÞ gðfÞ þ h þ tkl þ q2ðfÞ dkl vl;k oh oh oqðfÞ oq       oWðfÞ _ oWðsÞ oWðsÞ oWðsÞ oWðsÞ _ ðsÞ ðsÞ ðsÞ ^ ðsÞ Ek  Q kl þ qðsÞ Ek;l þ t kl  qðsÞ F kK þ qðsÞ q dkl vl;k  PkðsÞ þ qðsÞ þ qðfÞ ^ ðsÞ oF lK oq oEk oEk oEk;l ! " # ! oðqðfÞ WðfÞ Þ oðqðsÞ WðsÞ Þ oWðfÞ oWðsÞ q0ðsÞ o- _ oWðsÞ c ðsÞ þ qðsÞ q0ðsÞ  qðsÞ  Dkl  qðfÞ Ek;l wl þ c  vk þ wk wk ^ ðsÞ JðsÞ ^ ðsÞ 2 oDkl J ðsÞ oDkl oEk oq oqðfÞ oq ;k

þ Rk wk P 0:

ð21Þ D_ kl , vðl;kaÞ ,

_ This inequality is linear in the variables h, E_ k and E_ k;l . Thus, in order that this inequality be valid for all independent variations of these independent constitutive variables, the coefficients of them must vanish. Hence, we obtain

  g ¼  oWðsÞ ; g ¼  oWðfÞ ;  ðsÞ ðfÞ oh oh    oW oW oW  ðsÞ ðsÞ þ qðfÞ oEðfÞ ; Q kl ¼ qðsÞ oE ðsÞ ;  P k ¼  qðsÞ oEðsÞ k k k;l  h i  ðfÞ  t ¼  q2 oWðfÞ þ q q oWðfÞ d ; kl  kl ðfÞ ðfÞ oqðfÞ oq   ðsÞ ^ ðsÞ ooWq^ ðsÞ dkl :  t kl ¼ qðsÞ ooFWðsÞ F kK  qðsÞ q ðsÞ

ð22Þ

lK

Two contributions arise naturally from the remaining parts of the inequality (21) and one can see that in (24) mass transfer interact with damage rate and damage gradient dissipations.

" #   oðqðfÞ WðfÞ Þ oWðfÞ c c P 0; Rk þ wk  qðfÞ El;k wk þ 2 oEl oqðfÞ ! !   oðqðsÞ WðsÞ Þ oWðsÞ q0ðsÞ o- _ oWðsÞ ðsÞ  qðsÞ c P 0:  Dkl þ qðsÞ q0ðsÞ vk  ^ ðsÞ J ðsÞ ^ ðsÞ oDkl J ðsÞ oDkl oq oq

ð23Þ ð24Þ

;k

One can observes also that for an electro-chemically reacting charged mixture the solid component constitutionally has a pressure term due to mass transfer, while the fluid phase has a pressure contributed from the variation of charge density of the fluid continuum [1,2]. The implication of the axiom of objectivity, which states that the constitutive functions should remain form invariant under the full group of transformation of the spatial frame of references, is that the constitutive functions should depend on the following independent set [1,2]

  EKL ¼ 1 ðF kK F kL  dKL Þ; F KL ¼ Ek;l F kK F lL ;  2   EK ¼ Ek F kK ; J E ¼ 12 Ek Ek ; W K ¼ wk F kK   DKL ¼ Dkl F kK F lL ; ^ ðsÞ , h and DKL. as well as variables q, q(f), q

ð25Þ

1178

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

Thus, the free energy density functions have the following form

  WðsÞ ¼ WðsÞ ðEKL; F KL; EK ; q ^ ðsÞ ; DKL ; hÞ;    WðfÞ ¼ WðfÞ ðqðfÞ ; J E ; q; hÞ:

ð26Þ

Introducing (26) into (21) and (22) and making some manipulations that are given in [1], one can obtain the following constitutive equations

 ðsÞ  t ¼ q oWðsÞ F kK F lL  PðsÞ El  2Q ðsÞ El;m  q q ^ oW^ ðsÞ dkl ;  kl ðsÞ oEKL ðsÞ ðsÞ oq k km ðsÞ     ðsÞ oWðfÞ WðsÞ oWðsÞ ðsÞ  P k ¼  qðsÞ ooE ; Q F þ q E ¼  q kL k ðfÞ oJE ðsÞ oF KL F kK F lL ; kl L

ð27Þ

where the presence of damage affect the stress tensor of the solid continuum, the polarization density vector of solid matrix and the solid matrix part of the quadrupole tensor density. 4. Boundary condition One of the outstanding issues, subject of a debate, in mixture theory is the specification of boundary conditions, which has been addressed in [50,62,63] for mixtures in general and for flow of fluids through non-linear elastic solids such as rubber in particular. In a very important paper and with respect to problems that deal with diffusion of fluids through non-linear elastic materials, Rajagopal et al. [50] developed a novel scheme to split the total stress for a class of problems in which the boundary of the mixture is assumed to be in a state of saturation. This additional boundary condition is obtained by the variation of the Gibbs free energy of dilution being zero. With respect to the above pioneer work, Rajagopal [64] considers this aspect as crucial and has recently stated that: ‘‘No discussion of constitutive theory for bodies would be complete which does not discuss the constitutive relations that hold at the boundary of two bodies”. In fact, the basic question is whether conditions such as traction-free surfaces should apply to the individual traction vectors associated with the stress tensors of the solid and fluid phases, or to the traction vector of the mixture. Similar questions exist with regard to the velocity boundary conditions [65]. Massoudi [65] has discussed the importance of and the need for (additional) boundary conditions in Mixture theory and an overview of the model due to Rajagopal and Massoudi [66–68], which is appropriate for the flow of a linearly viscous fluid infused with solid particles, is investigated for additional boundary condition that arises due to higher gradients of density (or volume fraction). The challenging issue of how to ‘split’ the total stress or the total velocity at the boundary is also discussed. An attempt is also made in order to derive a saturation boundary condition. Quite often these boundary conditions are not derived from first principles; instead they are given as ad-hoc assumptions, or they are simply specified as mathematical conveniences. For the present contribution such a task is not simple and heavy. Thus, in our opinion, specific experiments must be used in order to specify these necessary additional boundary conditions before any theoretical attempt. 5. Small strain approximation including scalar damage Damage can be defined in a various number of ways depending on the analysis involved, and there are many physical measures of damage including microscopy (number/density of cracks), ultrasonic waves, density changes, micro-hardness measurements, changes in electrical resistance, and changes in relative amplitude of stress and strain or reduction in elastic modulus during fatigue. The damage law can be: (a) calculated using the remaining life approach, (b) from monotonic loading-unloading experimental test, (c) derived directly by means of a damage criterion, damage consistency condition and loading/unloading conditions. For commodity reasons and in order to state some analytical results, we consider the scalar damage variable case such that previous Eqs. (7) and (81) become

Dkl ¼ ddkl ;

-ðDij Þ ¼ ð1  dÞ and cðsÞ ¼ cðsÞ þ

q0ðsÞ od J ðsÞ

ot

 ðsÞ þ vk d;k :

ð28Þ

In the linear constitutive theory the difference between the material and the spatial coordinates may be disregarded. We have also assumed that the displacement gradient, electric field, the gradient of electric field, deviations in mass and electric charge densities are small and the second order effects may be neglected [1,2]. Thus, the strain energy density functions may be expressed as follows [1,2,49]

  W ¼ ð1  dÞW ðE F E ; q ðsÞ KL; KL; K ^ ðsÞ ; hÞ;  ðsÞ   WðfÞ ¼ WðfÞ ðqðfÞ ; J E ; q; hÞ;

ð291 Þ

where

1 2

1 2

q0ðsÞ WðsÞ ¼ W0 þ a1 q^ ðsÞ þ AK q^ ðsÞ EK þ BKL q^ ðsÞ EKL þ C KL q^ ðsÞ F KL  AKL EK EL  DKLM F KL EM  C KLM EKL EM þ BKLMN EKL EMN 1  C KLMN EKL F MN  DKLMN F KL F MN : 2

ð292 Þ

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

1179

The constitutive equations and material coefficients are restricted by material symmetry conditions, and the bone matrix is known to possess hexagonal polar structure. The reader can refer to the Refs. [1,2,28,38,52,57] for the details that lead to some properties and simplifications. Noting the following relation

F kK ¼ dkK þ dkM U M;K ;

ð30Þ

where dkK is the coordinate shifter and UM,K are the components of the infinitesimal displacement vector, the linearized partial stress tensors and polarization vector take the following form

    ðfÞ  t kl ¼  p0 þ ðqðfÞ  q0ðfÞ Þ @ p00 þ ðq  q0 Þ @@qp0 dkl ; @ qðfÞ  0   ðsÞ  t kl ¼ ð1  dÞ½a1 q ^ ðsÞ dkl þ ðBKL q ^ ðsÞ  C KLM EM þ BKLMN EMN  C KLM F MN ÞdkK dlL ;   b ðsÞ ^ ðsÞ þ AKL EL þ C MNK EMN þ DMNK F MN ÞdkK þ bEk ;  P k ¼ ð1  dÞðAK q   b ðsÞ ^ ðsÞ þ DKLM EM þ C MNKL EMN þ DKLMN F MN ÞdkK dlL ;  Q kl ¼ ð1  dÞðC KL q   @ W @W  p0 ¼ q20ðfÞ @ q ðfÞ þ q0ðfÞ q0 @qðfÞ : 0

ð31Þ

0ðfÞ

As the differences between the material and the spatial coordinates are disregarded in the linear theory, the balance equations are stated from (6), (9), (14) and (15) as

 ðsÞ  t þ R ¼ q0 o2 ul ; oq^ ðsÞ  q0 vk d ¼ c; l  kl;k ðsÞ ot 2 ðsÞ ðsÞ ;k ot   t ðfÞ þ q E  R ¼ q0 ovl ; oq ðfÞ þ q0 vðfÞ ¼ c;  kl;k l 0 l ðfÞ ot ðfÞ k;k ot  ðsÞ  oðe0 Ek þP Þ ðfÞ  ðe E þ PðsÞ Þ ¼ q k ; þ q v ¼ 0; 0 k

k

;k

ot

ð32Þ

0 k

 ¼ q  q0 and /(x, t) is the quasi-static electric potential function from which we derive the electric  ðfÞ ¼ qðfÞ  q0ðfÞ , q where q field vector Ek = /,k. The number of differential equations is sufficient to determine the mechanical and electrical fields completely. The proper boundary and/or initial conditions may be posed for each specific problem of concern. For the linear constitutive relations, the proper forms of Rk and c, with material coefficients submitted to symmetry properties, should be [1]

    ð1Þ ð2Þ ð3Þ  Rk ¼ vð0Þ W M þ vKM EM þ vKMN EMN þ vKMN F MN dkK ; KM    c ¼ c1 q  þ C K W K þ DK EK þ T KL EKL þ PKL F KL : ^ ðsÞ þ c2 q  ðfÞ þ c3 q

ð33Þ

6. Illustrative example and concluding remarks The propagation of a longitudinal wave in wet damaged bone is investigated as an illustrative example with respect to previous work [1] with the following simplifying assumptions:  The mechanical processes in the bone are so slow that the momentum transfer term may be neglected.  The main factors in bone remodelling are the electric field and the strain of the bone matrix, thus, the contribution of other factors, in the first approximation, may be neglected. Therefore, the mass transfer term may be taken to depend linearly on the electric field, the gradient of the electric field and the strain. Furthermore, the contribution of reactive mass density to stress tensors and polarization vector are neglected.  As the remodelling process is known as very slow, the contribution of the damage rate is neglected whereas the damage gradient effect is taken into account. Then, with the aid of above assumptions, the linear theory leads to the following set of constitutive equations

 ðfÞ  t ¼ ½p0 þ a1 q dkl ;  ðfÞ þ a2 q  kl  ðsÞ  t ¼ ð1  dÞ½C klm Em þ Bklmn emn  C klm Em;n ;  kl  ðsÞ  P ¼ ð1  dÞðAkl El þ C mnk r mn þ Dmnk Em;n ÞdkK þ bEk ;  k  ðsÞ  Q kl ¼ ð1  dÞðDklm Em þ C mnkl Emn þ Dklmn Em;n Þ;   R ¼ 0andc ¼ D E þ T E þ P E ; k k k kl kl kl k;l

ð34Þ

where ekl ¼ 12 ðuk;l þ ul;k Þ is the infinitesimal strain tensor. Assuming that all the phenomena occur along the bone axis (x3 = z), in one-dimensional case we have

  u ¼ u ¼ vðfÞ ¼ vðfÞ ¼ 0; u ¼ uðz; tÞ; vðfÞ ¼ vðz; tÞ;  1 2 3 1 2 3  q ¼q ðz; tÞ ^ ðsÞ ðz; tÞ; q  ðfÞ ¼ q  ðfÞ ðz; tÞ; q ^ ðsÞ ¼ q

ð35Þ

1180

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

and the remaining are all zero. Introducing (35) into (34) one can read

 ðfÞ  t ¼ ½p þ a q ; 0 1  ðfÞ þ a2 q  33  ðsÞ  t ¼ ð1  dÞB3333 u0 þ ð1  dÞC 333 /0 þ ð1  dÞC 3333 /00 ;  33  ðsÞ  P ¼ ð1  dÞC 333 u0  ½ð1  dÞA33 þ b/0  ð1  dÞD333 /00 ;  3  ðsÞ  Q 33 ¼ ð1  dÞC 333 u0  ð1  dÞD333 /0  ð1  dÞD3333 /00 ;   c ¼ T u0  D /0  P /00 ; 33 3 33

ð36Þ

where the prime denotes the differentiation with respect to z. Introducing (35) into (32) the following set of ordinary differential equations is stated

 0 0 0 0 ov  a1 q   ðfÞ þ a2 q  q0 / þ qðfÞ ot ¼ 0;  oq 0 00 ðfÞ 0 0 0   ot þ qðfÞ v þ T 33 u  D3 /  P33 / ¼ 0;  ^  oqðsÞ 0 0 00 ¼ 0;  ot  T 33 u0 þ D3 / þ P33 /  q0ðsÞ d ou ot  2  00 000  ð1  dÞ½B3333 u00 þ C 333 / þ C 3333 /   q0ðsÞ oot2u ¼ 0;    ð1  dÞ½C 333 u00  C 3333 u000 þ D3333 /IV   ½e0 þ b þ ð1  dÞA33 /00  q ¼ 0;  h i 000 0  0 00  ð1  dÞ C 333 ou  C 3333 ouot þ D3333 o/ot  ½e0 þ b þ ð1  dÞA33  o/ þ qv0 ¼ 0; ot ot

ð37Þ

^ ðsÞ , q  ðfÞ , q, /, u and v. from which we will determine all the unknown functions q Now, we shall seek a harmonic wave type of solution to the field equations and set

^ ðsÞ ; q  ðfÞ ; q; /; u; vÞ ¼ ðR0ðsÞ ; R0ðfÞ ; Q 0 ; U0 ; U 0 ; V 0 Þ exp½iðkz  xtÞ; ðq

ð38Þ

where x is the angular frequency, k is the wave number, R0ðsÞ ,. . ., V0 are the complex amplitudes to be determined from the field equations and the boundary conditions. Introduction of (38) into the set of differential Eq. (37) yields the following homogeneous algebraic equations

  a1 kR0 þ a2 kQ 0 þ q kU0  q0 xV 0 ¼ 0; 0  ðfÞ ðfÞ   xR0  kðD þ iP kÞU þ T kU þ q0 kV ¼ 0; 3 33 0 33 0 0  ðfÞ ðfÞ   xR0  kðD þ iP kÞU þ ðT k  q0 xd0 ÞU ¼ 0; 3 33 0 33 0  ðsÞ ðsÞ  2  k2 0 2 ðð1  dÞ½C 333 þ iC 3333 kÞU0 þ ðqðsÞ x  ð1  dÞB3333 k ÞU 0 ¼ 0;    2 2 2  k ðe0 þ b þ ð1  dÞ½A33 þ D3333 k ÞU0  k ðð1  dÞ½C 333 þ iC 3333 kÞU 0  Q 0 ¼ 0;  2 2  ix½ikðe þ b þ ð1  dÞ½A þ D 0 33 3333 k ÞU0 þ ðikð1  dÞ½C 333 þ k C 3333 ÞU 0  þ q0 V 0 ¼ 0:

ð39Þ

The Eq. (39)3 expresses the reactive mass density in terms of the displacement of the solid bone matrix, the external electric field and damage. In order to have a non-zero solution to the field quantities, the determinant of the coefficient matrix must vanish. For the specific case where the initial charge density vanishes, i.e. q0 = 0, the dispersion relation reduces to

A0 x4 þ A1 x2 þ A2 ¼ 0;

ð40Þ

where the coefficients A0, A1, A2 are given by

  A0 ¼ q0 ½e0 þ b þ ð1  dÞðA33 þ k2 D3333 Þ; ðsÞ   2 2 2 2  A ¼ k2 f½e þ b þ ð1  dÞðA þ k2 D 0 0 33 3333 Þða1 qðsÞ þ ð1  dÞB3333 Þ þ ð1  dÞ ðC 333 þ k C 3333 Þg;  1  4 2 2 2 2  A ¼ a k f½e þ b þ ð1  dÞðA þ k D 2 1 0 33 3333 Þðð1  dÞB3333 Þ þ ð1  dÞ ðC 333 þ C 3333 Þg:

ð41Þ

In this case we have two waves propagating in the medium: one of them is associated with the fluid and the other with the damaged solid bony matrix. The phase velocity n ¼ xk may be obtained from (40) as

 2  n ¼ a1 ;  1   n2 ¼ ð1dÞB3333 þ  2 q0ðsÞ q0

ðsÞ

ð1dÞ2 ðC 2333 þk2 C 23333 Þ

ðe0 þbþð1dÞðA33 þk2 D3333 ÞÞ

:

ð42Þ

Assuming that the applied external field is electrical origin, i.e. U0 is given, from (39)4 we have 2

U0 ¼

ðC 333 þ ikC 3333 Þk 2

ðq0ðsÞ x2  ð1  dÞk B3333 Þ

Ud ;

where the applied external electrical field is affected by damage as Ud = (1  d)U0.

ð43Þ

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

1181

Thus, from (36)5 and (39)3 both the mass transfer term and reactive mass density are derived, respectively as 3

2

c ¼ ðk P 33  ikD3 Þ/ðz; tÞ þ

q^ ðsÞ ¼

k

x

ik T 33 ðC 333 þ ikC 3333 Þ 2

ðq0ðsÞ x2  ð1  dÞk B3333 Þ

ðD3 þ ikP 33 Þ/ðz; tÞ þ

k

/d ðz; tÞ;

ðC 333 þ ikC 3333 Þk

2

x ðq0ðsÞ x2  ð1  dÞk2 B3333 Þ

ð44Þ

q0ðsÞ od k oz

!  T 33 /d ðz; tÞ;

ð45Þ

where /d(z, t)=(1-d)/(z, t). For a given (applied) external electrical field and measured or supposed distribution of damage, one can note that the increase or the decrease the damage level (repair) lead to: (a) the increase or a decrease of the phase velocity wave in the medium associated with the bone solid matrix without disturbing the fluid phase velocity (Eq. (42)2), (b) the attenuation of the amplitude of the bone matrix displacement (Eq. (43)) and the mass transfer term (Eq. (44)), (c) the attenuation of the reactive mass density affected by both damage and damage gradient (Eq. (45)) which does not appear in the mass transfer term. As a tool for clinical needs, the damage effect can be investigated by comparing the phase velocities (and/or displacement amplitude) of both the ‘‘virgin” bone with respect to the damaged one. Even if the electro-physiological mechanisms involved in the electrical stimulation of bone damage repair remain largely unknown and not roughly investigated, it is well argued that: (a) a physical exercise therapy program designed for bone healing is based on deformations and motions of bones under the application of a mild amount of stress, and (b) an electromagnetic field applied to bones hastens the healing process. These processes involve interactions of electromagnetic fields and mechanical deformations of porous solids, namely bones. Thus, the author believes the theory presented here is relevant to the discussion of the difficult problem of damage-bone remodelling subjected to electromechanical stimulation [59,69], or other biomechanical fields such as cardiac cell electro-mechanics [60,61]. These are some of the motivations behind the present theoretical paper even if experimental data remains, in the end, the final test for the its validity.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

H. Demiray, S. Dost, The effect of quadrupole on bone remodelling, Int. J. Eng. Sci. 3 (1996) 257–268. H. Demiray, Electro-mechanical remodelling of bones, Int. J. Eng. Sci. 9 (1983) 1117–1126. A. Ascenzi, Biomechanics and Galileo Galilei, J. Biomech. 26 (1993) 95–100. J.L. Wolff, The Law of Bone Remodelling, Springer, Berlin, 1986. S.C. Cowin, D.M. Hegedus, Bone remodelling I: theory of adaptive elasticity, J. Elasticity 6 (1976) 313–325. D.H. Hegedus, S.C. Cowin, Bone remodelling II: small strain adaptive elasticity, J. Elasticity 6 (1976) 337–352. R. Huiskes, Simulation of Self-Organization and Functional Adaptation in Bone, Springer-Verlag, Berlin, 1997. T.C. Lee, L. Noelke, G.T. McMahon, J.P. Mulville, D. Taylor, Functional adaptation in bone, in: P. Pedersen, M.P. Bendsoe (Eds.), Synthesis in Bio Solid Mechanics, Kluwer Academic Publishers, Dortrecht, 1998. R.T. Hart, D.T. Davy, Theories of bone modelling and remodelling, in: S.C. Cowin (Ed.), Bone Mechanics, CRC Press, 1989 (Chapter 11). M. Zidi, S. Ramtani, Bone remodelling theory applied to the study of n unit-elements model, J. Biomech. 32 (1999) 743–747. J. Currey, The Mechanical Adaptations of Bones, Princeton University Press, 1984. H.M. Frost, Presence of microscopic cracks in vivo in bone, Bull. Henry Ford Hospital 8 (1960) 25–35. D.R. Carter, W.C. Hayes, Compact bone fatigue damage-I. Residual strength and stiffness, J. Biomech. 10 (1977) 325–337. P.J. Prendergast, D. Taylor, Prediction of bone adaptation using damage accumulation, J. Biomech. 27 (1994) 1067–1076. D.B. Burr, R.B. Martin, M.B. Schaffler, E.L. Radin, Bone remodelling in response to in vivo fatigue microdamage, J. Biomech. 18 (1985) 189–200. D.B. Burr, C. Milgrom, R.D. Boyd, W.L. Higgins, G. Radinel, Experimental stress fractures of the tibia-biological and mechanical aetiology in rabbits, J. Bone Joint Surg. 72B (1990) 370–375. S. Mori, D.B. Burr, Increased intracortical remodelling following fatigue damage, Bone 14 (1993) 103–109. S. Ramtani, M. Zidi, Damaged-bone remodelling theory: thermodynamical approach, Mech. Res. Commun. 26 (1999) 701–708. S. Ramtani, M. Zidi, A theoretical model of the effect of continuum damage on a bone adaptation model, J. Biomech. 34 (2001) 471–479. T.C. Lee, Detection and accumulation of microdamage in bone, M.D. Thesis, University of Dublin, 1997. D.B. Burr, C.H. Turner, P. Naick, M.R. Forwood, R. Pidaparti, En bloc staining of bone under load does not improve dye diffusion into microcracks, J. Biomech. 31 (1998) 285–288. N.L. Fazzalari, M.R. Forwood, B.A. Manthey, K. Smith, P. Holesik, Three-dimensional confocal images of microdamage in cancellous bone, Bone 23 (1998) 373–378. R. Justus, J.H. Luft, A mechanochemical hypothesis for bone remodelling induced by mechanical stress, Calc. Tissue Res. 5 (1970) 222–235. C.A.L. Bassett, A.A. Pilla, R.J. Pawluk, A non-operative salvage of surgically resistant pseudarthrosis and non-unions by pulsing electromagnetic fields: a preliminary report, Clin. Orthop. 124 (1977) 128–143. C.T. Brighton, Z.B. Friedenberg, E.I. Mitchell, R.E. Booth, Treatment of non-union with constant direct current, Clin. Orthop. 124 (1977) 106–123. R.O. Becker, The significance of bioelectric potentials, Bioelectroch. Bioener. 1 (1974) 187–199. C.A.L. Bassett, R.O. Becker, Generation of electric potentials in bone in response to mechanical stress, Science 137 (1962) 1063–1064. E. Fukada, I. Yasuda, On the piezoelectric effect of bone, J. Phys. Soc. (Jpn) 12 (1957) 1158–1162. G. Aschero, F. Gizdulich, F. Mango, M. Romano, Converse piezoelectric effect detected in fresh cow bone, J. Biomechanics. 29 (9) (1996) 1169–1174. Beck, Belinda, The relationship of streaming potential magnitude to strain and periosteal modelling, Ph.D. Dissertation University of Oregon, 1996. D. Gross, W. Williams, Streaming potential and the electromechanical response of physiologically-moist bone, J. Biomechanics. 15 (4) (1982) 277–295. C.T. Hung, F.D. Allen, S.R. Pollack, C.T. Brighton, What is the role of the convective current density in the real-time response of cultured bone cells to fluid flow?, J Biomechanics. 29 (11) (1996) 1403–1409. M.W. Johnson, R.A. Chakkalakal, J.L. Katz, Comparison of the electromechanical effects in wet and dry bone, J. Biomechanics. 15 (1980) 881–885.

1182

S. Ramtani / International Journal of Engineering Science 46 (2008) 1173–1182

[34] F. McDonald, W.J.B. Houston, An in vivo assessment of muscular activity and the importance of electrical phenomena in bone remodelling, J. Anatomy. 172 (1990) 165–175. [35] J. Rubin, K.J. Mcleod, L. Titus, M.S. Nanes, B.D. Catherwood, C.T. Rubin, Formation of osteoclast-like cells is suppressed by low frequency, low intensity electric fields, J. Orthopaed. Res. 14 (1996) 7–15. [36] G. Reinish, Piezoelectric properties of bone as functions of moisture content, Nature 253 (1975) 626–627. [37] J.C. Anderson, C. Eriksson, Electrical properties of wet collagen, Nature 218 (1968) 167–169. [38] J.C. Anderson, C. Eriksson, Piezoelectric properties of dry and wet bone, Nature 227 (1970) 491–492. [39] W.S. Williams, Sources of piezoelectricity in tendon and bone, CRC Crit. Rev. Bioeng. 2 (1974) 95–117. [40] R.S. Lakes, The role of gradient effects in the piezoelectricity of bone, IEEE Trans. Biomed. Eng. BME27 (1980) 282–283. [41] D.A. Chakkalakal, M.W. Johnson, R.A. Harper, J.L. Katz, Dielectric properties of fluid saturated bone, IEEE Trans. Biomed. Eng. BME-27 (1980) 95–100. [42] R.S. Lakes, J.L. Katz, Dielectric relaxation in cortical bone, J. Appl. Phys. 48 (1977) 808–811. [43] E. Korostoff, Stress generated potentials in bone: relationship to piezoelectricity of collagen, J. Biomech. 10 (1979) 41–44. [44] A. Bur, Measurements of the dynamic piezoelectric properties of bone as a function of temperature and humidity, J. Biomech. 1 (1976) 495–507. [45] E. Korostoff, A Linear piezoelectric model for characterizing stress generated potentials in bone, J. Biomech. 12 (1979) 335–347. [46] A. Gjelsvik, Bone remodelling and piezoelectricity II, J. Biomech. 6 (1973) 187–193. [47] N. Guzelsu, A piezoelectric model for dry bone and tissue, J. Biomech. 11 (1978) 257–267. [48] M. Johnson, D. Chakkalakal, R.A. Harper, J.L. Katz, Comparison of the electromechanical effects in wet and dry bone, J. Biomech. 13 (1980) 437–442. [49] D.A. Chakkalakal, M.W. Johnson, Electrical properties of compact bone, Clin. Orthopaed. Relat. Res. (161) (1981) 133–145. (November–December). [50] K.R. Rajagopal, L. TaoMechanics of Mixtures, Series on Advances in Mathematics for Applied Sciences, vol. 351995, World Scientific Publishing Co. Pvt. Ltd. [51] H. Demiray, N. Güzelsu, A mixture model for wet bones – I theory, Int. J. Eng. Sci. 15 (1977) 707–718. [52] H. Demiray, N. Güzelsu, Electromechanical properties and related models of bone tissues, Int. J. Eng. Sci. 17 (1979) 813–851. [53] D.A. Chakkalakal, Mechanoelectric transduction in bone, J. Mater. Res. 4 (1989) 1034. [54] S. Ramtani, M. Zidi, Damaged-bone adaptation under steady homogeneous stress, ASME J. Biomech. Eng. 124 (2002) 1–6. [55] H. Demiray, A.C. Eringen, Continuum theory of slightly ionized plasma, diamagnetic effects, Plasma Phys. 15 (1973) 903–920. [56] R.M. Bowen, Theory of mixtures, in: A.C. Eringen (Ed.), Continuum Physics, vol. 3, Academic Press, 1976, pp. 1–127. [57] R.A. Grot, Relativistic continuum physics: electromagnetic interactions, in: A.C. Eringen (Ed.), Continuum Physics, Academic Press, vol. 3, 1976, pp. 129–219. [58] T.C. Lee, E.R. Myers, W.C. Hayes, Fluorescence-aided detection of microdamage in compact bone, J. Anat. 193 (1998) 179–184. [59] S. Ment, Effects of seven days of continuous capacitive electrical stimulation on bone growth around titanium implants in the rat tibia, Master Degree Thesis, Mc Gill University, Montreal, 1999. [60] D.P. Nickerson, N.P. Smith, P.J. Hunter, A Model of Cardiac cellular electromechanics. The integrated heart: modelling cardiac structure and function, Philos. Trans.: Math., Phys. Eng. Sci. 359 (1783) (2001) 1159–1172. [61] Peter Kohl, Frederick Sachs, Mechanoelectric feedback in cardiac cells. The integrated heart: modelling cardiac structure and function, Philos. Trans.: Math., Phys. Eng. Sci. 359 (783) (2001) 1173–1185. [62] K.R. Rajagopal, A.S. Wineman, M.V. Gandhi, On boundary conditions for a certain class of problems in mixture theory, Int. J. Eng. Sci. 24 (1986) 1453– 1463. [63] L. Tao, K.R. Rajagopal, On boundary conditions in mixture theory, in: K.R. Rajagopal (Ed.), Recent Advances in Elasticity and Viscoelasticity, World Scientific Publishing, Singapore, 1995, pp. 130–149. [64] K.R. Rajagopal, On implicit constitutive theories, Appl. Math. 48 (2003) 279–319. [65] M. Massoudi, Boundary conditions in mixture theory and in CFD applications of higher order models, Comput. Math. Appl. 53 (2007) 156–167. [66] G.M. Johnson, M. Massoudi, K.R. Rajagopal, Flow of a fluid–solid mixture between flat plates, Chem. Eng. Sci. 46 (1991) 1713–1723. [67] G.M. Johnson, M. Massoudi, K.R. Rajagopal, Flow of a fluid infused with solid particles through a pipe, Int. J. Eng. Sci. 29 (1991) 649–661. [68] M. Massoudi, K.R. Rajagopal, T.X. Phuoc, On the fully developed flow of a dense particulate mixture in a pipe, Powder Tech. 104 (1999) 258–268. [69] A.C. Eringen, Electromagnetic theory of microstretch elasticity and bone modelling, Int. J. Eng. Sci. 42 (2004) 231–242.