* B5

by Co.H Tran  & Phong . T . Ngo , University of Natural Sciences  , HCMC  Vietnam -
          -   ,  &
                                                        Copyright  2005
                                                    Sat , May 15  2005  
 **  Abstract  :  The Lame's  differential equation  is solved by  the finite-difference  method .
 ** Subjects:  Viscoelasticity , Hereditary Solid Mechanics , The Differential equation  .

NOTE: This worksheet demonstrates the use of Maple for calculating the solution of Lame's  differential equation .

The authors expect that this worksheet will only be used for teaching and educational purposes  ..

Co.H Tran - Phong .T .Ngo   -
Use  Maple 9.5  
 . All rights reserved.  Copying or transmitting of this material without the permission of the authors is not allowed .



The Lame's  equation of the plane-deformation  problem in the cylinder made of orthotropic viscoelastic composite material does not have constant modules .
The modules Et , Er and Ert will be replaced with  the functions Er(t) , Et(t)  and Ert(t) respectively .

                                                                                             [Maple OLE 2.0 Object]

a*.   The plane-deformation problem  :  Bai toan bien dang phang  cua ong tru truc huong composite dan nhot :
We examine an orthotropic viscoelastic composite material cylinder which has the horizontal section within limit of 2 circles : r  =  a  , r  =  b  (  a  <  b  )  . Choosing the cylindrical  coordinates  r  ,   
[Maple OLE 2.0 Object]  ,  z  (  the axial  z  is along with the cylinder ) . The components of stress and  deformation  are functions of r , t   respectively .

