- Email: [email protected]
A mixture theory for heat-induced alterations in hydration and mechanical properties in soft tissues
International Journal of Engineering Science 39 (2001) 1535±1556
www.elsevier.com/locate/ijengsci
A mixture theory for heat-induced alterations in hydration and mechanical properties in soft tissues L. Tao, J.D. Humphrey, K.R. Rajagopal * Mechanical Engineering and Biomedical Engineering, Texas A&M University, College Station, TX 77843-3123, USA Received 18 October 2000; accepted 20 December 2000
Abstract We propose a framework to describe the mechanical behavior of thermally damaged biological soft tissues based on a melding of ideas from classical continuum thermodynamics and irreversible thermodynamics. Relevant assumptions are delineated, details of the constitutive formulation are provided, and a set of interfacial conditions are derived that can be used as the boundary conditions for certain diusion processes. The utility of the framework is illustrated by considering two problems that correspond to tractable experiments: transverse diusion through a ®nitely deformed membrane, and radial diusion in a thermally damaged cylindrical tissue. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Hydration; Mixture theory; Soft tissues; Thermal damage
1. Introduction Advances in laser, microwave, radio-frequency, and similar technologies continue to encourage the search for new applications of heat (i.e., thermal energy) to treat disease and injury. Examples of heat-treatment can be found in cardiology, dermatology, gynecology, oncology, ophthalmology, orthopedics, and urology, to name a few areas. The eects of increased temperature levels on cells and tissue are manifold, including the increased expression of heat-shock proteins, outright cell death, the denaturation of proteins, and localized edema. Associated with these changes, of course, are possible changes in the electrical, optical, mechanical, and thermal properties of the treated tissue. To date, however, there has been little attempt to quantify rigorously such changes in properties even though this information is essential for the analysis
*
Corresponding author. Tel.: +1-979-845-1251; fax: +1-979-845-3081. E-mail address: [email protected] (K.R. Rajagopal).
0020-7225/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 1 9 - 2
1536
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
and design of clinical heating strategies. Rather than address a speci®c clinical application or a particular data set from the literature, the purpose of this paper is to outline a possible approach for mathematically modeling the changes in hydration and non-contractile mechanical properties that arise due to heat-treatment of a soft collagenous tissue. Towards this end, we shall model a typical soft tissue as a ¯uid-infused solid that experiences ®nite deformations and exhibits nonlinear behavior as a continuum mixture. Moreover, we assume that severe heattreatment results in the coexistence of varying degrees of native and denatured collagen. Hence, a theoretical framework for modeling such a system must account for the behavior of a material consisting of at least three constituents, one of which (mobile ¯uid) is capable of motion with respect to the other two, which in turn behave as a constrained mixture (native and damaged tissue). Most heat-treated tissues and organs consist of parenchymal tissue (collagen, muscle, etc.) surrounded by delimiting membranes; examples range from tendons to muscular organs such as the heart and uterus. As such, these materials are necessarily non-homogeneous and the mechanical properties vary with the material point. Thermally-induced changes likewise vary from one material point to another. Moreover, given that the delimiting membranes serve as semipermeable barriers, heat-induced alterations in permeability must be accounted for as well. Although a complex problem, we can draw upon the many advances in the continuum theory of mixtures (e.g., see [3,4,20,24]) to formulate a general ®rst approach. We note, however, that certain modi®cations will be needed, and hence we will meld some ideas from classical continuum (e.g., [25]) and irreversible (e.g., [8,9,11,12]) thermodynamics. Our overall philosophy is that initial (empirical) observations are needed to delineate the general characteristics of the problem that must be modeled. Once done, one can then establish an appropriate (albeit perhaps simpli®ed) theoretical framework to guide the design and interpretation of speci®c experiments, which in turn yield constitutive relations for evaluation. For purposes herein, we shall rely primarily on observations from heat-treated tendons. Reasons for this are twofold: ®rst, much of the relevant data in the literature has come from rat tail tendon or chordae tendineae; second, being thin cylindrical structures consisting primarily of collagen, these tissues admit tractable initial-boundary value problems for study both analytically and in the laboratory. 2. Background Collagen is the most abundant protein in the body; it exists in at least 15 dierent types, each similar but individually specialized for a speci®c function. Herein, we will be concerned primarily with the ®brillar type I collagen and the meshed type IV collagen. The former is the primary constituent in the ``core'' of a tendon, endowing the tissue with marked tensile strength, and the latter is the primary constituent of the surrounding delimiting membrane or ``sheath''. Collagen molecules are composed of three polypeptide chains, each having some 1300±1700 amino acid residues, the majority of which (1000±1400) are organized into a central triple helical structure [2]. The latter is on the order of 285 nm long and 1.5 nm in diameter and consists of a repeating triplet of amino acids, (G±X±Y)n , where G is glycine and X and Y are often proline or hydroxyproline. The collagen molecule is stabilized by extensive hydrogen
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1537
bonds, with associated water playing a key role. In ®brillar types, the collagen molecules form a quarter-staggered arrangement (which gives rise to the characteristic 67 nm bandedness), with lysine and hydroxylysine residues forming stable intra- and inter-molecular cross-links particularly in the non-helical end domains. In contrast, type IV collagen is stabilized via extensive disulphide bonding. As a result, the activation energy of type IV collagen is about twice that of type I [15]. The eects of heat on collagen can be ``reversible'' or ``irreversible''. For example, moderate heating can result in a local unfolding within the protein that is reversed upon the restoration of normal temperatures. Severe heating results in a time-dependent irreversible transformation of the native triple helical structure into a more random (coiled) structure [15,17]. It is thought that the latter transformation occurs primarily via the breaking of long sequences of hydrogen bonds that stabilize the triple helix; particular sub-domains along the molecule (e.g., between residues 877 and 936; [15]) may be more susceptible to the breaking of such sequences. Heat-induced breakage of reducible cross-links [14] may also play a role. Regardless of the precise mechanisms, on a gross scale denaturation results in an irreversible shrinkage of the tissue, which for us becomes a convenient metric of thermal damage. This shrinkage may result due to exposure of hydrophobic regions of the molecule that are initially within the center of the helical structure. Regardless, Hormann and Schlebusch [10] showed that gross shrinkage provides much the same information on denaturation kinetics as many biochemical assays (though not the same detail on the underlying mechanisms of course). Denaturation also results in a gross change in the hydration of the collagen, which is thought to be due ®rst to the liberation of water and then a subsequent increased absorption of water via water bridges. Again, independent of precise molecular mechanisms, on a gross scale heating and post-heating recovery involves changes in hydration. Much work has been done on the diusion of a ¯uid through an elastic solid ([5,13,20,22], to mention just a few). These studies have been mainly concerned with the application of mixture theory within the context of speci®c models, however, and none consider the unique features of bio-thermomechanical problems. In this paper we revisit the diusion problem, speci®cally the diusion of water through a thermally damaged biological soft tissue that experiences ®nite deformations. We shall treat the body as a mixture of mobile water with weak hydrogen bonding and a porous solid with tightly bonded water molecules. Further, upon damage, some hydrogen bonds are broken during the heat treatment, and the associated region of the native tissue ``transforms'' into an altered material. The theory of interacting continua [1,4,20,25] provides a natural setting for the study of such problems. As noted above, we will infuse certain ideas from irreversible thermodynamics into the classical theory of interacting continua. We feel that such a thermodynamic setting provides a reasonable basis for studying the problem at hand. While it is somewhat formal in its approach, in view of the lack of a clear understanding of the basic mechanisms (physical and physiological) that govern the problem, such a methodology can lend insight into the various issues that have to be considered further in the theory. We will state explicitly the assumptions that are introduced in the analysis, so that we can recognize the status and relevance of these assumptions, knowledge of which may lead to possible modi®cations and improvements in further studies. A feature of the present analysis is the formulation of a set of interfacial jump conditions that can be adopted as boundary conditions in describing some speci®c diusion processes in tractable experiments.
1538
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
3. Preliminaries The body consists of mobile water, loosely bonded with the tissue, and a porous solid having tightly bonded water molecules. The tissue may be partitioned further into two sub-components: one associated with that part wherein the hydrogen bonds are broken due to thermal damage and another that is the intact native tissue. It is expected that the damaged and the native constituents intermingle (over representative volume elements), with inter-conversion occurring from damaging and healing. We shall suppose that the solid constituents undergo at a point the same motions in a homogenized sense. Within the context of mixture theory, a set of thermodynamic quantities can be assigned to each component: quantities associated with the mobile water and the tissue will be denoted by the superscripts f and t, respectively, and quantities associated with the damaged and native subcomponents will be denoted, respectively, by the superscripts d and n. Thus, the following equations describe the balance of mass, linear momentum, and energy for each component a, oqa div qa va ma ; ot X
ma 0;
a f ; n; d;
1
mf 0;
2
n;d
o a a T q v div qa va va div Ta pa ; ot
a f ; t;
3
qt qn qd ; X
4
pa 0;
5
f ;t
o a a 1 a2 1 a2 a a a q jv j div q jv j v div Ta va ot 2 2
qa pa va sa ;
a f ; t; 6
X
pa va sa 0:
7
f ;t
Here, qa , va , Ta , a and qa are, respectively, the mass density, velocity, Cauchy stress tensor, internal energy density, and heat ¯ux associated with the ath component; ma , pa and sa denote, respectively, the mass density conversion rate, the internal momentum supply, and the energy supply to the ath component resulting from the interactions between the components.
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1539
It should be emphasized that no electro-magnetic eect is included explicitly here, though hydrogen bonding can occur as follows. The external electro-magnetic ®eld is assumed negligible, the component f is neutral, and the component t is also neutral since the charge of the porous solid and the charge of the bonded water molecules supposedly neutralize each other due to the bonding. If the charge per unit mass of the solid is treated as constant in time, the charge distribution c at time t throughout the solid is determined completely by the charge distribution c0 in any reference con®guration j and the deformation of the solid. Therefore, the charge c is not needed explicitly in the formulation provided that the eects of c0 on ma , Ta , pa and so on are properly accounted for. Next, we infuse into our framework a thermodynamic notion that is widely used in nonequilibrium thermodynamics. The entropy density ga for component a can be introduced, for example, under the hypothesis of local equilibrium, and a balance equation can be proposed that has an analogous structure to the balance equations of (1), (3) and (6) [8,27], o a a q g div qa ga va div Ja ra ; ot
a f ; t;
8
where Ja and ra are, respectively, the entropy ¯ux and the internal entropy production rate per unit volume. The separation of the right-hand side of (8) into two parts leads to a natural representation due to the restrictions imposed by the second law of thermodynamics for isolated systems, as demonstrated below. Several observations are noteworthy here. First, by adopting the above equations, we tacitly restrict ourselves to motions without possible rapid and random ¯uctuations about some sensible mean value; as such, no equations need be posited for ¯uctuations. Examples of such ¯uctuations are the presence of turbulence or discontinuous distributions of a component a, which can be severe if the characteristic dimension of the internal geometric structure of the physical components, like pore size, is large. Neglecting ¯uctuations in the present formulation is justi®ed because the diusion process is slow (on the order of hours), and the water and tissue components are distributed regularly within an internal geometric structure that has a small characteristic dimension l. Second, the above equations are local. Introduction of continuously dierentiable ®elds qa , va , a J , ga , ra and so on is based on the assumption that each component can be modeled as a continuum and the hypothesis of local equilibrium, that is, these ®elds are averaged over some characteristic volume V or area, large compared with the characteristic dimension l. In formulating governing equations, independent of primitive quantities that are self-evident but unde®ned, we nonetheless have to motivate some inferential method for quantifying them. For instance, whereas temperature cannot be measured directly, we indirectly infer it by recognizing that the change in the length of a mercury column bears a one to one relationship with changes in temperature. The diculty with postulating new governing equations such as (8) is that, at this moment, we have no means to motivate a constitutive structure for quantities such as Ja or an inferential basis for measuring them. Nonetheless, as we shall see, such formalism does lead to ways for developing certain meaningful restrictions on the class of processes that we are interested in. With this understanding, we proceed to derive a restriction on ra by applying (8) to an isolated mixture U with U V,
1540
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Z
o r dv ot U
Z X U
qa ga dv P 0;
f ;t
r :
X
ra ;
9
f ;t
where the inequality required by classical thermodynamics guarantees a non-decreasing entropy in an isolated body. Next, taking into account the relative arbitrariness of U ( V), treating V as small enough due to the smallness of l, and assuming that the form of r is independent of whether the mixture in motion is isolated or not, we get X ra P 0; 10 r f ;t
for the motion of the mixture in general. Third, the hypothesis of local equilibrium may allow us to apply the concept of ga and (8) to a motion of the mixture away from a speci®c equilibrium. The reasoning is that in each local volume and during a short interval of time, the local motion is not far away from a corresponding equilibrium state, and ga and (8) can be de®ned according to the hypothesis of local equilibrium. If the motion is far away from any special state of equilibrium, then, unlike in irreversible thermodynamics, we do not necessarily need to introduce linear forms for the constitutive formulation based on any single equilibrium state, as is usual in rational and extended thermodynamics. Fourth, ga as well as a will be treated as a derived quantity completely determined by quantities such as mass density, velocity, deformation gradient, temperature and so on. That is, ga as well as a are treated as a state variable as in classical thermodynamics. Consequently, Eqs. (8) and (10) can be linked directly to the balance equations (1)±(7), and some restrictions can be imposed on the forms of ma , qa , pa , Ta , and Ja following the standard practice in the thermodynamics of continua [8,16,25]. 1 We are aware of the uncertainties associated with such a scheme, or for that matter any other scheme adopting an entropy balance equation to obtain the structure of undetermined quantities contained in the other balance equations. For instance, there is ¯exibility, to some extent, in proposing the forms of the entropy ¯ux and the (non-negative) internal entropy source, which is due partly to our lack of knowledge on how to quantify and/or measure entropy in a complex material undergoing a complicated motion. This uncertainty thus requires that the results from this scheme be tested experimentally, if possible. Finally, consider some issues related to the solid component. We have assumed that the two sub-components of the tissue have individual mass densities but undergo the same motions, i.e., qt qn qd ;
vn vd vt :
11
Then, from Eqs. (1) and (2)1 , we get oqt div qt vt 0: ot 1
12
When the ¯uctuations are not negligible, (10) may not play the same role as in rational thermodynamics or irreversible thermodynamics on the closure of the above equations since an additional assumption is needed, as in the ¯ows of bubbly liquids [23].
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1541
If we introduce a mass fraction for the damaged solid, m
qd ; qt
13
Eqs. (1) and (12) yield o md t v grad m t : ot q
14
Eq. (12) can be integrated to yield qt det F qt0 :
15
Here, qt0 is the eective mass density of the tissue component in the original reference con®guration j0 of the tissue, and F is the deformation gradient of the solid with respect to j0 . This j0 might be identi®ed with the state of the tissue in which qd 0 or qd is negligible compared with qn and it is (almost) stress free. As long as there is no damage, the response of the material is elastic. As the thermal damage takes place, qd increases, and the response is no longer purely elastic. The sub-component of damaged constituent will have a dierent reference con®guration jd from which its response is elastic. Motivated by the idea of multiple natural con®gurations [18], we introduce a reference con®guration j that is intermediate between j0 and jd and may evolve due to thermal damage; we can then decompose the deformation gradient F into a part Fj that is associated with an elastic response and another that is associated with the evolution of natural con®gurations. If the con®guration j is determined by a quantity G, a mapping from j0 to j relates F and Fj through Fj FG 1 ;
16
plus an equation for the evolution of G. 4. Speci®c model To apply the above framework to study the diusion of water through a thermally damaged tissue, we de®ne here q :
X
qa ;
f ;t
1X a a q; : q f ;t La : grad va ;
1X a a 1X a a q v ; g : qg; q f ;t q f ;t X q : qa ; ua : va v;
v :
f ;t
d o : v grad: dt ot
Then Eqs. (1)±(7), (8) and (14) yield
17
1542
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
dqa dt
ua grad qa
dm dt
ut grad m
qa tr La ma ;
md ; qt
q
d X tr Ta La dt f ;t
q
X dg Ja div dt f ;t
18
19
pa va
div qa qa a ua ;
qa ga ua r:
20
21
We now assume that ``both'' components ± mobile water and tissue ± have the same absolute temperature ®eld T, and the entropy and the internal energy of the mixture are completely determined by the primary ®elds qf , qd , T, Fj , and G and the reference con®guration j fg; g fg; g qf ; m; T ; Fj ; G; j
22
which means that the entropy and the internal energy of the mixture are modeled similar to those for an elastic process with deformation Fj with respect to j for ®xed qf and m. The dissipative response due to the viscosity of water will be accounted for through the internal momentum supply pf to be introduced and discussed later. We should mention that this treatment is intended to apply only to the diusion process considered here, and other quantities might have to be included if dierent processes of the mixture are to be considered. For example, g and may depend on the velocity gradient Lf if visco-elasticity is signi®cant in the process under consideration. We restrict ourselves to a special class of motions for which the assumption (22) is reasonable. Introducing the Helmholtz potentials wa : a
1X a a qw; q f ;t
23
T g w qf ; m; T ; Fj ; G; j :
24
T ga ;
w :
we have w
Substituting (24) into (21) and using (18) and (20) results in
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
"
1543
!
X 1 a ow dT ow a a a a T r T div J q q w u q g tr Tt qFj T Lt T oT dt oFj f ;t # ow f f f ow f f f f ow T qq I L p q gradw q f gradq q gradT vt f oq oq oT X 1 ow md ow T T ow T qa qa a ua grad tr q F G q j T om qt oFj oG f ;t o : vt grad GT P0: ot
vf
25
Moreover, we will assume that the tissue and the mobile water are intrinsically mechanically incompressible (though thermally compressible) and that the volume additivity constraint holds, i.e., qf qt 1; qfR qtR
26
where qfR and qtR ( qt0 = 1 b0 ) are the physical mass densities of the water and tissue in j0 , respectively, and b0 is the volume fraction of mobile water in j0 . This relation can be partially justi®ed by the fact that the tissue component includes plentiful water molecules bonded tightly with the solid. With the presence of these bonded water molecules, the tissue component and the mobile water component may avoid `void ®lling' at the molecular level, and consequently volume additivity holds. Next, in view of the range of temperature changes over the processes modeled here, it is appropriate to treat qfR and qtR as constant. It then follows from (26), (1), and (15) that t q grad tR vt q
f
v
qt qf tr tR Lt fR Lf q q
0:
27
Since Eq. (25) is expected to hold under the above constraint, a Lagrange multiplier p is introduced to enforce the constraint thus leading to (" X 1 a ow dT a a t a a fR ow J q q w u q g tr T p q q f T r T div T oT dt oq f ;t # ) f qt ow qf q tR I qFj T Lt Tf p fR I Lf pf pgrad fR qf gradw q q q oFj X 1 ow md f ow t f a a a a gradT v v T q q q u grad q oT T om qt f ;t ow o T T ow t 28 tr q Fj G v grad GT P0: oFj oG ot
1544
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Here we have used p qfR q ow=oqf instead of p in the mathematical treatment for the reason to be clear below. With the inclusion of qf in w, the Helmholtz potential should be independent of det Fj G det F on the basis of (15) and (26), i.e., w w qf ; m; T ; Fj ; G; j with det Fj G absent. Eq. (28) can be simpli®ed if we choose Ja
1 a q qa wa ua T
29
and g
ow ; oT
30
which are the usual assumptions that are made. This eliminates the divergence term and the time derivative of the temperature term, thus Eq. (28) yields, (" # ) t q ow qf t t f fR ow I qFj T L T p fR I Lf T r tr T p q q f oq qtR q oFj f X q ow pf p grad fR qf grad w qf qa qa a ua grad T vt vf T oT q f ;t d 1 ow m ow ow o grad tr q FTj G T 31 q vt grad GT P 0: T om qt oFj oG ot Inequality (31) places some restrictions on the structures of the functions pf , Tt , Tf , q, md and G. Notice that there is no unique way to arrange the cross terms involved in the velocity dierence and the gradient of temperature; this may be immaterial if the cross eects are handled properly through the appropriate choice of constitutive relations for pf and q. Another issue is that (29) and (30) are dictated here by the elimination of the divergence and the time rate of change of temperature and they should be considered as assumptions. 2 Relation (29) will be used later in generating the interface conditions in the isothermal process, which may provide a way to check its validity indirectly through experiments. Now we propose a set of constitutive relations, which reduce to that proposed by Shi et. al. [22] in the case of isothermal processes without damage, t q ow t fR ow p q q f I qFj T ; 32 T tR oq q oFj
2 It is possible, for example, to add a term like Ea ua to the right-hand side of (29) for Ja where Ea is a quantity to be determined. This speci®c addition will aect the forms of pf , Tt and Tf , but will not aect the balance of linear momentum equations in the case of symmetric Ea for isothermal processes, a result similar to that of Shi et. al. [22] which was proposed from a dierent perspective.
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Tf
p
qf I; qfR
f q p p grad fR q f
X f ;t
md
33
qf grad w qf
ow grad T b vt oT
vf ;
1 q q u C grad ; T a
cqt
1545
34
a a a
35
ow ; om
36
o ow T t G v grad G b FTj ot oFj
T
ow : oG
37
Here, c P 0, b P 0, b and C are, respectively, the drag coecient and heat conductivity, general positive de®nite second-order tensors required by the inequality (31) and to be determined by experiments. The reason not to treat b simply as a multiple of the unit tensor is due to the possible anisotropy of the pore structure of the tissue. Physically, b characterizes the momentum exchange (due to viscosity) between the viscous water, which moves relative to the solid constituents, and the tissue; it is important to retain this term even though the diusion is slow since the `channels' for the diusion are very small and the contact area per unit volume between the water and tissue is very large. This term should play such a role more eectively than a relation for Tf which includes the viscous eect as an internal source term in the equations of motion, since the standard form of Tf , in terms of the velocity gradient in the ¯uid, may not describe the eects due to interactions in the pores adequately. A simple example to illustrate this inadequacy is to consider a diusion velocity ®eld of the form (vf1 , vf2 , vf3 vf x1 ; 0; 0. This velocity ®eld would not yield f if the standard velocity gradient form were employed, yet it would give rise any viscous term in T12 to a viscous drag locally. Next, associated with (36) and (37), an activation criterion like a yield condition in plasticity could be introduced if necessary that is based on temperature. Finally, it is worth mentioning that we can also obtain the above relations on the basis of maximizing the entropy production r, following Rajagopal and Srinivasa [19]. The form of w in the above relations can be restricted further by invariance requirements. Here, we require Galilean invariance [21,25], instead of the more strigent and commonly used assumption of material frame indierence (MFI). (The issue regarding MFI will be discussed in detail elsewhere.) This invariance requires that the form of w be invariant under the coordinate transformation x Qx ct a; where Q and c are constant with QT Q QQT I. Therefore,
1546
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
w qf ; m; T ; Fj ; G; j w qf ; m ; T ; Fj ; G ; j w qf ; m; T ; QRj Uj ; G; j;
38
with the help of qf qf , m m, T T , G G, Fj QFj , Fj Rj Uj , where Rj is proper orthogonal and Uj is symmetric and positive de®nite. Now pick Q RTj t0 for some instant t0 and (38) yields w qf ; m; T ; Fj ; G; j w qf ; m; T ; RTj t0 Rj Uj ; G; j
8t;
and thus w qf ; m; T ; Fj ; G; j w qf ; m; T ; Uj ; G; j at t0 : The arbitrariness of t0 results in w depending on j, qf , m, T, Cj U2j FTj Fj and G but not through det Cj GGT det Fj G2 ). Additional restrictions on the form of w can be obtained if the symmetry of j is taken into account. As a proper approximation for the processes modeled here, we may assume that the previously introduced reference con®gurations j0 and jd have the same symmetry (G0 Gd ) or that jd has a dierent symmetry (but G0 Gd ). For example, if the tissue is assumed to be orthotropic with symmetry planes X1 constant and X2 constant in j, w may be taken as [7], 2
2
2
w qf ; m; T ; Cj ; G; j w qf ; m; T ; Cj11 ; Cj22 ; Cj33 ; Cj13 ; Cj32 ; Cj21 . . . ; j:
39
Alternatively, if the tissue is assumed to be transversely isotropic with its preferred direction along the X3 direction in j0 , w will take the following form, 2
2
w qf ; m; T ; Cj ; G; j w qf ; m; T ; I; II; Cj33 ; Cj13 Cj32 . . . ; j;
40
where I : tr Cj ;
II :
1h 2 tr Cj 2
i 2 tr Cj :
41
The ®nal, speci®c form of w needs to be ®xed through experiments. To have a better understanding of the role of w in the process treated here, consider the case of isothermal slow diusion, which occurs during the ``recovery'' process post-heating [6]. With this premise, the inertial term in the balance of linear momentum for the water component f is likely negligible, and Eq. (3) (with a f ) reduces to vf
vt
qf 1 b grad p qfR w : fR q
42
Eqs. (33) and (34), with grad T neglected, have been used in the derivation. This relation indicates that there are two parts aecting the isothermal slow diusion of water: the ®rst part is associated with the gradient of the spherical stress of the water component (p), similar to the commonly used Darcy's law; the second part is associated with the gradient of w, similar to Fick's law, which accounts for variations in the eective mass density qf (concentration) and the
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1547
deformation F. Also, the inverse of the drag coecient b 1 is part of the diusion coecient, which further shows the necessity for treating b as a general positive de®nite second-order tensor in the case of anisotropic tissues. By identifying that the spatial variations of both the spherical part of the stress of the water component and the Helmholtz potential are the driving forces for the diusion, we need to develop proper boundary conditions, an issue to be discussed in the next section. 5. Interfacial conditions In this section, we discuss interfacial jump conditions that are necessary to study the diusion process. We will deal with three classes of interfaces: (1) mixtures present on both sides; (2) a mixture on one side and pure water on the other side; and (3) a mixture on one side and water within a rigid porous solid on the other side. Consider an interface S in a mixture. On each side of S, denoted by the superscripts and respectively, there may be mobile water (denoted by the superscript f) and tissue (denoted by the superscript t). The two tissue components may dier from each other or one component may be absent. The interface links the two tissue components perfectly if both are present, that is, x Xt ; t x x x x Xt ; t always holds, if it held initially, where Xt denotes the location of particles from both tissues. Further, we assume that there are no physical properties attached to the interface except those derivable from the components on both sides of S while x ! x , which implies that there is no accumulation of mass, linear momentum, energy or entropy on the interface. Nonetheless, this interface is a material surface since it is composed of the same particles either from the contact of two types of tissues or from one side of the interface being occupied by a tissue and the other side free. Consider a ®eld quantity /a , associated with component a on either side of S, which can be characterized by the following balance equation: o a a q / div qa /a va div pa ga ; ot
43
where pa and ga are, respectively, the ¯ux and the source related to /a . Choosing a control volume around S, applying the balance of /a for the mixture to the control volume, and shrinking the control volume, we obtain [20,26] f f f q / v
vt
pf pt
n qf /f vf
vt
pf pt
n;
44
where n is the normal to the interface. This relation provides the required interfacial jump conditions by assigning proper physical quantities for /a and pa . Class 1. For an interface with tissues present on both sides, picking, respectively, (i) /f 1, T T f p pt 0 (balance of mass) and (ii) /f vf , pf Tf , pt Tt (balance of linear momentum), we obtain from (44) f f q v
vt
n qf vf
vt
n;
vt vt ;
45
1548
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
" f f
q v v
f
t
v n
X
#
"
aT
f f
T n
q v v
f
t
X
v n
f ;t
# aT
T n
:
46
f ;t 2
Further, picking, respectively, (iii) /f f jvf j =2, pf Tf vf qf , pt Tt vt qt (balance of energy) and (iv) /f gf , pf Jf qf qf wf uf =T , pt Jt qt qt wt ut =T (entropy equation), and combining the two resultant jump conditions to eliminate the heat ¯ux qf qt , we obtain " # X 1 2 qa qa a ua Ta va n vf vt qf jvf j 2 f ;t " # X 1 2 qf jvf j qa qa a ua Ta va n vf vt 2 f ;t
47
and " 1 f 2 f vf q w jv j 2
v
t
X
# a a
Tv
f ;t
" 1 f 2 f n q w jv j vf 2
v
t
X
# a a
Tv
n:
48
f ;t
Class 2. Similarly we can obtain the jump conditions across the interface, on the side ()) where the mixture is present and on the other side () occupied by the pure water, taking into account that Tf P I where P is the water pressure (the viscous eect in the pure water region is ignored since the ¯ow of water is slow and uniform and thus the velocity gradient is very small),
qf vf
" qf vf vf
vt
n qfR vfR
t
v n
X
vt #
TaT n
n;
vt vt ;
qfR vfR vfR
49
vt n P n ;
50
f ;t
" # X 1 qf jvf j2 vf vt qa qa a ua Ta va n 2 f ;t 1 fR 2 fR fR fR t fR fR q jv j n: v q Pv v 2 " 1 f 2 f vf q w jv j 2
v
t
X f ;t
# a a
Tv
fR
n q
1 w jvfR j2 2 fR
51 fR
v
t
v Pv
fR
n:
52
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1549
Here fR and wfR are, respectively, the internal energy and the Helmholtz potential for pure water. Class 3. In the case that on the side ()) of the interface the mixture is present and on the side there is water and a static rigid porous solid through which the water ¯ows, we need to modify the above relations. We will assume that the rigid solid has such a large porosity a0 that it has negligible eect on the ¯ow of the water compared to the tissue, but the pores are small enough that the geometry of the interface is completely determined by the rigid solid. Also, we assume that the rigid solid will exert little eect on the tissue along the directions tangent to the interface. We can take qf a0 qfR , f fR , wf wfR , Tf a0 P I, qf a0 qfR and qt qs where qs is the heat ¯ux associated with the rigid porous solid, and obtain, by modifying (49)±(52) f f n a0 qfR vfR n; q v
53
" # 1 f 2 f X a a a f a a a q q u T v n q jv j v 2 f ;t 1 2 a0 qfR fR jvfR j vfR a0 qfR qs a0 P vfR n; 2 1 f 2 f f q w jv j v 2
f f
T v
54
1 fR 2 fR fR fR fR n a0 q w jv j v a0 P v n: 2
55
The relation associated with the balance of linear momentum is not useful here since the stress inside the rigid porous solid is indeterminate. We should mention here that, in the case of isothermal processes, the interfacial conditions of (47), (51) and (54) associated with the balance of energy will be left out. 6. Applications 6.1. Diusion through a membrane Consider the slow isothermal diusion of water through a thin membrane of original thickness H without thermal damage. Suppose that the membrane is orthotropic, with X1 0 and X2 0 as the symmetry planes, and 2
2
2
w w qf ; T ; C; j0 A qf ; T ; C11 ; C22 ; C33 ; C13 ; C32 ; C21 ; j0 :
56
While laid on a rigid porous plate at X3 0, the membrane undergoes the deformation x1 kX1 ;
x2 kX2 ;
x3 f X3 ; t;
where k is a constant. The water moves according to
57
1550
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
vf 0; 0; vf X3 ; t ;
58
under the pressure dierence PH P0 P 0, where PH and P0 are, respectively, the pressure of water immediately above and below the membrane. A simple calculation yields 2 2 3 2 3 k 0 0 0 0 k of : 59 F 4 0 k 0 5; C 4 0 k2 0 5; g : oX3 0 0 g 0 0 g2 Balance of mass and the volume additivity constraint give 1 qt qt0 k2 g ;
60
:
61
f
q q
fR
1
qt0 2 kg qtR
1
Balance of linear momentum for the mixture as a whole and for the water component yields, o f t T33 T33 0; ox3
62
f 1 o q p qfR w b33 fR vf ox3 q
of ot
0:
63
Here we have neglected inertia and assumed that Tt , Tf and pf are described by (32)±(34) and (56) with b13 b23 0. The boundary conditions for these two equations can be derived from (49) to (55) and the ¯ow conditions mentioned above: f jX3 0 0;
p qfR w p qfR w f t T33 T33
64 P0 qfR wfR ;
65
X3 H
PH qfR wfR ;
66
X3 H
67
X3 0
PH :
Solving the above equations, we get qt ow ow 2 q w q tR f 2g fR wfR q oC33 X3 0 q oq
PH P0 ; qfR
68
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Z
H
0
b33
qf qfR
1
f
v
of ot
1551
PH
g dX3
P0 :
69
Moreover, taking into account the thinness of the membrane, we can obtain the mass ¯ux from (69) F0
q v
"
f f
X3 0
qf 2 gHb33
# X3 0
PH P0 ; qfR
70
with g to be determined from (68) and qt and qf from (60) and (61). We may interpret gHb33 in (70) as the drag coecient for the membrane as a surface. 6.2. Radial diusion in a membrane-covered cylindrical tissue Consider a cylindrical tissue with a membranous covering, with the core tissue occupying the region 0; RI and the membrane occupying the region RI ; RO in the reference con®guration j0 . The core is assumed to be a mixture of a transversely isotropic tissue and water with its preferred direction along the axial direction Z in j0 . Suppose that the membrane is orthotropic with symmetry planes H and Z. Initially the tissue is at temperature T0 . Then, suppose an axial load W is applied to the tissue which is then immersed into a water bath at temperature Th (> T0 ) and pressure P0 . Due to the overall thinness (RO small), let the whole tissue reach the equilibrium temperature Th rather quickly. This heating causes thermal damage to the core tissue, but the damage to the membrane is relatively small and may be neglected (this restriction can be relaxed easily in our formulation). Also, because of this heating, there is water diusion through the membrane and into or out of the tissue (both during and following heating). Thus, we have w w qf ; m; T ; Cj ; G; j
71
for the core tissue and w w qf ; T ; CRR ; CHH ; CZZ ; CRZ 2 ; CZH 2 ; CHR 2 ; j0
72
for the membrane. We shall assume that the temperature ®eld is approximated by T
T0 ; Th ;
t < 0 and t > th ; t 2 0; th ;
73
that the tissue and the membrane undergo the deformation r f R; t;
h H;
z k tZ;
that water diuses according to
R 2 0; RI [ RI ; RO ;
74
1552
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
vfr vf R; t;
vfh 0;
_ vfz kZ;
R 2 0; RI [ RI ; RO ;
75
and that the core tissue is damaged due to heating, thus giving rise to 2 3 g R; t 0 0 G4 0 h R; t 0 5 0 0 k R; t
76
in a cylindrical coordinate system. It follows that for R 2 0; RI 2 of
0
oR
F4 0 0
f R
0
2
3
0 0 5; k
1 of g oR
6 Fj 4 0 0
2
3 0 7 0 5;
0 f hR
6 Cj 6 4
k k
0
1 of g oR
0 0
2
0
f 2 hR 0
0
3
7 7 0 5; k 2
77
k
and for R 2 RI ; RO 2 of
0
oR
F4 0 0
f R
0
3 0 0 5; k
2 6 C4
of 2 oR 0 0
0 f 2
R
0
0
3
7 0 5: k2
78
Next, a straight-forward but lengthy calculation yields 1 f of q q k ; R oR t
t0
f
fR
" q q
1
1# qt0 f of ; k R oR qtR
q qf qt ;
vf
qfR qf
p p^ R; t
79
80
81 ! qt of 1 k_ f ; qtR ot 2 k 1 fR 2 q kkZ ; 2
in 0; RI [ RI ; RO , and
82
83
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
om ot
c
1553
ow ; om
84
o p^ qfR w oR
q
fR
of ovf oR ot
ovf br of oR qf oR
f
v
of ot
;
85
f 1 t of o2 f of ovf ovf f of q of of fR q f q q tR v br v q oR ot oR qfR oR ot2 ot oR ot " " # 2 2 t o q ow 1 of ow 2q of 1 of ow fR fR q w q q tR f 2q oR g oR oCjRR f oR g oR oCjRR q oq t
f hR
2
ow oCjHH
# 0;
86
" 2 og 2 of ow b 3 ot g oR oCjRR
# ow ; oGjRR
87
" 2 oh 2 f ow b 3 ot h R oCjHH
# ow ; oGjHH
88
2 ok 2k ow b 3 ot k oCjZZ
ow ; oGjZZ
89
in 0; RI , along with the axial, boundary and interfacial conditions # 2 ow qt k ow of p^ q q f f 2q dR oq qtR k oCjZZ oR 0 of W 2 ow 2qk RO RI ; f oCZZ oR R 2p
Z
RI
"
fR
p^
ow q q f oq fR
qt qtR 90
O
f jR0 0;
91
f jRR f jRR f jRR ; I
I
O
92
1554
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1 f 2 f of v p^ q w v 2 ot RRI fR 1 fR 2 q of fR fR fR of P0 q w v v br v f 2 ot qf oR R f O f of ov of ov qfR RO RI ; vf oR ot ot oR R fR
of ot
93
O
"
2 # t q ow 1 of ow 1 of of qfR q tR f 2q qfR w vf 2 vf qf vf vf g oR oCjRR 2 ot ot q oq RI ( " # 2 2 1 fR 2 2q of of ow f ow v qfR wfR 2 f oR oR oCRR R oCHH R O ) t 2 qfR of of of ovf of ovf f fR q f t of o f br v f RO RI ; q tR v q ot ot oR oR ot2 q oR q oR ot
94
RO
"
2 qt ow of ow 1 fR q w vf 2 q q tR f 2q oR oCRR 2 q oq 1 fR 2 fR fR q w ; v 2 R
of v ot
fR
f
f f
q v
v
f
of ot
# RO
95
O
p^ q
vfR
fR
1 w vf 2 2
of v ot
f
RO
1 k_ : fj 2 k RRO
P0 q
fR
1 w vfR 2 2 fR
of v ot
fR
R O
;
96
97
tR fR We should mention in fR that the above derivation we have used the approximations q q 2 =2 . Also, an approximate treatment has been adopted for the motion and jp^ R; tj q kkZ of the membrane due to its thinness as in the previous case.The appropriate initial conditions are
f R; 0 f0 R; k 0 k0 ;
o f R; 0 0; ot
_ k 0 0;
98 99
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
m R; 0 0; g R; 0 1;
1555
100 h R; 0 1;
k R; 0 1:
101
To simulate the diusion process, we need to know the forms of w and br as well as the initial con®guration. This must await appropriate experiments, which may now be designed given this framework. 7. Summary We have developed a mathematical model for the diusion of water through a biological soft tissue that has suered thermal damage, based on a melding of ideas from classical continuum thermodynamics and irreversible thermodynamics. A set of interfacial conditions were derived that furnish boundary conditions that may be relevant to several classes of problems. Two speci®c processes of diusion, one for transverse diusion through a membrane and one for radial diffusion in a cylindrical tissue covered by a membrane, were considered because the predictions of the theory can be tested by experiments. Acknowledgements This work was supported in part by NSF grant (BES-9896166), NIH grant (HL-56135), and ARP grant (036327). References [1] R.J. Atkin, R.E. Craine, Continuum theories of mixtures: basic theory and historical development, Quart. J. Mech. Appl. Math. 29 (1976) 209±244. [2] S. Ayad, R.P. Boot-Handford, M.J. Humphries, K.E. Kadler, C.A. Shuttleworth, The Extracellular Matrix Facts Book, Academic Press, New York, 1994. [3] A. Bedford, D.S. Drumheller, Theories of immiscible and structured mixtures, Int. J. Eng. Sci. 21 (1983) 863±960. [4] R.M. Bowen, Theory of Mixture, in: A.C. Eringen (Ed.), Continuum Mechanics, Academic Press, New York, 1976. [5] R.M. Bowen, Incompressible porous media models by use of the theory of mixture, ASME Int. J. Eng. Sci. 18 (1980) 1129±1148. [6] S.S. Chen, N.T. Wright, J.D. Humphrey, Heat-induced changes in the mechanics of a collagenous tissue: isothermal, isotonic shrinkage, ASME J. Biomech. Eng. 120 (1998) 382±388. [7] A.E. Green, W. Zerna, Theoretical Elasticity, Oxford University Press, Oxford, 1968. [8] S.R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover, New York, 1984. [9] R. Haase, Thermodynamics of Irreversible Processes, Dover, New York, 1990. [10] H. Hormann, H. Schlebusch, Reversible and irreversible denaturation of collagen ®bers, Biochem. 10 (1971) 932± 937. [11] D. Jou, J. Casas-Vazquez, G. Lebon, Extended Irreversible Thermodynamics, Springer, New York, 1993. [12] A. Katchalsky, P.F. Curran, Nonequilibrium Thermodynamics in Biophysics, Harvard University Press, Cambridge, MA, 1965.
1556
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
[13] W.M. Lai, J.S. Hou, V.C. Mow, A triphasic theory for the swelling and deformation behaviors of articular cartilage, ASME J. Biomech. Eng. 113 (1991) 245±258. [14] M. Le Lous, F. Flandin, D. Herbage, J.C. Allain, In¯uence of collagen denaturation on the chemorheological properties of skin, assessed by dierential scanning calorimetry and hydrothermal isometric tension measurement, Biochim. Biophys. Acta 717 (1982) 295±300. [15] C.A. Miles, T.V. Burjanadze, A.J. Bailey, The kinetics of the thermal denaturation of collagen in unrestrained rat tail tendon determined by dierential scanning calorimetry, J. Mol. Biol. 245 (1995) 437±446. [16] I. Muller, T. Ruggeri, Rational Extended Thermodynamics, Springer, New York, 1998. [17] P.L. Privalov, Stability of proteins: proteins which do not present a single cooperative system, in: Advances in Protein Chemistry, vol. 35, Academic Press, New York, 1982, pp. 55±104. [18] K.R. Rajagopal, Multiple con®gurations in continuum mechanics, Reports of the Institute for Computational and Applied Mechanics, University of Pittsburgh, 1995, 6pp. [19] K.R. Rajagopal, A.R. Srinivasa, Mechanics of the inelastic behavior of materials. Part II: Inelastic response, Int. J. Plasticity 14 (1998) 969±995. [20] K.R. Rajagopal, L. Tao, Mechanics of Mixtures, World Scienti®c, New Jersey, 1995. [21] L.I. Sedov, A Course in Continuum Mechanics, vol. II, Wolters-Noordho, Groningen, 1972. [22] J.J. Shi, K.R. Rajagopal, A.S. Wineman, Application of the theory of interacting continua to the diusion of a ¯uid through a non-linearly elastic solid, Int. J. Eng. Sci. 19 (1981) 871±889. [23] L. Tao, K.R. Rajagopal, Flows of bubbly liquids, Applied Mech. Eng. 4 (1999) 407±442. [24] L.R.G. Treloar, The Physics of Rubber Elasticity, Oxford University Press, Oxford, 1958. [25] C. Truesdell, Rational Thermodynamics, Springer, New York, 1984. [26] C. Truesdell, R. Toupin, in: W. Flugge (Ed.), The Classical Field Theories, Handbuch der Physik, vol. III/1, Springer, New York. [27] W. Yourgrau, A. van der Merwe, G. Raw, Treatise on Irreversible and Statistical Thermodynamics, Dover, New York, 1982.
www.elsevier.com/locate/ijengsci
A mixture theory for heat-induced alterations in hydration and mechanical properties in soft tissues L. Tao, J.D. Humphrey, K.R. Rajagopal * Mechanical Engineering and Biomedical Engineering, Texas A&M University, College Station, TX 77843-3123, USA Received 18 October 2000; accepted 20 December 2000
Abstract We propose a framework to describe the mechanical behavior of thermally damaged biological soft tissues based on a melding of ideas from classical continuum thermodynamics and irreversible thermodynamics. Relevant assumptions are delineated, details of the constitutive formulation are provided, and a set of interfacial conditions are derived that can be used as the boundary conditions for certain diusion processes. The utility of the framework is illustrated by considering two problems that correspond to tractable experiments: transverse diusion through a ®nitely deformed membrane, and radial diusion in a thermally damaged cylindrical tissue. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Hydration; Mixture theory; Soft tissues; Thermal damage
1. Introduction Advances in laser, microwave, radio-frequency, and similar technologies continue to encourage the search for new applications of heat (i.e., thermal energy) to treat disease and injury. Examples of heat-treatment can be found in cardiology, dermatology, gynecology, oncology, ophthalmology, orthopedics, and urology, to name a few areas. The eects of increased temperature levels on cells and tissue are manifold, including the increased expression of heat-shock proteins, outright cell death, the denaturation of proteins, and localized edema. Associated with these changes, of course, are possible changes in the electrical, optical, mechanical, and thermal properties of the treated tissue. To date, however, there has been little attempt to quantify rigorously such changes in properties even though this information is essential for the analysis
*
Corresponding author. Tel.: +1-979-845-1251; fax: +1-979-845-3081. E-mail address: [email protected] (K.R. Rajagopal).
0020-7225/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 1 9 - 2
1536
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
and design of clinical heating strategies. Rather than address a speci®c clinical application or a particular data set from the literature, the purpose of this paper is to outline a possible approach for mathematically modeling the changes in hydration and non-contractile mechanical properties that arise due to heat-treatment of a soft collagenous tissue. Towards this end, we shall model a typical soft tissue as a ¯uid-infused solid that experiences ®nite deformations and exhibits nonlinear behavior as a continuum mixture. Moreover, we assume that severe heattreatment results in the coexistence of varying degrees of native and denatured collagen. Hence, a theoretical framework for modeling such a system must account for the behavior of a material consisting of at least three constituents, one of which (mobile ¯uid) is capable of motion with respect to the other two, which in turn behave as a constrained mixture (native and damaged tissue). Most heat-treated tissues and organs consist of parenchymal tissue (collagen, muscle, etc.) surrounded by delimiting membranes; examples range from tendons to muscular organs such as the heart and uterus. As such, these materials are necessarily non-homogeneous and the mechanical properties vary with the material point. Thermally-induced changes likewise vary from one material point to another. Moreover, given that the delimiting membranes serve as semipermeable barriers, heat-induced alterations in permeability must be accounted for as well. Although a complex problem, we can draw upon the many advances in the continuum theory of mixtures (e.g., see [3,4,20,24]) to formulate a general ®rst approach. We note, however, that certain modi®cations will be needed, and hence we will meld some ideas from classical continuum (e.g., [25]) and irreversible (e.g., [8,9,11,12]) thermodynamics. Our overall philosophy is that initial (empirical) observations are needed to delineate the general characteristics of the problem that must be modeled. Once done, one can then establish an appropriate (albeit perhaps simpli®ed) theoretical framework to guide the design and interpretation of speci®c experiments, which in turn yield constitutive relations for evaluation. For purposes herein, we shall rely primarily on observations from heat-treated tendons. Reasons for this are twofold: ®rst, much of the relevant data in the literature has come from rat tail tendon or chordae tendineae; second, being thin cylindrical structures consisting primarily of collagen, these tissues admit tractable initial-boundary value problems for study both analytically and in the laboratory. 2. Background Collagen is the most abundant protein in the body; it exists in at least 15 dierent types, each similar but individually specialized for a speci®c function. Herein, we will be concerned primarily with the ®brillar type I collagen and the meshed type IV collagen. The former is the primary constituent in the ``core'' of a tendon, endowing the tissue with marked tensile strength, and the latter is the primary constituent of the surrounding delimiting membrane or ``sheath''. Collagen molecules are composed of three polypeptide chains, each having some 1300±1700 amino acid residues, the majority of which (1000±1400) are organized into a central triple helical structure [2]. The latter is on the order of 285 nm long and 1.5 nm in diameter and consists of a repeating triplet of amino acids, (G±X±Y)n , where G is glycine and X and Y are often proline or hydroxyproline. The collagen molecule is stabilized by extensive hydrogen
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1537
bonds, with associated water playing a key role. In ®brillar types, the collagen molecules form a quarter-staggered arrangement (which gives rise to the characteristic 67 nm bandedness), with lysine and hydroxylysine residues forming stable intra- and inter-molecular cross-links particularly in the non-helical end domains. In contrast, type IV collagen is stabilized via extensive disulphide bonding. As a result, the activation energy of type IV collagen is about twice that of type I [15]. The eects of heat on collagen can be ``reversible'' or ``irreversible''. For example, moderate heating can result in a local unfolding within the protein that is reversed upon the restoration of normal temperatures. Severe heating results in a time-dependent irreversible transformation of the native triple helical structure into a more random (coiled) structure [15,17]. It is thought that the latter transformation occurs primarily via the breaking of long sequences of hydrogen bonds that stabilize the triple helix; particular sub-domains along the molecule (e.g., between residues 877 and 936; [15]) may be more susceptible to the breaking of such sequences. Heat-induced breakage of reducible cross-links [14] may also play a role. Regardless of the precise mechanisms, on a gross scale denaturation results in an irreversible shrinkage of the tissue, which for us becomes a convenient metric of thermal damage. This shrinkage may result due to exposure of hydrophobic regions of the molecule that are initially within the center of the helical structure. Regardless, Hormann and Schlebusch [10] showed that gross shrinkage provides much the same information on denaturation kinetics as many biochemical assays (though not the same detail on the underlying mechanisms of course). Denaturation also results in a gross change in the hydration of the collagen, which is thought to be due ®rst to the liberation of water and then a subsequent increased absorption of water via water bridges. Again, independent of precise molecular mechanisms, on a gross scale heating and post-heating recovery involves changes in hydration. Much work has been done on the diusion of a ¯uid through an elastic solid ([5,13,20,22], to mention just a few). These studies have been mainly concerned with the application of mixture theory within the context of speci®c models, however, and none consider the unique features of bio-thermomechanical problems. In this paper we revisit the diusion problem, speci®cally the diusion of water through a thermally damaged biological soft tissue that experiences ®nite deformations. We shall treat the body as a mixture of mobile water with weak hydrogen bonding and a porous solid with tightly bonded water molecules. Further, upon damage, some hydrogen bonds are broken during the heat treatment, and the associated region of the native tissue ``transforms'' into an altered material. The theory of interacting continua [1,4,20,25] provides a natural setting for the study of such problems. As noted above, we will infuse certain ideas from irreversible thermodynamics into the classical theory of interacting continua. We feel that such a thermodynamic setting provides a reasonable basis for studying the problem at hand. While it is somewhat formal in its approach, in view of the lack of a clear understanding of the basic mechanisms (physical and physiological) that govern the problem, such a methodology can lend insight into the various issues that have to be considered further in the theory. We will state explicitly the assumptions that are introduced in the analysis, so that we can recognize the status and relevance of these assumptions, knowledge of which may lead to possible modi®cations and improvements in further studies. A feature of the present analysis is the formulation of a set of interfacial jump conditions that can be adopted as boundary conditions in describing some speci®c diusion processes in tractable experiments.
1538
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
3. Preliminaries The body consists of mobile water, loosely bonded with the tissue, and a porous solid having tightly bonded water molecules. The tissue may be partitioned further into two sub-components: one associated with that part wherein the hydrogen bonds are broken due to thermal damage and another that is the intact native tissue. It is expected that the damaged and the native constituents intermingle (over representative volume elements), with inter-conversion occurring from damaging and healing. We shall suppose that the solid constituents undergo at a point the same motions in a homogenized sense. Within the context of mixture theory, a set of thermodynamic quantities can be assigned to each component: quantities associated with the mobile water and the tissue will be denoted by the superscripts f and t, respectively, and quantities associated with the damaged and native subcomponents will be denoted, respectively, by the superscripts d and n. Thus, the following equations describe the balance of mass, linear momentum, and energy for each component a, oqa div qa va ma ; ot X
ma 0;
a f ; n; d;
1
mf 0;
2
n;d
o a a T q v div qa va va div Ta pa ; ot
a f ; t;
3
qt qn qd ; X
4
pa 0;
5
f ;t
o a a 1 a2 1 a2 a a a q jv j div q jv j v div Ta va ot 2 2
qa pa va sa ;
a f ; t; 6
X
pa va sa 0:
7
f ;t
Here, qa , va , Ta , a and qa are, respectively, the mass density, velocity, Cauchy stress tensor, internal energy density, and heat ¯ux associated with the ath component; ma , pa and sa denote, respectively, the mass density conversion rate, the internal momentum supply, and the energy supply to the ath component resulting from the interactions between the components.
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1539
It should be emphasized that no electro-magnetic eect is included explicitly here, though hydrogen bonding can occur as follows. The external electro-magnetic ®eld is assumed negligible, the component f is neutral, and the component t is also neutral since the charge of the porous solid and the charge of the bonded water molecules supposedly neutralize each other due to the bonding. If the charge per unit mass of the solid is treated as constant in time, the charge distribution c at time t throughout the solid is determined completely by the charge distribution c0 in any reference con®guration j and the deformation of the solid. Therefore, the charge c is not needed explicitly in the formulation provided that the eects of c0 on ma , Ta , pa and so on are properly accounted for. Next, we infuse into our framework a thermodynamic notion that is widely used in nonequilibrium thermodynamics. The entropy density ga for component a can be introduced, for example, under the hypothesis of local equilibrium, and a balance equation can be proposed that has an analogous structure to the balance equations of (1), (3) and (6) [8,27], o a a q g div qa ga va div Ja ra ; ot
a f ; t;
8
where Ja and ra are, respectively, the entropy ¯ux and the internal entropy production rate per unit volume. The separation of the right-hand side of (8) into two parts leads to a natural representation due to the restrictions imposed by the second law of thermodynamics for isolated systems, as demonstrated below. Several observations are noteworthy here. First, by adopting the above equations, we tacitly restrict ourselves to motions without possible rapid and random ¯uctuations about some sensible mean value; as such, no equations need be posited for ¯uctuations. Examples of such ¯uctuations are the presence of turbulence or discontinuous distributions of a component a, which can be severe if the characteristic dimension of the internal geometric structure of the physical components, like pore size, is large. Neglecting ¯uctuations in the present formulation is justi®ed because the diusion process is slow (on the order of hours), and the water and tissue components are distributed regularly within an internal geometric structure that has a small characteristic dimension l. Second, the above equations are local. Introduction of continuously dierentiable ®elds qa , va , a J , ga , ra and so on is based on the assumption that each component can be modeled as a continuum and the hypothesis of local equilibrium, that is, these ®elds are averaged over some characteristic volume V or area, large compared with the characteristic dimension l. In formulating governing equations, independent of primitive quantities that are self-evident but unde®ned, we nonetheless have to motivate some inferential method for quantifying them. For instance, whereas temperature cannot be measured directly, we indirectly infer it by recognizing that the change in the length of a mercury column bears a one to one relationship with changes in temperature. The diculty with postulating new governing equations such as (8) is that, at this moment, we have no means to motivate a constitutive structure for quantities such as Ja or an inferential basis for measuring them. Nonetheless, as we shall see, such formalism does lead to ways for developing certain meaningful restrictions on the class of processes that we are interested in. With this understanding, we proceed to derive a restriction on ra by applying (8) to an isolated mixture U with U V,
1540
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Z
o r dv ot U
Z X U
qa ga dv P 0;
f ;t
r :
X
ra ;
9
f ;t
where the inequality required by classical thermodynamics guarantees a non-decreasing entropy in an isolated body. Next, taking into account the relative arbitrariness of U ( V), treating V as small enough due to the smallness of l, and assuming that the form of r is independent of whether the mixture in motion is isolated or not, we get X ra P 0; 10 r f ;t
for the motion of the mixture in general. Third, the hypothesis of local equilibrium may allow us to apply the concept of ga and (8) to a motion of the mixture away from a speci®c equilibrium. The reasoning is that in each local volume and during a short interval of time, the local motion is not far away from a corresponding equilibrium state, and ga and (8) can be de®ned according to the hypothesis of local equilibrium. If the motion is far away from any special state of equilibrium, then, unlike in irreversible thermodynamics, we do not necessarily need to introduce linear forms for the constitutive formulation based on any single equilibrium state, as is usual in rational and extended thermodynamics. Fourth, ga as well as a will be treated as a derived quantity completely determined by quantities such as mass density, velocity, deformation gradient, temperature and so on. That is, ga as well as a are treated as a state variable as in classical thermodynamics. Consequently, Eqs. (8) and (10) can be linked directly to the balance equations (1)±(7), and some restrictions can be imposed on the forms of ma , qa , pa , Ta , and Ja following the standard practice in the thermodynamics of continua [8,16,25]. 1 We are aware of the uncertainties associated with such a scheme, or for that matter any other scheme adopting an entropy balance equation to obtain the structure of undetermined quantities contained in the other balance equations. For instance, there is ¯exibility, to some extent, in proposing the forms of the entropy ¯ux and the (non-negative) internal entropy source, which is due partly to our lack of knowledge on how to quantify and/or measure entropy in a complex material undergoing a complicated motion. This uncertainty thus requires that the results from this scheme be tested experimentally, if possible. Finally, consider some issues related to the solid component. We have assumed that the two sub-components of the tissue have individual mass densities but undergo the same motions, i.e., qt qn qd ;
vn vd vt :
11
Then, from Eqs. (1) and (2)1 , we get oqt div qt vt 0: ot 1
12
When the ¯uctuations are not negligible, (10) may not play the same role as in rational thermodynamics or irreversible thermodynamics on the closure of the above equations since an additional assumption is needed, as in the ¯ows of bubbly liquids [23].
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1541
If we introduce a mass fraction for the damaged solid, m
qd ; qt
13
Eqs. (1) and (12) yield o md t v grad m t : ot q
14
Eq. (12) can be integrated to yield qt det F qt0 :
15
Here, qt0 is the eective mass density of the tissue component in the original reference con®guration j0 of the tissue, and F is the deformation gradient of the solid with respect to j0 . This j0 might be identi®ed with the state of the tissue in which qd 0 or qd is negligible compared with qn and it is (almost) stress free. As long as there is no damage, the response of the material is elastic. As the thermal damage takes place, qd increases, and the response is no longer purely elastic. The sub-component of damaged constituent will have a dierent reference con®guration jd from which its response is elastic. Motivated by the idea of multiple natural con®gurations [18], we introduce a reference con®guration j that is intermediate between j0 and jd and may evolve due to thermal damage; we can then decompose the deformation gradient F into a part Fj that is associated with an elastic response and another that is associated with the evolution of natural con®gurations. If the con®guration j is determined by a quantity G, a mapping from j0 to j relates F and Fj through Fj FG 1 ;
16
plus an equation for the evolution of G. 4. Speci®c model To apply the above framework to study the diusion of water through a thermally damaged tissue, we de®ne here q :
X
qa ;
f ;t
1X a a q; : q f ;t La : grad va ;
1X a a 1X a a q v ; g : qg; q f ;t q f ;t X q : qa ; ua : va v;
v :
f ;t
d o : v grad: dt ot
Then Eqs. (1)±(7), (8) and (14) yield
17
1542
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
dqa dt
ua grad qa
dm dt
ut grad m
qa tr La ma ;
md ; qt
q
d X tr Ta La dt f ;t
q
X dg Ja div dt f ;t
18
19
pa va
div qa qa a ua ;
qa ga ua r:
20
21
We now assume that ``both'' components ± mobile water and tissue ± have the same absolute temperature ®eld T, and the entropy and the internal energy of the mixture are completely determined by the primary ®elds qf , qd , T, Fj , and G and the reference con®guration j fg; g fg; g qf ; m; T ; Fj ; G; j
22
which means that the entropy and the internal energy of the mixture are modeled similar to those for an elastic process with deformation Fj with respect to j for ®xed qf and m. The dissipative response due to the viscosity of water will be accounted for through the internal momentum supply pf to be introduced and discussed later. We should mention that this treatment is intended to apply only to the diusion process considered here, and other quantities might have to be included if dierent processes of the mixture are to be considered. For example, g and may depend on the velocity gradient Lf if visco-elasticity is signi®cant in the process under consideration. We restrict ourselves to a special class of motions for which the assumption (22) is reasonable. Introducing the Helmholtz potentials wa : a
1X a a qw; q f ;t
23
T g w qf ; m; T ; Fj ; G; j :
24
T ga ;
w :
we have w
Substituting (24) into (21) and using (18) and (20) results in
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
"
1543
!
X 1 a ow dT ow a a a a T r T div J q q w u q g tr Tt qFj T Lt T oT dt oFj f ;t # ow f f f ow f f f f ow T qq I L p q gradw q f gradq q gradT vt f oq oq oT X 1 ow md ow T T ow T qa qa a ua grad tr q F G q j T om qt oFj oG f ;t o : vt grad GT P0: ot
vf
25
Moreover, we will assume that the tissue and the mobile water are intrinsically mechanically incompressible (though thermally compressible) and that the volume additivity constraint holds, i.e., qf qt 1; qfR qtR
26
where qfR and qtR ( qt0 = 1 b0 ) are the physical mass densities of the water and tissue in j0 , respectively, and b0 is the volume fraction of mobile water in j0 . This relation can be partially justi®ed by the fact that the tissue component includes plentiful water molecules bonded tightly with the solid. With the presence of these bonded water molecules, the tissue component and the mobile water component may avoid `void ®lling' at the molecular level, and consequently volume additivity holds. Next, in view of the range of temperature changes over the processes modeled here, it is appropriate to treat qfR and qtR as constant. It then follows from (26), (1), and (15) that t q grad tR vt q
f
v
qt qf tr tR Lt fR Lf q q
0:
27
Since Eq. (25) is expected to hold under the above constraint, a Lagrange multiplier p is introduced to enforce the constraint thus leading to (" X 1 a ow dT a a t a a fR ow J q q w u q g tr T p q q f T r T div T oT dt oq f ;t # ) f qt ow qf q tR I qFj T Lt Tf p fR I Lf pf pgrad fR qf gradw q q q oFj X 1 ow md f ow t f a a a a gradT v v T q q q u grad q oT T om qt f ;t ow o T T ow t 28 tr q Fj G v grad GT P0: oFj oG ot
1544
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Here we have used p qfR q ow=oqf instead of p in the mathematical treatment for the reason to be clear below. With the inclusion of qf in w, the Helmholtz potential should be independent of det Fj G det F on the basis of (15) and (26), i.e., w w qf ; m; T ; Fj ; G; j with det Fj G absent. Eq. (28) can be simpli®ed if we choose Ja
1 a q qa wa ua T
29
and g
ow ; oT
30
which are the usual assumptions that are made. This eliminates the divergence term and the time derivative of the temperature term, thus Eq. (28) yields, (" # ) t q ow qf t t f fR ow I qFj T L T p fR I Lf T r tr T p q q f oq qtR q oFj f X q ow pf p grad fR qf grad w qf qa qa a ua grad T vt vf T oT q f ;t d 1 ow m ow ow o grad tr q FTj G T 31 q vt grad GT P 0: T om qt oFj oG ot Inequality (31) places some restrictions on the structures of the functions pf , Tt , Tf , q, md and G. Notice that there is no unique way to arrange the cross terms involved in the velocity dierence and the gradient of temperature; this may be immaterial if the cross eects are handled properly through the appropriate choice of constitutive relations for pf and q. Another issue is that (29) and (30) are dictated here by the elimination of the divergence and the time rate of change of temperature and they should be considered as assumptions. 2 Relation (29) will be used later in generating the interface conditions in the isothermal process, which may provide a way to check its validity indirectly through experiments. Now we propose a set of constitutive relations, which reduce to that proposed by Shi et. al. [22] in the case of isothermal processes without damage, t q ow t fR ow p q q f I qFj T ; 32 T tR oq q oFj
2 It is possible, for example, to add a term like Ea ua to the right-hand side of (29) for Ja where Ea is a quantity to be determined. This speci®c addition will aect the forms of pf , Tt and Tf , but will not aect the balance of linear momentum equations in the case of symmetric Ea for isothermal processes, a result similar to that of Shi et. al. [22] which was proposed from a dierent perspective.
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Tf
p
qf I; qfR
f q p p grad fR q f
X f ;t
md
33
qf grad w qf
ow grad T b vt oT
vf ;
1 q q u C grad ; T a
cqt
1545
34
a a a
35
ow ; om
36
o ow T t G v grad G b FTj ot oFj
T
ow : oG
37
Here, c P 0, b P 0, b and C are, respectively, the drag coecient and heat conductivity, general positive de®nite second-order tensors required by the inequality (31) and to be determined by experiments. The reason not to treat b simply as a multiple of the unit tensor is due to the possible anisotropy of the pore structure of the tissue. Physically, b characterizes the momentum exchange (due to viscosity) between the viscous water, which moves relative to the solid constituents, and the tissue; it is important to retain this term even though the diusion is slow since the `channels' for the diusion are very small and the contact area per unit volume between the water and tissue is very large. This term should play such a role more eectively than a relation for Tf which includes the viscous eect as an internal source term in the equations of motion, since the standard form of Tf , in terms of the velocity gradient in the ¯uid, may not describe the eects due to interactions in the pores adequately. A simple example to illustrate this inadequacy is to consider a diusion velocity ®eld of the form (vf1 , vf2 , vf3 vf x1 ; 0; 0. This velocity ®eld would not yield f if the standard velocity gradient form were employed, yet it would give rise any viscous term in T12 to a viscous drag locally. Next, associated with (36) and (37), an activation criterion like a yield condition in plasticity could be introduced if necessary that is based on temperature. Finally, it is worth mentioning that we can also obtain the above relations on the basis of maximizing the entropy production r, following Rajagopal and Srinivasa [19]. The form of w in the above relations can be restricted further by invariance requirements. Here, we require Galilean invariance [21,25], instead of the more strigent and commonly used assumption of material frame indierence (MFI). (The issue regarding MFI will be discussed in detail elsewhere.) This invariance requires that the form of w be invariant under the coordinate transformation x Qx ct a; where Q and c are constant with QT Q QQT I. Therefore,
1546
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
w qf ; m; T ; Fj ; G; j w qf ; m ; T ; Fj ; G ; j w qf ; m; T ; QRj Uj ; G; j;
38
with the help of qf qf , m m, T T , G G, Fj QFj , Fj Rj Uj , where Rj is proper orthogonal and Uj is symmetric and positive de®nite. Now pick Q RTj t0 for some instant t0 and (38) yields w qf ; m; T ; Fj ; G; j w qf ; m; T ; RTj t0 Rj Uj ; G; j
8t;
and thus w qf ; m; T ; Fj ; G; j w qf ; m; T ; Uj ; G; j at t0 : The arbitrariness of t0 results in w depending on j, qf , m, T, Cj U2j FTj Fj and G but not through det Cj GGT det Fj G2 ). Additional restrictions on the form of w can be obtained if the symmetry of j is taken into account. As a proper approximation for the processes modeled here, we may assume that the previously introduced reference con®gurations j0 and jd have the same symmetry (G0 Gd ) or that jd has a dierent symmetry (but G0 Gd ). For example, if the tissue is assumed to be orthotropic with symmetry planes X1 constant and X2 constant in j, w may be taken as [7], 2
2
2
w qf ; m; T ; Cj ; G; j w qf ; m; T ; Cj11 ; Cj22 ; Cj33 ; Cj13 ; Cj32 ; Cj21 . . . ; j:
39
Alternatively, if the tissue is assumed to be transversely isotropic with its preferred direction along the X3 direction in j0 , w will take the following form, 2
2
w qf ; m; T ; Cj ; G; j w qf ; m; T ; I; II; Cj33 ; Cj13 Cj32 . . . ; j;
40
where I : tr Cj ;
II :
1h 2 tr Cj 2
i 2 tr Cj :
41
The ®nal, speci®c form of w needs to be ®xed through experiments. To have a better understanding of the role of w in the process treated here, consider the case of isothermal slow diusion, which occurs during the ``recovery'' process post-heating [6]. With this premise, the inertial term in the balance of linear momentum for the water component f is likely negligible, and Eq. (3) (with a f ) reduces to vf
vt
qf 1 b grad p qfR w : fR q
42
Eqs. (33) and (34), with grad T neglected, have been used in the derivation. This relation indicates that there are two parts aecting the isothermal slow diusion of water: the ®rst part is associated with the gradient of the spherical stress of the water component (p), similar to the commonly used Darcy's law; the second part is associated with the gradient of w, similar to Fick's law, which accounts for variations in the eective mass density qf (concentration) and the
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1547
deformation F. Also, the inverse of the drag coecient b 1 is part of the diusion coecient, which further shows the necessity for treating b as a general positive de®nite second-order tensor in the case of anisotropic tissues. By identifying that the spatial variations of both the spherical part of the stress of the water component and the Helmholtz potential are the driving forces for the diusion, we need to develop proper boundary conditions, an issue to be discussed in the next section. 5. Interfacial conditions In this section, we discuss interfacial jump conditions that are necessary to study the diusion process. We will deal with three classes of interfaces: (1) mixtures present on both sides; (2) a mixture on one side and pure water on the other side; and (3) a mixture on one side and water within a rigid porous solid on the other side. Consider an interface S in a mixture. On each side of S, denoted by the superscripts and respectively, there may be mobile water (denoted by the superscript f) and tissue (denoted by the superscript t). The two tissue components may dier from each other or one component may be absent. The interface links the two tissue components perfectly if both are present, that is, x Xt ; t x x x x Xt ; t always holds, if it held initially, where Xt denotes the location of particles from both tissues. Further, we assume that there are no physical properties attached to the interface except those derivable from the components on both sides of S while x ! x , which implies that there is no accumulation of mass, linear momentum, energy or entropy on the interface. Nonetheless, this interface is a material surface since it is composed of the same particles either from the contact of two types of tissues or from one side of the interface being occupied by a tissue and the other side free. Consider a ®eld quantity /a , associated with component a on either side of S, which can be characterized by the following balance equation: o a a q / div qa /a va div pa ga ; ot
43
where pa and ga are, respectively, the ¯ux and the source related to /a . Choosing a control volume around S, applying the balance of /a for the mixture to the control volume, and shrinking the control volume, we obtain [20,26] f f f q / v
vt
pf pt
n qf /f vf
vt
pf pt
n;
44
where n is the normal to the interface. This relation provides the required interfacial jump conditions by assigning proper physical quantities for /a and pa . Class 1. For an interface with tissues present on both sides, picking, respectively, (i) /f 1, T T f p pt 0 (balance of mass) and (ii) /f vf , pf Tf , pt Tt (balance of linear momentum), we obtain from (44) f f q v
vt
n qf vf
vt
n;
vt vt ;
45
1548
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
" f f
q v v
f
t
v n
X
#
"
aT
f f
T n
q v v
f
t
X
v n
f ;t
# aT
T n
:
46
f ;t 2
Further, picking, respectively, (iii) /f f jvf j =2, pf Tf vf qf , pt Tt vt qt (balance of energy) and (iv) /f gf , pf Jf qf qf wf uf =T , pt Jt qt qt wt ut =T (entropy equation), and combining the two resultant jump conditions to eliminate the heat ¯ux qf qt , we obtain " # X 1 2 qa qa a ua Ta va n vf vt qf jvf j 2 f ;t " # X 1 2 qf jvf j qa qa a ua Ta va n vf vt 2 f ;t
47
and " 1 f 2 f vf q w jv j 2
v
t
X
# a a
Tv
f ;t
" 1 f 2 f n q w jv j vf 2
v
t
X
# a a
Tv
n:
48
f ;t
Class 2. Similarly we can obtain the jump conditions across the interface, on the side ()) where the mixture is present and on the other side () occupied by the pure water, taking into account that Tf P I where P is the water pressure (the viscous eect in the pure water region is ignored since the ¯ow of water is slow and uniform and thus the velocity gradient is very small),
qf vf
" qf vf vf
vt
n qfR vfR
t
v n
X
vt #
TaT n
n;
vt vt ;
qfR vfR vfR
49
vt n P n ;
50
f ;t
" # X 1 qf jvf j2 vf vt qa qa a ua Ta va n 2 f ;t 1 fR 2 fR fR fR t fR fR q jv j n: v q Pv v 2 " 1 f 2 f vf q w jv j 2
v
t
X f ;t
# a a
Tv
fR
n q
1 w jvfR j2 2 fR
51 fR
v
t
v Pv
fR
n:
52
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1549
Here fR and wfR are, respectively, the internal energy and the Helmholtz potential for pure water. Class 3. In the case that on the side ()) of the interface the mixture is present and on the side there is water and a static rigid porous solid through which the water ¯ows, we need to modify the above relations. We will assume that the rigid solid has such a large porosity a0 that it has negligible eect on the ¯ow of the water compared to the tissue, but the pores are small enough that the geometry of the interface is completely determined by the rigid solid. Also, we assume that the rigid solid will exert little eect on the tissue along the directions tangent to the interface. We can take qf a0 qfR , f fR , wf wfR , Tf a0 P I, qf a0 qfR and qt qs where qs is the heat ¯ux associated with the rigid porous solid, and obtain, by modifying (49)±(52) f f n a0 qfR vfR n; q v
53
" # 1 f 2 f X a a a f a a a q q u T v n q jv j v 2 f ;t 1 2 a0 qfR fR jvfR j vfR a0 qfR qs a0 P vfR n; 2 1 f 2 f f q w jv j v 2
f f
T v
54
1 fR 2 fR fR fR fR n a0 q w jv j v a0 P v n: 2
55
The relation associated with the balance of linear momentum is not useful here since the stress inside the rigid porous solid is indeterminate. We should mention here that, in the case of isothermal processes, the interfacial conditions of (47), (51) and (54) associated with the balance of energy will be left out. 6. Applications 6.1. Diusion through a membrane Consider the slow isothermal diusion of water through a thin membrane of original thickness H without thermal damage. Suppose that the membrane is orthotropic, with X1 0 and X2 0 as the symmetry planes, and 2
2
2
w w qf ; T ; C; j0 A qf ; T ; C11 ; C22 ; C33 ; C13 ; C32 ; C21 ; j0 :
56
While laid on a rigid porous plate at X3 0, the membrane undergoes the deformation x1 kX1 ;
x2 kX2 ;
x3 f X3 ; t;
where k is a constant. The water moves according to
57
1550
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
vf 0; 0; vf X3 ; t ;
58
under the pressure dierence PH P0 P 0, where PH and P0 are, respectively, the pressure of water immediately above and below the membrane. A simple calculation yields 2 2 3 2 3 k 0 0 0 0 k of : 59 F 4 0 k 0 5; C 4 0 k2 0 5; g : oX3 0 0 g 0 0 g2 Balance of mass and the volume additivity constraint give 1 qt qt0 k2 g ;
60
:
61
f
q q
fR
1
qt0 2 kg qtR
1
Balance of linear momentum for the mixture as a whole and for the water component yields, o f t T33 T33 0; ox3
62
f 1 o q p qfR w b33 fR vf ox3 q
of ot
0:
63
Here we have neglected inertia and assumed that Tt , Tf and pf are described by (32)±(34) and (56) with b13 b23 0. The boundary conditions for these two equations can be derived from (49) to (55) and the ¯ow conditions mentioned above: f jX3 0 0;
p qfR w p qfR w f t T33 T33
64 P0 qfR wfR ;
65
X3 H
PH qfR wfR ;
66
X3 H
67
X3 0
PH :
Solving the above equations, we get qt ow ow 2 q w q tR f 2g fR wfR q oC33 X3 0 q oq
PH P0 ; qfR
68
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
Z
H
0
b33
qf qfR
1
f
v
of ot
1551
PH
g dX3
P0 :
69
Moreover, taking into account the thinness of the membrane, we can obtain the mass ¯ux from (69) F0
q v
"
f f
X3 0
qf 2 gHb33
# X3 0
PH P0 ; qfR
70
with g to be determined from (68) and qt and qf from (60) and (61). We may interpret gHb33 in (70) as the drag coecient for the membrane as a surface. 6.2. Radial diusion in a membrane-covered cylindrical tissue Consider a cylindrical tissue with a membranous covering, with the core tissue occupying the region 0; RI and the membrane occupying the region RI ; RO in the reference con®guration j0 . The core is assumed to be a mixture of a transversely isotropic tissue and water with its preferred direction along the axial direction Z in j0 . Suppose that the membrane is orthotropic with symmetry planes H and Z. Initially the tissue is at temperature T0 . Then, suppose an axial load W is applied to the tissue which is then immersed into a water bath at temperature Th (> T0 ) and pressure P0 . Due to the overall thinness (RO small), let the whole tissue reach the equilibrium temperature Th rather quickly. This heating causes thermal damage to the core tissue, but the damage to the membrane is relatively small and may be neglected (this restriction can be relaxed easily in our formulation). Also, because of this heating, there is water diusion through the membrane and into or out of the tissue (both during and following heating). Thus, we have w w qf ; m; T ; Cj ; G; j
71
for the core tissue and w w qf ; T ; CRR ; CHH ; CZZ ; CRZ 2 ; CZH 2 ; CHR 2 ; j0
72
for the membrane. We shall assume that the temperature ®eld is approximated by T
T0 ; Th ;
t < 0 and t > th ; t 2 0; th ;
73
that the tissue and the membrane undergo the deformation r f R; t;
h H;
z k tZ;
that water diuses according to
R 2 0; RI [ RI ; RO ;
74
1552
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
vfr vf R; t;
vfh 0;
_ vfz kZ;
R 2 0; RI [ RI ; RO ;
75
and that the core tissue is damaged due to heating, thus giving rise to 2 3 g R; t 0 0 G4 0 h R; t 0 5 0 0 k R; t
76
in a cylindrical coordinate system. It follows that for R 2 0; RI 2 of
0
oR
F4 0 0
f R
0
2
3
0 0 5; k
1 of g oR
6 Fj 4 0 0
2
3 0 7 0 5;
0 f hR
6 Cj 6 4
k k
0
1 of g oR
0 0
2
0
f 2 hR 0
0
3
7 7 0 5; k 2
77
k
and for R 2 RI ; RO 2 of
0
oR
F4 0 0
f R
0
3 0 0 5; k
2 6 C4
of 2 oR 0 0
0 f 2
R
0
0
3
7 0 5: k2
78
Next, a straight-forward but lengthy calculation yields 1 f of q q k ; R oR t
t0
f
fR
" q q
1
1# qt0 f of ; k R oR qtR
q qf qt ;
vf
qfR qf
p p^ R; t
79
80
81 ! qt of 1 k_ f ; qtR ot 2 k 1 fR 2 q kkZ ; 2
in 0; RI [ RI ; RO , and
82
83
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
om ot
c
1553
ow ; om
84
o p^ qfR w oR
q
fR
of ovf oR ot
ovf br of oR qf oR
f
v
of ot
;
85
f 1 t of o2 f of ovf ovf f of q of of fR q f q q tR v br v q oR ot oR qfR oR ot2 ot oR ot " " # 2 2 t o q ow 1 of ow 2q of 1 of ow fR fR q w q q tR f 2q oR g oR oCjRR f oR g oR oCjRR q oq t
f hR
2
ow oCjHH
# 0;
86
" 2 og 2 of ow b 3 ot g oR oCjRR
# ow ; oGjRR
87
" 2 oh 2 f ow b 3 ot h R oCjHH
# ow ; oGjHH
88
2 ok 2k ow b 3 ot k oCjZZ
ow ; oGjZZ
89
in 0; RI , along with the axial, boundary and interfacial conditions # 2 ow qt k ow of p^ q q f f 2q dR oq qtR k oCjZZ oR 0 of W 2 ow 2qk RO RI ; f oCZZ oR R 2p
Z
RI
"
fR
p^
ow q q f oq fR
qt qtR 90
O
f jR0 0;
91
f jRR f jRR f jRR ; I
I
O
92
1554
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
1 f 2 f of v p^ q w v 2 ot RRI fR 1 fR 2 q of fR fR fR of P0 q w v v br v f 2 ot qf oR R f O f of ov of ov qfR RO RI ; vf oR ot ot oR R fR
of ot
93
O
"
2 # t q ow 1 of ow 1 of of qfR q tR f 2q qfR w vf 2 vf qf vf vf g oR oCjRR 2 ot ot q oq RI ( " # 2 2 1 fR 2 2q of of ow f ow v qfR wfR 2 f oR oR oCRR R oCHH R O ) t 2 qfR of of of ovf of ovf f fR q f t of o f br v f RO RI ; q tR v q ot ot oR oR ot2 q oR q oR ot
94
RO
"
2 qt ow of ow 1 fR q w vf 2 q q tR f 2q oR oCRR 2 q oq 1 fR 2 fR fR q w ; v 2 R
of v ot
fR
f
f f
q v
v
f
of ot
# RO
95
O
p^ q
vfR
fR
1 w vf 2 2
of v ot
f
RO
1 k_ : fj 2 k RRO
P0 q
fR
1 w vfR 2 2 fR
of v ot
fR
R O
;
96
97
tR fR We should mention in fR that the above derivation we have used the approximations q q 2 =2 . Also, an approximate treatment has been adopted for the motion and jp^ R; tj q kkZ of the membrane due to its thinness as in the previous case.The appropriate initial conditions are
f R; 0 f0 R; k 0 k0 ;
o f R; 0 0; ot
_ k 0 0;
98 99
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
m R; 0 0; g R; 0 1;
1555
100 h R; 0 1;
k R; 0 1:
101
To simulate the diusion process, we need to know the forms of w and br as well as the initial con®guration. This must await appropriate experiments, which may now be designed given this framework. 7. Summary We have developed a mathematical model for the diusion of water through a biological soft tissue that has suered thermal damage, based on a melding of ideas from classical continuum thermodynamics and irreversible thermodynamics. A set of interfacial conditions were derived that furnish boundary conditions that may be relevant to several classes of problems. Two speci®c processes of diusion, one for transverse diusion through a membrane and one for radial diffusion in a cylindrical tissue covered by a membrane, were considered because the predictions of the theory can be tested by experiments. Acknowledgements This work was supported in part by NSF grant (BES-9896166), NIH grant (HL-56135), and ARP grant (036327). References [1] R.J. Atkin, R.E. Craine, Continuum theories of mixtures: basic theory and historical development, Quart. J. Mech. Appl. Math. 29 (1976) 209±244. [2] S. Ayad, R.P. Boot-Handford, M.J. Humphries, K.E. Kadler, C.A. Shuttleworth, The Extracellular Matrix Facts Book, Academic Press, New York, 1994. [3] A. Bedford, D.S. Drumheller, Theories of immiscible and structured mixtures, Int. J. Eng. Sci. 21 (1983) 863±960. [4] R.M. Bowen, Theory of Mixture, in: A.C. Eringen (Ed.), Continuum Mechanics, Academic Press, New York, 1976. [5] R.M. Bowen, Incompressible porous media models by use of the theory of mixture, ASME Int. J. Eng. Sci. 18 (1980) 1129±1148. [6] S.S. Chen, N.T. Wright, J.D. Humphrey, Heat-induced changes in the mechanics of a collagenous tissue: isothermal, isotonic shrinkage, ASME J. Biomech. Eng. 120 (1998) 382±388. [7] A.E. Green, W. Zerna, Theoretical Elasticity, Oxford University Press, Oxford, 1968. [8] S.R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover, New York, 1984. [9] R. Haase, Thermodynamics of Irreversible Processes, Dover, New York, 1990. [10] H. Hormann, H. Schlebusch, Reversible and irreversible denaturation of collagen ®bers, Biochem. 10 (1971) 932± 937. [11] D. Jou, J. Casas-Vazquez, G. Lebon, Extended Irreversible Thermodynamics, Springer, New York, 1993. [12] A. Katchalsky, P.F. Curran, Nonequilibrium Thermodynamics in Biophysics, Harvard University Press, Cambridge, MA, 1965.
1556
L. Tao et al. / International Journal of Engineering Science 39 (2001) 1535±1556
[13] W.M. Lai, J.S. Hou, V.C. Mow, A triphasic theory for the swelling and deformation behaviors of articular cartilage, ASME J. Biomech. Eng. 113 (1991) 245±258. [14] M. Le Lous, F. Flandin, D. Herbage, J.C. Allain, In¯uence of collagen denaturation on the chemorheological properties of skin, assessed by dierential scanning calorimetry and hydrothermal isometric tension measurement, Biochim. Biophys. Acta 717 (1982) 295±300. [15] C.A. Miles, T.V. Burjanadze, A.J. Bailey, The kinetics of the thermal denaturation of collagen in unrestrained rat tail tendon determined by dierential scanning calorimetry, J. Mol. Biol. 245 (1995) 437±446. [16] I. Muller, T. Ruggeri, Rational Extended Thermodynamics, Springer, New York, 1998. [17] P.L. Privalov, Stability of proteins: proteins which do not present a single cooperative system, in: Advances in Protein Chemistry, vol. 35, Academic Press, New York, 1982, pp. 55±104. [18] K.R. Rajagopal, Multiple con®gurations in continuum mechanics, Reports of the Institute for Computational and Applied Mechanics, University of Pittsburgh, 1995, 6pp. [19] K.R. Rajagopal, A.R. Srinivasa, Mechanics of the inelastic behavior of materials. Part II: Inelastic response, Int. J. Plasticity 14 (1998) 969±995. [20] K.R. Rajagopal, L. Tao, Mechanics of Mixtures, World Scienti®c, New Jersey, 1995. [21] L.I. Sedov, A Course in Continuum Mechanics, vol. II, Wolters-Noordho, Groningen, 1972. [22] J.J. Shi, K.R. Rajagopal, A.S. Wineman, Application of the theory of interacting continua to the diusion of a ¯uid through a non-linearly elastic solid, Int. J. Eng. Sci. 19 (1981) 871±889. [23] L. Tao, K.R. Rajagopal, Flows of bubbly liquids, Applied Mech. Eng. 4 (1999) 407±442. [24] L.R.G. Treloar, The Physics of Rubber Elasticity, Oxford University Press, Oxford, 1958. [25] C. Truesdell, Rational Thermodynamics, Springer, New York, 1984. [26] C. Truesdell, R. Toupin, in: W. Flugge (Ed.), The Classical Field Theories, Handbuch der Physik, vol. III/1, Springer, New York. [27] W. Yourgrau, A. van der Merwe, G. Raw, Treatise on Irreversible and Statistical Thermodynamics, Dover, New York, 1982.