A non-linear viscoelastic constitutive equation for soft biological tissues, based upon a structural model

OOZl-9290/80,0701-0559

J. Bmnechmics Vol 13. pp. 559-W. @ Pergamon Press Ltd. 1980 Printed in Great Britain

$02.00/0

A NON-LINEAR VISCOELASTIC CONSTITUTIVE EQUATION FOR SOFT BIOLOGICAL TISSUES, BASED UPON A STRUCTURAL MODEL* W. F. DECRAEMER, M. A. MAES,V. J. VANHUYSEand P. VANPEPERSTRAFTE Laboratory for Experimental Physics, Rijksuniversitair Centrum Antwerpcn, Grocncnborserlaan, B-2020 Antwerpcn, BcIgium

171,

Aktract - A viscoelastic model based on the material structure of soft biological tissue is proposed and a corresponding non-linear viscoclastic constitutive equation is derived. It is shown how the relaxation spectrum H(r) may be derived using only the fundamental Fourierterm of the stress corresponding with harmonic strain oscillation. The ntmzrical determination of the parameters of the model, based upon a maximum likelihood method applied to a bivariatc normal distribution, is shor’tlydiscussed.

1.THE VBCOELAsIlC

MODEL

In a previous paper (Decraemer et al. 1980) we have shown that the ‘quasi-elastic’ stress-strain relation (measured under uniaxial tension) of some soft biological tissues (skin, papillary muscle, vein, fascia, tympanic membrane), can be described by a model based on the material structure. We assumed a large number (N) of purely elastic fibres (modulus k) embedded in a gelatin like liquid, all having the same cross-sectional area S but tierent initial length I,, normally distributed around a mean p with a standard deviation s. According to their initial length I,, we consider the fibres more or less folded. As those materials show a viscoelastic behaviour (even non-linear) instead of purely elastic, this model has to be adapted in order to describe viscoelastic phenomena. We assume therefore the existence of internal frictional forces between fibres or between fibres and the material in which they are embedded. We introduce this damping in the model by assuming that the fibres all have identical linear viscoelastic properties described by the relaxation function G(t). At the end of the paper we will illustrate the fitting of the model using some experimental data obtained in the course of a comparative test of different ways of preserving human tympanic membranes for later use as tympanic membrane prosthesis. The complete description of the experimental apparatus and the prep aration of the materials will be presented in a separate paper together with the results of this test_

moment T we give a length increase d&f) to a strip of material of initial length lo and actual length 1.Hereby we exert a force d&r) to stretch all fibres with initial length Ii < 1. Fibres with Ii > I+ dl will be somewhat unfolded, not yet stretched; we assume that this requires a negligible small force. dF(r) = G(O)

In this expression G(0) takes over the role of the elastic modulus k in the purely elastic constitutive equation (Dccraemer et aI., 1980). For t > T, the force increase dF(r) due to the length increase d](r), will have relaxed to ‘1 = ‘(*) d I(7)

-

dF(t,r) = G(t - 7) J

-

Ii

t

li

4-b

-(~-4)‘/2*

,j&.

(2)

x$ze

Integrating over the whole time history we find the total force F(t) at any time t 2 T to be

r I

F(r) = N

J-a

WI W - ~1~ U‘ -e-@-4~z~dl,dt.

’

(3)

Putting IL THE CONSTITUTIVEEQUATION A constitutive equation corresponding to this model can be deduced in the following way. Suppose at the

and making use of Leibnitz’s formula for the derivative of an integral, we find that dF=(r) -p-dr

* Receiwd 13 December 1978; in revised form 31 October 1979. 559

NG(0) d/(r) &s

dr

W.

560

F. DECRAEMER,M. A. MAE& V. J. VANHUYSEand P. VANHPERSTRAETE III. THE CONSTITUTW E EQUATION FOR HARMONIC SI-RAIN OSCILLATIONS

From (5) and (3), it follows that

F(t) =

GO- 4 dW3r)l G(O)7’

d7

(6)

Introducing (7) and dividing by A,, (the sample cross-section at rest) we finally come to the constitutive equation expressed in terms of Lagrangian stress (a) and strain (E) corresponding to our model

