====================================================================

  programmes Maple  -  24/03/06

  Méthode des éléments finis : thermique

     Yves DEBARD

     Institut Universitaire de Technologie du MANS
     Département Génie Mécanique et Productique
     Avenue Olivier Messiaen
     72085 LE MANS CEDEX 9

====================================================================

lin_3n

# calcul des matrices élémentaires

restart:with(linalg):

# représentation de la géométrie et jacobien

x:=(1+xi)*L/2;J:=L/2;

# fonctions d'interpolation

N:=[xi*(-1+xi)/2,1-xi*xi,xi*(xi+1)/2];
dNxi:=[-1/2+xi,-2*xi,xi+1/2];

# matrice de conductivité

B:=scalarmul(dNxi,1/J);
k:=Matrix(3,3,(i,j)->int(B[i]*B[j]*lambda*A*J,xi=-1..1),shape=symmetric);

# matrice de capacité

m:=Matrix(3,3,(i,j)->int(N[i]*N[j]*A*rho*cP*J,xi=-1..1),shape=symmetric);

# vecteur flux dû à une source volumique d'intensité q

f:=vector(3,i->int(N[i]*A*q*J,xi=-1..1));

--------------------------------------------------

nonlin_1

# températures imposées et lambda(T)

# résolution par la méthode de substitution

restart:with(linalg):

substitution:=1:
conductivite_constante:=2:

methode:=substitution;

L:=0.1:A:=1:

# températures imposées
T1:=20:T3:=100:

# interpolation
assume(xi,real):
x:=(1+xi)*L/2:
J:=L/2:
N:=vector([xi*(xi-1)/2,1-xi*xi,xi*(xi+1)/2]):
B:=vector([(2*xi-1)/L,(-4*xi)/L,(2*xi+1)/L]):

lambda:=proc(T::vector)
  local Txi:
  Txi:=dotprod(T,N):
  20+(120-20)/100*(Txi):
end:

T:=vector([T1,(T1+T3)*0.5,T3]);

precision:=0.01:
max_iterations:=20:
for iter from 1 to max_iterations do
  K:=matrix(3,3,(i,j)->int(A*lambda(T)*B[i]*B[j]*J,xi=-1..1)):
  Residu:=-K[2,1]*T[1]-K[2,3]*T[3]-K[2,2]*T[2]:
  if iter=1 then Kt1:=K[2,2]:fi:
  if methode=substitution then Kt:=K[2,2] else Kt:=Kt1:fi:
  dT:=Residu/Kt:
  T[2]:=T[2]+dT;
  print("iteration = ",iter);
  print("Kt",Kt);print("résidu",Residu);print("dT = ",dT);print("T2 = ",T[2]);
  if abs(dT)<precision then break:fi:
  if iter=max_iterations then print("non convergence en ",iter,"iterations");fi:
od:

dT2:=T[2]-0.5*(T[1]+T[3]);
F:=multiply(K,T):flux1:=F[1];flux3:=F[3];
plot([x,dotprod(T,N),xi=-1..1],
labels=["x en m","T en °C"],title="Température dans le mur");
 

---------------------------------------

nonlin_2

# température imposée et rayonnement en milieu infini

restart:

newton_raphson:=1:
newton_raphson_modifiee:=2:

methode:=newton_raphson;

stefan:=5.67e-8;TCelsius:=273.15;

L:=0.3;
T1:=300;
Tinfini:=800;
lambda:=40;
emissivite:=0.6;

T2:=0;
precision:=0.01:
max_iterations:=20:
if methode=newton_raphson_modifiee 
  then Kt:=lambda/L+4*emissivite*stefan*T2^3;fi:
for iter from 1 to max_iterations do
  if methode=newton_raphson
    then Kt:=lambda/L+4*emissivite*stefan*T2^3:fi;
  residu:=emissivite*stefan*(Tinfini^4-T2^4)-lambda/L*(T2-T1):
  dT2:=residu/Kt:
  T2:=T2+dT2;
  print("iteration = ",iter);
  print("Kt = ",Kt);
  print("Résidu = ",residu);
  print("dT2 = ",dT2);
  print("T2 = ",T2);
  if abs(dT2)<precision then break:fi:
  if iter=max_iterations then print("non convergence en ",iter,"iterations");fi:
od:

------------------------------------

nonlin_3
 
# rayonnement + convection + lambda(T)

# résolution par une méthode mixte : substitution + Newton-Raphson

restart:with(linalg):

stefan:=5.67e-8:TCelsius:=273.15:
L:=0.5:A:=1:
# convection
h:=100:Tfluide:=20+TCelsius:
# radiation
emissivite:=0.8:Tinfini:=500+TCelsius:

# interpolation
assume(xi,real):
x:=(1+xi)*L/2:
J:=L/2:
N:=vector([xi*(xi-1)/2,1-xi*xi,xi*(xi+1)/2]):
B:=vector([(2*xi-1)/L,(-4*xi)/L,(2*xi+1)/L]):

TNodales_celsius:=proc(T::vector)
   vector(3,i->T[i]-TCelsius):
 end:

lambda:=proc(T::vector)
  local Txi:
  Txi:=dotprod(T,N):
  5+(20-5)/500*(Txi-TCelsius):
end:

F:=proc(T::vector)
  vector([A*emissivite*stefan*(Tinfini^4-T[1]^4),0,A*h*(Tfluide-T[3])]):
end:

T:=vector([Tinfini,(Tinfini+Tfluide)*0.5,Tfluide]):TNodales_celsius(T);

precision:=0.01:
max_iterations:=20:
for iter from 1 to max_iterations do
  K:=matrix(3,3,(i,j)->int(A*lambda(T)*B[i]*B[j]*J,xi=-1..1)):
  Residu:=matadd(F(T),multiply(K,T),1,-1):
  # matrice tangente
  K[1,1]:=K[1,1]+A*4*emissivite*stefan*T[1]^3:
  K[3,3]:=K[3,3]+A*h:
  dT:=linsolve(K,Residu):
  T:=matadd(T,dT);
  print("iteration = ",iter);print("dT = ",dT);
  print("T = ",TNodales_celsius(T));
  if norm(dT)<precision then break:fi:
  if iter=max_iterations then print("non convergence en ",iter,"iterations");fi:
od:

dT2:=T[2]-0.5*(T[1]+T[3]);
flux1:=F(T)[1];flux3:=F(T)[3];
plot([x,dotprod(TNodales_celsius(T),N),xi=-1..1],
labels=["x en m","T en °C"],title="Température dans le mur");