Application Center - Maplesoft
  Submit your work
Application Center Applications Nonlinear Viscoelastic Behaviour of Brain Tissue Preview

App Preview:

Nonlinear Viscoelastic Behaviour of Brain Tissue

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



 

November 2015

 

                   Nonlinear Viscoelastic Behaviour of  Brain Tissue

 

                                   Univ.-Prof. Dr.-Ing. habil. Josef  BETTEN

                                             RWTH Aachen University

                 Mathematical Models in Materials Science and Continuum Machanics      

                                                    Augustinerbach 4-20

                                              D-52056  A a c h e n ,  Germany

 

                                              <betten@mmw.rwth-aachen.de>

 

Abstract

 

In this worsheet the relaxation of  Brain Tissue has been calculated to experimental data by two nonlinear model functions:

a five parameters PRONY-Series and compared with a three parameters Sqrt(t)-Law due to BETTEN, Creep Mechanics, 3rd edtion, 2008 Springer-Verlag Berlin / Heidelberg.

   The experiments have been carried out by P. AMEDIEU  and published in his dissertation Contribution à la biomécaniqe des tissus mous intracrániens,  2004 Université de Picardie Jules Verne, Faculte de Medicine, Amiens Cedex France

           

Relaxation

 

> 

restart:

> 

with(Statistics):

> 

X:=vector([0,25,50,100,150,200,250]);

X := Vector[row](7, {(1) = 0, (2) = 25, (3) = 50, (4) = 100, (5) = 150, (6) = 200, (7) = 250})

(1)
> 

whattype(X);

symbol

(2)
> 

Y:=vector([1,0.57,0.49,0.426,0.398,0.383,0.374]);

Y := Vector[row](7, {(1) = 1, (2) = .57, (3) = .49, (4) = .426, (5) = .398, (6) = .383, (7) = .374})

(3)
> 

whattype(Y);

symbol

(4)
> 

DATA:=seq([X[i],Y[i]],i=1..7);

DATA := [0, 1], [25, .57], [50, .49], [100, .426], [150, .398], [200, .383], [250, .374]

(5)

 

 

      PRONY-Series:

 

> 

R(t):=evalf(NonlinearFit(A+B*exp(-C*t)+D*exp(-E*t),X,Y,t),4);

R(t) := .3674+.3937*exp(-0.9820e-1*t)+.2389*exp(-0.1388e-1*t)

(6)
> 

R(0):=evalf(subs(t=0,R(t)));        R(infinity):=simplify(subs(t=infinity,R(t)));

R(0) := 1.0000

R(infinity) := .3674

(7)

 

      Sqrt(t)-Law:

 

> 

r(t):=evalf(NonlinearFit(a+b*exp(-c*sqrt(t)),X,Y,t),4);

r(t) := .3540+.6461*exp(-.2196*t^(1/2))

(8)
> 

r(0):=simplify(subs(t=0,r(t)));     r(infinity):=simplify(subs(t=infinity,r(t)));

r(0) := 1.0001

r(infinity) := .3540

(9)
> 

alias(th=thickness,co=color):

> 

p[1]:=plot(R(t),t=0..250,co=black,th=2,axes=boxed):

> 

p[2]:=plot(r(t),t=0..250,0..1,co=black,th=3,                     title="Five-Parameter PRONY & Three-Parameter Sqrt(t)"):

> 

p[3]:=plot(0.354,t=0..250,co=black,linestyle=4):

> 

plots[display](seq(p[k],k=1..3));

> 

# L[2]-Distance-Norm between PRONY R(t) and Sqrt(t)-Law r(t):   

> 

L[2]:=sqrt((1/250)*Int((R(tau)-r(tau))^2,tau=0..250))= evalf(sqrt((1/250)*int((R(t)-r(t))^2,t=0..250)),4);                                                                                       

L[2] := (1/50)*(10^(1/2)*(Int((R(tau)-r(tau))^2, tau = 0 .. 250))^(1/2)) = 0.1623e-1