The function u’(w Fe/A,,) describes the immediate stress response, it is a purely elastic relation (index e) between stress and strain. The tune function G, [with

G,(O)= 1] is called the ‘reduced relaxation function’; it describes the relative variation of the stress following a step variation in strain. We write the constitutive equation (8) in this particular way to show that it is similar to the linear viscoelastic constitutive equation, in which Ehas been replaced by a function of &,namely a’(e). This is a wellknown method to introduce non-linear&r between u and e (Locket, 1972and Fung, 1972).It implies that the relaxation function G(c,t) may be written as a product of a function of s oniy and a function oft only. Indeed it is on the basis of a hypothesis of such a separation of the variables E(or 1 = e + 1) and t that (Fun& 1972) deduces a constitutive equation for soft tissues, similar to quation (8). In deriving equation (8) we worked quite ditTerently: we started from a mechanical model and derived the corresponding constitutive equation. This provided at the same time an analytical expression for d[Cr. equation (4X divided by A~]. The splitting of the relaxation function must always be seen as an approximation introduced to simplify the manipulation of the viscoelastic constitutive equation. This approximation will be all the better if the strain interval is not too extended. Pinto (1972) shows on the basis of stress relaxation data obtained on soft tissue (rabbit papillary muscle) that the separation of variables holds fairly well for E values up to 0.30. For experiments under small harmonic strain oscillations (with amplitudes cdup to a few%) on a constant prestrain level & the approximation is reasonable for measurements at the same E, but it is poorer when the measurements are taken at greater E,levels. At the end of section III we will discuss a criterium for estimating the goodness of the approximation. In the remaining sections we are focused on the viscoelastic behaviour under harmonic strain oscillation. In this light we will discuss the determination of the relaxation spectrum and the fitting of the model to experimental data.

As in most of our experiments we used harmonic strain oscillations, we describe the constitutive quation under this special strain history. Due to non-linearities of soft biological material, a harmonic strain oscillation of amplitude E,,, super posed on an initial strain E*

&PO

tso

e=.q+.sdeiru

t>o

(9)

gives rise to an anharmonic 6. Consequently, the corresponding stress will also be anharmonic. Expanding d and u as Fourierseries we may show (see Appendix A) that there exist a linear viscoelastic relations between corresponding terms of u and 6, all with the’same relaxation function G,(t). Consequently we only need one Fourierterm in order to determine G,, e.g. the easiest to measure, ul. Taking into account that & = 0, because U’(E)represents a conservative, purely elastic relation, we get from (A3) and (AlO) ~,(e,ri$e@++~~~~~~I = G+r(w)6i(++“.

(IO)

The complex modulus G,’ is often written as G: (w) = G; (w) + iG: (w)

(11)

where G; and G: represent respectively the storage and the loss modulus. From (A9)we find that G: and G: are given by G:(w) = G, + cm[G,(z) + - G,] sin wzdoz, I (12.1) and G:(w) =

[G,(z) - G,] cos wz dc~z.

(12.2)

From (10) and (11) we see that G; and G: are related to the amplitude ul and phase angle 4r of the fundamental term of u(t) as follows

Ul(E,W) G:(4 = -cos~1kw) b;(E)

(13.1)

cr(e, 0) G:’(w) = -sin&(e,w). 6(e)

(13.2)

Equations (13) also provide a necessary condition for the separation of the variables Eand t in the relaxation function. Indeed, dividing (13.1)by (13.2)we learn that cotg & -and thus &, which may be measured rather easily - must remain independent of e. By an analogue reasoning we find that this qually holds for the phase angles of the higher Fourier components of u. As the amplitudes of these components are only small fractions of the ground component uI (aI >>ut > uj ; higher components may be neglected)

the condition that r+i must be e-independent is the

561

A non-linear viscoclastic constitutive equation

most important. If this is not fulfilled, there is no doubt that the model will only give a poor approximation of the reality.

H(r) = -

dlnw

H(r) = IV. DETERMlNATlON OF THE RELAXATIONFUNCIION G(r) BY MEANS OF ITS RELAXATiONSPECTRUMH(r)

Expressing G(t) in the usual way as a function of its relaxation spectrum H(r), we have +C? H(r) e-“I d In T. G(r) = G, + (14) I0 Now the determination of G(t) can be replaced by the determination of H(r) and G,. In order to calculate the relaxation spectra we may apply the mechanism of the approximation formulae for linear viscoelastic material, using either the storage modulus G’ or the loss modulus G”. In Fig. 1 we represent u&,, and 4i in function of the logarithm of the strain oscillation frequency v, In v (or In o, with w = 271v) for an experiment on fresh human tympanic membrane; Fig. l(a) illustrates that for soft biological tissue, ei (In w) is an almost linear function (E, s,, held constant) of In w and from Fig. l(b) we learn that r$i remains small (& c 6”). Due to the small value of &, the uncertainty on @i has leas inknce on the accumcy of the calculation of G: with equation (11) than on G:’ so we will choose an approximation formula for H(r) which makes use of G:. Ferry (1970) reviews several of such approxhnation formulae. The most convenient for our purpose are the following formulae proposed by N. W. Tschoegl (see Ferry, 1970).

64-

111

f

I

dG

l/2 dZG

-

’ (dlnw)2

I/W=&

g& - (dl/2d2G

l/w-r,J2

(15.1)

(15.2)

Equation (15.1) holds when the slope of H(r) is positive, (15.2) when negative. From Fig. l(b) we learn that 0, remains sufficiently small to approximate cos$i by 1. Referring to (13.1) we see that by this approximation G; becomes proportional to ui (E, w) and thus also linearly related to ln w in the considered range. For a linear G&n w) function, the second derivative d2G;/(d In w)~, becomes zero and H,(T) takes the constant value dG&n w)/d In w, noted as R,. H,(r) = dG;/d In w = R,. (16) It is clear that this equation only holds in a r-interval corresponding to the experimental c+interval. The direct determination of the relaxation spectrum leads to a spectrum with constant intensity between limits rl and ?2 ; 71 corresponds with the highest experimental w and z2 with the lowest one. As relaxation phenomena described by time constants r outside this range have not yet or no longer influence on the steady state viscoelastic behaviour of the material, we may choose a zero intensity for the spectrum outside the interval (fi, TV).So we have K(r) = Ro H,(r) = 0

for r1 < r < 72 elsewhere

This choice corresponds to the hypothesis made by Fung (1972) - and adopted by Pinto (1972) - who states that ‘any fairly flat and broad continuous

11 E IM I

I’

I

2-

I

0 0.01

0.1

(17)

I 1

v4nz1

I 10

100

w

Fig. 1. Experimental data mprarenting CT& (Fig. la) and 41 (Fig. lb) in function ofla v for human tympanic membrane.

W. F.

562

DECRAEMER,M. A.

MAES.V. J. VANHUYSEand P.

spectrum’ will be able to describe relaxation phenomena similar to those registered for soft biological tissue. The block distribution H(T)has also been used for other relaxation phenomena due to internal damp ing, e.g. dielectrical, magnetical (see Neubert, 1963). Introducing the following block spectrum H(7) for the ordinary relaxation function G(r)

H(T)=c

for ri Q T d tz

H(r)= 0

elsewhere

1 _ .

knowledge of R,,. Calculation of G: (no) with the block spectrum yields 1 G:(nw) = 1 + : [ln(l + n2 0~~ri) 1 + clnlf I 71 -

ln( 1 + n2 o2 r:)]

. (20) I Pinto (1972) discusses how the parameters c, ri and r2 may be determined if out of relaxation experiments the reduced relaxation function G,(r) is known. Next we will show how these parameters may be calculated on the basis of stress-strain data for small harmonical strain oscillations. + ic[tg-’ nor2 - tg-‘nwr,]

(18)

the ‘intensity’ R,, of the corresponding spectrum for the reduced relaxation function G,(t) shows up to be

R. =

VANPEWRSTRAETE

(19)

1 + clnz

We will use in the following discussion the more fundamental parameterset (c, rr, r2) as this implies

V. NUMERICAL DETERMINATION OF THE PARAMETERS OF THE MODEL

Writing the elastic stress component

u; (E) as

4 (4 = s&(s)

Fig. 2(a). Experimental

theoretical

data of I I (u&,,) ~0s~~ and q = (a,/s,)sin 4, in function of Inv and their fitting for hutpan tympanic membrane. The parameter values of the adaptation tic: 01 = 2.60 10’ N/m’, 7, = 10m4s, r2 = 20.38 s and c = 0.0547. lO’N/d

Fig

2(b). Experimental data fitting from Fig. 2(a), alternatively represented by means of CT&, and phase angle +1 in function of In v.

(21)

A non-linear viscoelastic constitutive equation

and making use of (20) we may rewrite equation (10) as a,(@e@~

=

Ed

a 1 + cln?

‘51 x

1 +f[ln(l {

+~~r:)-ln(l

+w*r:)]

+ ic[tg- ’ 0.x2 - tg-’ c&J

(22) I

or

563

Ferry, J. D. (1970) Viscoelnstic Properties ofPolymers. Wiley & Sons, New York. Fung, Y. C. (1972) Stress-strain-history relation of soft tissues in simple elongation. Biomechics. Prentice Hall, New Jersey. Locke& F. J. (1972) Nonlinear Yiscoelastic Solids. Academic Press London. Neubert, H. K. P. (1963) A simple model representing internal damping in solid materials. Aeronaut. Quart. 187-210. Pinto, J. G. (1972) Mechanical properties of the papillary muscle in the passive and active state. Ph.D. Dissertation. University of California, San Diego. University microfilms. Ann Arbor, Michigan.

(23) APPENDIX

where r=

[a,Wcos $JI]/%

A

For a harmonica1 strain variation

(24)

and 4 = Cbr(e~)sin &I/Q. (25) In our experiments & and cdare held constant while the frequency is varied. Consequently a&Ed) stays constant. Referring to (21) we emphasize that tx is a parameter which plays the role of an incremental elastic modulus. The parameters a, c, rr and rz appear in both the real part (r) and the imaginary part (4). Therefore we fit the model to experimental data using a method of maximum likelihood applied to a bivariate normal distribution (for details see Appendix B). An example of a model fitting to experimental data obtained on fresh human tympanic membrane is given in Fig. 2. Tbe fit reveals to be better for T = (cJed) cosd~ thanforq = (ul/ad)sin#r [&isthephaseangle between the fundamental component of stress and a, see equation (A3)]. This is due to the fact that 4r is a small phase angle and therefore less accurately measurable. The same fitting procedure has been used with success on other soft biological tissues, s.a. human fascia and human vein. To determine the Fourierspectrum of u we had to restrict ourselves to 21 points per period as for manual data collecting this seemed to be a good compromise between accuracy and time consumption. For the time being we are switching over to computer controlled data acquisition, which will enable us to use more data points per period and to calculate with much more precision not only the amplitude and the phase of the fundamental stress component but evenso of the second harmonic (higher harmonics are of less importance due to their very small amplitudes) in order to check the model also for this Fourierterm. Preliminary measurements point in the direction that the relaxation spectrum, derived from the fundamental stress component, may also fit the second stress harmonic. A closer look upon this aspect of the model has been undertaken

E = & + cdeior

the functions d and (r may be, for sufficient large I, expanded as Fourierseries

n=1.2.

and U[E(‘)I = %(E

(A3)

The coefficients, which are functions of 4 and ~1,are noted as U,(E) and <(E). Bringing (A2) and (A3) in the constitutive equation (8), results in

(A4)

+

W

z n-1.

. .

-

7)

..J-=

d

,f(E)elt~+#$i.O)I d7 (AS) dr i The first integral represents a stress relaxation, corresponding with a step in u’ equal to tie&) occurring at I = 0, yielding x

r > 0.

G,(rk%e)

For steady state conditions this term takes the constant value large r

G, G(c)

646)

where G, represents the end-value of G,(r)- a monotonously decreasing function - for large r. In the other integral terms of (AS) we split up G,(r) in the constant term G, and a time varying term G,(t), so we may write the general term of the summation as

Making use of the substitution r we have for u(.F,r)

REFERENCES

Decraemer, W. F., Maes, M. A. and Vanhuyse, V. J. (1980) An elastic stress-strain relation for soft biological tissues based on a structural model. J. Biomechanics. To be published.

(Al)

7 =

cc

u(~,r)==d~G,(c)+

1 n=1.7...

z

in the last integral 0

[G,+inw

I 0

d

564

W.F.DECRAEMEII,M.A.MAES,V.J.VANHUYS~~~~~P.

The expression between the first pair of square brackets may be written as G:(W), where =

Wi,][e b,bq, - QWI

[ri- WJ12 _ b&h 2 %.

r= G:(o)

VANPEPERSTMIZE

J G&e-'"

G, + i

doz.

L

(A9)

’

0

+

k, - Q(41’] 1

For steady state conditions the expression for u reduced to

O’I!

is tiima,.

(B3)

J

To minimize X we used alternatively a grid search and a method of steepest descent.

0

NOMEI’KLATURE

‘43

APi'EM)IX B

The probability, per unit of r and of q, to find an experimental set of values (r,,qJ is given by

; F’

G(t) G,(t)

1 Pi =

(Bl) 2x0+&

- P:

~~,(r, - WdlCq-

-- 1 exp 2(1 - p:,

;:

H(f)

QWI I,

up,cl,

;

sample cross-section at rest = @W/a* force exerted on specimen instantaneous force response relaxation function reduced relaxation function storage modulus loss modulus relaxation spectrum fibre rest-length freqt==y (W probability of finding a set of M experimental value3 (r, = (01/@os4tr q1 - (o&&Wt) complex moduks intensity of the relaxation spectrum If(r) for r,CrGt,

where[Jhh QWI

ami (q, uq,)wremt

themeamand the

standard deviations of the marginal distributiona, and pi the correlation coe&isnt. The total probability of &ding a set of M values (rkq,) is given by

Ro s

t n

p = fj P,.

S32) 6

(=I

According to the principle of maximum likelihood, the best estimates for the parameters a, c. it, 72 are the values for which P ismaximal, i.e. when

d 4; P

intensity of the relaxation spsctrum &(7(r) for 51 < r < 7’2 standard deviation of l,distribution time (s) pipnrwta to be minimixed in order to obtain best parameter values 9 c, rt, r2 Straia

for harmonica1 e variation we note 8 = & + sd eio” phase angle of nth Fourier term for u phase angle of nth Fourier term for a’ mean of I,-diitribution instantaneous strees rapome Fourierc~ciertt of order n for o Fouriercdfficient of order n for 6