paramath

      
 FINITE DIFFERENCE METHOD AND
THE  LAME'S EQUATION  IN HEREDITARY SOLID MECHANICS .                                          
by Co.H Tran  & Phong . T . Ngo , University of Natural Sciences  , HCMC  Vietnam -
          -  coth123@math.com   ,  coth123@yahoo.com  &  ntphong_6@yahoo.com        
                                                        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  ..

 Copyright
Co.H Tran - Phong .T .Ngo   - FINITE DIFFERENCE METHOD AND
THE  LAME'S EQUATION  IN HEREDITARY SOLID MECHANICS .   Use  Maple 9.5  
 . All rights reserved.  Copying or transmitting of this material without the permission of the authors is not allowed .

A. THE  DISPLACEMENT DIFFERENTIAL EQUATION :

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 .

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  ,      ,  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  ,      ,  z  (  truc  z  huong doc theo ong tru . Các  thanh phan bien dang và ung suat tuong ung la     la cac ham theo r ,  t  .  

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

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 )

B. FINITE DIFFERENCE METHOD  :

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 )  

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 ) :      ;      .  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) )   :  

     ;                         

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 )      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     with         and  the expression of     ;   

C. NUMERICAL SOLUTION  :

>

restart:with(plots):with(PDETools):with(LinearAlgebra):  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; row[0]:=[rhs(bound_con[1]),seq(0,j=1..N-1)]; row[1]:=[coeff(lhs(e[1]),u[0,t]),coeff(lhs(e[1]),u[1,t]),coeff(lhs(e[1]),u[2,t]),seq(0,j=1..N-3)];  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; row[N]:=[seq(0,j=1..N-1),rhs(bound_con[2])]; row_matrix:=[row[0],seq(row[n],n=1..N-2),row[N]];A:=(row_matrix);        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

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.