(10)

   Error-Norm for PRONY  R(t)  to Experimental DATA:

 

> 

with(linalg):

> 

for i from 1 to 7 do              v[i]:=evalf(subs(t=DATA[i][1],R(t))-DATA[i][2]) od:

> 

V:=vector([seq(v[i],i=1..7)]);

V := Vector[row](7, {(1) = 0., (2) = 0.601329e-4, (3) = -0.3492771e-3, (4) = 0.104461655e-2, (5) = -0.81364654e-3, (6) = -0.71959670e-3, (7) = 0.833858634e-3})

(11)
> 

L[2]:=(1/sqrt(number_of_points))*Norm(V,2)=            evalf((1/sqrt(7))*norm(V,2),4);

L[2] := Norm(V, 2)/number_of_points^(1/2) = 0.6649e-3

(12)

 

   Error-Norm for Sqrt(t)-Law  r(t)  to Experimental DATA:

 

> 

for i from 1 to 7 do              w[i]:=evalf(subs(t=DATA[i][1],r(t))-DATA[i][2]) od:

> 

W:=vector([seq(w[i],i=1..7)]);

W := Vector[row](7, {(1) = 0.1e-3, (2) = -0.5014264e-3, (3) = 0.7490916e-3, (4) = -0.12314622e-3, (5) = -0.12193427e-3, (6) = -0.5662576e-4, (7) = 0.6091550e-4})

(13)
> 

l[2]:=(1/sqrt(number_of_points))*Norm(W,2)=        evalf((1/sqrt(7))*norm(W,2),4);

l[2] := Norm(W, 2)/number_of_points^(1/2) = 0.3503e-3

(14)
> 

Q:=evalf(L[2]/l[2],4);

Q := Norm(V, 2)/Norm(W, 2) = 1.898

(15)
> 

 

The above error-norms  L[2], l[2], or  Q = L[2] / l[2] illustrate that the Sqrt(t)-Law with only three parameters is a better approximation as the PRONY-Series with five parameters to be determined. In this worksheet and in a lot of other examples of  stress relaxation or even in  creep problems BETTEN  has  analysed in more detail that the SQRT(t)-Law furnishes the best approximation.

 

 

  Stress and Structural Relaxation

 

Besides the stress relaxation another kind, namely the structural relaxation, is also important. This sort of  relaxation governs the time-dependent response of a fluid to a change in temperature. For instance, a liquid is held at temperature T[1] until the property p( t, T[1]) reaches its equilibrium value p(infinity, T[1]), then it is suddenly cooled to T[2]. The instantaneous change in p is proportional to the difference T[1] - T[2],  followed by relaxation toward the equilibrium value p(infinity, T[2]).

Based upon a lot of experiments on  glass , SCHERER, G. (1986), Relaxation in Glass and Composites, has shown that both the stress and structural relaxation in glass can be predicted by the relation

 

> 

restart:

 

> 

r(t):=exp(-(t/lambda)^b);

r(t) := exp(-(t/lambda)^b)

(16)

 

 

often called KOHLRAUSCH function, which is a modified form of  the MAXWELL  relaxation function with  b = 1. The exponent  b  for a vriety of glasses was found by SCHERER to be near the value of  b = 0.5, so that the assumption of the Sqrt(t)-Law, introduced by BETTEN, is justified. However, the relation  r(t) is valid only for stabilized glass, i.e., the glass is held at a given temperature until its properties do no longer change with time, then the load can be applied.

In unstabilized glasses the viscosity eta and other typical properties, e.g., the density, vary with time. Then, the relaxation function should be replaced by the formular

 

> 

restart:

> 

r(t):=exp(-(Int(G[0]/eta(tau),tau=0..t))^b);

r(t) := exp(-(Int(G[0]/eta(tau), tau = 0 .. t))^b)

(17)

 

where, in agreement with experimental results, the exponent b can again be assumed to  b = 0.5, as has been discussed by SCHERER in more detail.

 

>