Xet ong tru có tiet dien ngang gioi han boi 2 duong tron dong tam co ban kinh  r  =  a ,  r  = b   (  a  <  b  )  ,  ong tru duoc làm bang vat lieu có tính truc huong . Chon he toa do tru  r  ,    [Maple OLE 2.0 Object]  ,  z  (  truc  z  huong doc theo ong tru . Các  thanh phan bien dang và ung suat tuong ung la [Maple OLE 2.0 Object]    la cac ham theo r ,  t  .  

    The two components of deformation-tensor : ( 2 thanh phan cua tensor bien dang la : )    [Maple OLE 2.0 Object]
    and  the differential equation of equilibrium is : ( phuong trinh vi phan can bang )      
[Maple OLE 2.0 Object]    when t = 0  , boundary conditions  : ( khi t = 0 , cac dieu kien bien : )   [Maple OLE 2.0 Object]

b* .   The displacement - differential equation :    Phuong trinh vi phan chuyen vi :

The differential equation of the  cylinder displacement in the case of viscoelastic plane-deformation : ( Phuong trinh chuyen vi ong trong truong hop bien dang phang dan nhot ) [Maple OLE 2.0 Object]

Here   [Maple OLE 2.0 Object]  is   (t/To)^(2/5)/(100.+(t/To)^(1/2))
( To  :  const  )   


The boundary conditions  of the problem are given at two edges  ( Dieu kien bien cua bai toan duoc cho o 2 canh )  : r  =  a  and  r  =  b  .  ( a = 1 , b = 2 )    %? 

Now we choose the number of mesh points ( Ta chon so diem luoi ) N = 20  .  The interval over which we approximate this equation is ( Doan xap xi cua phuong trinh la ) [a, b] .  And the step size for this interval is ( Va kich thuoc buoc nhay cho doan nay la )   h := 1/20

The difference operators are ( Cac toan tu sai phan la ) Uj  and  Ujj  ,  And we have two boundary conditions  equations ( Va ta co 2 phuong trinh dieu kien bien ) :   e[0] := u[0,t] = 1   ;    e[20] := u[20,t] = -1  .  For determining the values at the interior mesh points  we obtain the N-1 equations ( De xac dinh cac gia tri cho cac diem trong , ta thu duoc N - 1 phuong trinh )  , then  by replacing   u '(x) and u ''(x) ( Va thay the u'(x) va u"(x) )   :  

U1 := proc (k, t) options operator, arrow; 1/2*(u[k+1,t]-u[k-1,t])/h end proc     ;    U2 := proc (k, t) options operator, arrow; (u[k+1,t]-2*u[k,t]+u[k-1,t])/(h^2) end proc                     

We arrange this system of N+1 equations in the form of matrix equation ( Sap xep he thong gom N+1 phuong trinh nay ) . The matrix of it has  N+1 rows( Ma tran chinh co N+1 hang ). The first row is fixed  with the boundary condition at r = a  ( Hang dau duoc xep cho dieu kien bien tai r = a ) . Obviously the last row is fixed with the boundary condition at r = b   ( Hien nhien hang cuoi cung duoc xep cho dieu kien bien tai  r = b ) .  Now, we join these rows  by listing them out , then construct the matrix symboled A . ( Lien ket cac hang nay lai , va xay dung nen ma tran A ) .

The unknown values  will be written as a vector ( cac gia tri chua biet se duoc viet dang vector )   u[j], j = 1 .. N   and the right hand side of the equations is a column vector  B ( va ve phai phuong trinh la 1 vector cot B ) . Solving the matrix equation  for u ( Giai phuong trinh ma tran tim nghiem u ) .  Then we find epsilon[theta](r,t) := u(r,t)/r    with      E[theta](r,t) := (100/((t/To)^.1)+1)*Ee   and  the expression of   sigma[theta](r,t) := u(r,t)*E[theta](r,t)/r  ;   


Use  Maple 9.5

 m:=(100.+(t/To)^(1/10))*(t/To)^(2/5)/(100.+(t/To)^(1/2)); To:=1;
lame_cyl:=r^2*diff(u(r,t),r$2)+r*diff(u(r,t),r)+m*u(r,t)=0; bound_con:=u(1,t)=1,u(2,t)=-1; a:=1; b:=2;
N:=20; h:=(b-a)/N;R:=k->a+k*h; U1:=(k,t)->(u[k+1,t]-u[k-1,t])/(2*h); U2:=(k,t)->(u[k+1,t]-2*u[k,t]+u[k-1,t])/h^2;
e[0]:=u[0,t]=rhs(bound_con[1]);e[N]:=u[N,t]= rhs(bound_con[2]);
for k from 1 to N-1 do e[k]:=eval(lame_cyl,{r=R(k),u(r,t)=u[k,t], diff(u(r,t),r)=U1(k,t),diff(u(r,t),r$2)=U2(k,t)} );od;
 for n from 2 to N-1 do row[n]:=[seq(0,j=1..n-2),coeff(lhs(e[n]),u[n-1,t]),coeff(lhs(e[n]),u[n,t]),coeff(lhs(e[n]),u[n+1,t]),seq(0,j=1..N-1-n)];od;
 U:=Vector([seq(u[j,t],j=1..N)]);B:=Vector([rhs(bound_con[1]),seq(rhs(e[j]),j=1..N-2),rhs(bound_con[2])]);                                                                 with(linalg):A.U = B;;;A1:=inverse(A);U:=linsolve(A,B); result:=[seq(evalf(subs(t=10,U(k,t)),1),k=0..N)]:
print("Ham epsilon[theta](1,t)");plot3d(-U[N-1]/r,r=1..1.000001,t=0..100 );print("Ham epsilon[theta](2,t)");plot3d(U[N-1]/r,r=2..2.000001,t=0..100);
u(r,t):=U[N-1];;Ee:=0.5;;epsilon[theta](r,t):=u(r,t)/r;E[theta](r,t) := (100/(t/To)^.1+1)*Ee; sigma[theta](r,t):=epsilon[theta](r,t)*E[theta](r,t) ;Ee:=0.5;sigma[a](t):=normal(subs(r=1,-sigma[theta](r,t)));sigma[b](t):=normal(subs(r=2,sigma[theta](r,t)));;;with(plottools):with(plots):plot([sigma[a](t),sigma[a](t)],t=0..10^10,style=[point,line],symbol=cross,color=[blue,black],thickness=[1],legend=[`sigma[an](t)`,`sigma[a](t)`],axes=box,symbolsize=[10]);plot([sigma[b](t),sigma[b](t)],t=0..10^10,style=[point,line],symbol=diamond,color=[red,black],thickness=[1],legend=[`sigma[bn](t)`,`sigma[b](t)`],axes=box,symbolsize=[10]);print("HAM u(r,t)");                                                                                                              smartplot3d(u(r,t));

Warning, the name changecoords has been redefined

m := (100.+(t/To)^(1/10))*(t/To)^(2/5)/(100.+(t/To)^(1/2))
To := 1
lame_cyl := r^2*diff(u(r,t),`$`(r,2))+r*diff(u(r,t),r)+(100.+t^(1/10))*t^(2/5)/(100.+t^(1/2))*u(r,t) = 0
bound_con := u(1,t) = 1, u(2,t) = -1
a := 1
b := 2
N := 20
h := 1/20
R := proc (k) options operator, arrow; a+k*h end proc
U1 := proc (k, t) options operator, arrow; 1/2*(u[k+1,t]-u[k-1,t])/h end proc
U2 := proc (k, t) options operator, arrow; (u[k+1,t]-2*u[k,t]+u[k-1,t])/h^2 end proc
e[0] := u[0,t] = 1
e[20] := u[20,t] = -1
e[1] := 903/2*u[2,t]-882*u[1,t]+861/2*u[0,t]+(100.+t^(1/10))*t^(2/5)/(100.+t^(1/2))*u[1,t] = 0
e[2] := 495*u[3,t]-968*u[2,t]+473*u[1,t]+(100.+t^(1/10))*t^(2/5)/(100.+t^(1/2))*u[2,t] = 0…..
[Maple Plot]
[Maple Plot]

[Maple Plot]
[Maple Plot]
[Maple Plot]

Disclaimer:  While every effort has been made to validate the solutions in this worksheet, the authors  are  not responsible for any errors contained and are not liable for any damages resulting from the use of this material.

Legal Notice : The copyright for this application is owned by the authors. Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the authors for permission if you wish to use this application in for-profit activities.

Today, there have been 33 visitors (70 hits) on this page!
=> Do you also want a homepage for free? Then click here! <=