Équation de diffusion 1D avec Mathematica

1. Fonction de calcul

Ci-dessous la fonction de calcul (télécharger).

n : nombre de points
type0 : type de la condition limite en x=0, "neumann" ou "dirichlet"
lim0 : valeur de la condition limite en x=0
type1 : type de la condition limite en x=1
lim1 : valeur de la condition limite en x=1
coef : coefficients de diffusion sous la forme {{x1,D1},{x2,D2},...{1,Dk}}
vecS : sources, liste à n éléments
vecU : état initial, liste à n éléments
temps : temps initial
dt : pas de temps
tf : temps final

  
diffusion[n_,type0_,lim0_,type1_,lim1_,coef_,vecS_,vecU_,temps_,dt_,tf_]:=Module[
    {t,u,dx,matA,matB,vecC,vecD,a,nd,alpha,j1,j2,k,j},
    t = temps;
    u = vecU;
    dx=1/(n-1);
    matA = SparseArray[{},{n,n}];
    matB = SparseArray[{},{n,n}];
    vecC = vecS*dt;
    a = 2*dt/dx^2;
    nd = Length[coef];
    vecD = Array[0,{n,1}];
    j1=1;
    For[k=1,k<=nd,k++,
        j2 = Floor[coef[[k]][[1]]*n];
        d = coef[[k]][[2]];
        For[j=j1,j<=j2,j++,
            vecD[[j]] = d
        ];
        j1 = j2;
    ]; 
    For[j=2,j<=n-1,j++,
        alpha = a*vecD[[j]];
        matA[[j,j-1]] = -alpha/2; matB[[j,j-1]] = alpha/2;
        matA[[j,j]] = 1+alpha; matB[[j,j]] = 1-alpha;
        matA[[j,j+1]] = -alpha/2; matB[[j,j+1]] = alpha/2; 
    ];
    If[type0=="dirichlet",
        matA[[1,1]] = 1; matA[[1,2]] = 0; vecC[[1]] = lim0;
    ];
    If[type0=="neumann",
        matA[[1,1]] = -1; matA[[1,2]] = 1; vecC[[1]] = dx*lim0;
    ];
    If[type1=="dirichlet",
        matA[[n,n-1]] = 0; matA[[n,n]] = 1; vecC[[n]] = lim1;
    ];
    If[type1=="neumann",
        matA[[n,n-1]] = 1; matA[[n,n]] = -1; vecC[[n]] = dx*lim1;
    ];
    For[k=1,k<=nd-1,k++,
        (*frontiere entre deux coefficients de diffusion differents*)
        j = Floor[coef[[k]][[1]]*n]-1; 
        matA[[j,j-1]] = -vecD[[j]]; 
        matA[[j,j]] = vecD[[j]]+vecD[[j+1]]; 
        matA[[j,j+1]] = -vecD[[j+1]];
        matB[[j,j-1]] = 0; matB[[j,j]] = 0; matB[[j,j+1]] = 0; vecC[[j]] = 0;
    ];
    While[t<tf,
        u = LinearSolve[matA,matB.u+vecC];
        t = t+dt;
    ];
    Return[{u,t}];
]

2. Exemples

2.a. Diffusion thermique dans une plaque

On considère une plaque (perpendiculaire à l'axe x) de conductivité thermique uniforme, soumise en x=0 à une température constante U=1 et en x=1 à une température constante U=0. Il n'y a aucune source thermique dans la plaque. Initialement la température est nulle sur l'intervalle [0,1].

Get["equationDiffusion.m"];
n = 200;
u = Table[0,{n}];
s = Table[0,{n}];
coef = {{1,1}};
t = 0;
{u1,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u,t,0.0001,0.001];
{u2,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u1,t,0.001,0.01];
{u3,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u2,t,0.01,0.1];
{u4,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u3,t,0.1,1];

ListLinePlot[{u1,u2,u3,u4},PlotRange->{{0,n},{0,1}},AxesLabel->{"x","T"}]

plotA.pdf

2.b. Diffusion thermique dans un système à trois couches

On considère un système isolé formé de deux plaques initialement à deux températures différentes, mises en contact thermique par une troisième plaque mince de conductivité plus faible.

n = 500;
u = Table[0,{n}];
s = Table[0,{n}];
For[j=1,j<=Floor[n/2],j++,u[[j]]=1];
coef = {{0.45,1},{0.55,0.05},{1,1}};
t = 0;
{u1,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u,t,0.00001,0.001];
{u2,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u1,t,0.001,0.01];
{u3,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u2,t,0.01,0.1];
{u4,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u3,t,0.1,1];

ListLinePlot[{u1,u2,u3,u4},PlotRange->{{0,n},{0,1}},AxesLabel->{"x","T"}]

plotB.pdf