TD 10 Hydrostatique et diffusion de la chaleur

.

Hydrostatique

On considère un fluide isotherme en équilibre soumis à un champ de pesanteur uniforme. Ecrire l'équation de l'hydrostatique. On considère que le gaz est parfait c'est à dire que P(z) = k[B]*T*rho(z) . Résoudre simultanément l'ensemble des équations (hydrostatique et équation d'etat). Sachant que pour z=0, la pression est ; Tracer rho(z)/rho(0) et P(z)/P(0) en fonction de z pour rho(0) =0.7, m*g = 1 et k[B]*T = 1.

> ehy:=diff(P(z),z)+m*g*rho(z);

ehy := diff(P(z),z)+m*g*rho(z)

> etat:=P(z)=k[B]*T*rho(z);

etat := P(z) = k[B]*T*rho(z)

> resh:=dsolve({ehy,etat,rho(0)=0.7},{rho(z),P(z)});

resh := {P(z) = 7/10*k[B]*T*exp(-m/T*g/k[B]*z), rho...

> rgp:=[subs(m=k[B]*T/g,rhs(resh[1])/(k[B]*T)),subs(m=k[B]*T/g,P[0]=0.7,rhs(resh[2]))];

rgp := [7/10*exp(-z), 7/10*exp(-z)]

> plot(rgp[1],z=0..5);

[Maple Plot]

On considère maintenant que les particules ont une taille finie. L' équation d'état exacte à une dimension est donnée par P(z) = k[B]*T*rho(z)/(1-rho(z)) . Reprendre les questions précédentes.

> ehy:=diff(P(z),z)+m*g*rho(z);

ehy := diff(P(z),z)+m*g*rho(z)

> etat2:=P(z)=k[B]*T*rho(z)/(1-rho(z));

etat2 := P(z) = k[B]*T*rho(z)/(1-rho(z))

> resh2:=dsolve({ehy,etat2,rho(0)=0.7},{rho(z),P(z)});

resh2 := {rho(z) = k[B]*T*LambertW(1/k[B]/T*exp(-m/...
resh2 := {rho(z) = k[B]*T*LambertW(1/k[B]/T*exp(-m/...

> rr:=subs(m=k[B]*T/g,k[B]=1/T,rhs(resh2[2]));

rr := LambertW(exp(-z+ln(7/3)+7/3))

> pp:=subs(m=k[B]*T/g,k[B]=1/T,rhs(resh2[1]));

>

pp := LambertW(exp(-z+ln(7/3)+7/3))/(LambertW(exp(-...

> plot([pp,rr],z=0..5);

[Maple Plot]

Tracer les courbes rho(z) pour le gaz parfait et le gaz réel sur un même graphe. De même pour P(z). Conclusion

> plot({rgp[1],rr},z=0..5);

[Maple Plot]

> plot([rgp[2],pp],z=0..5);

[Maple Plot]

Diffusion de la chaleur

On considre un milieu homogène la température est theta[0] au temps t=0 compris entre x=a et x=b dont la diffusivité de chaleur est lambda . Ce mileu est mis en contact avec un thermostat à la temprature theta[1] au point x=a. On désigne par theta(x,t) la tempraure l'abscisse x et l'instant t . En l'absence de source de chaleur, le champ theta(x,t) obéit à l'équation de la chaleur : lambda*diff(theta(x,t),x,x) = diff(theta(x,t),t) .

Ecrire une procédure (ou une fonction) qui à une fonctionapplique l'équation différentielle précédente

> eq:=theta->lambda*diff(theta(x,t),x,x) = diff(theta(x,t),t):

Mur semi-infini. On prend ici a=0 et b = infinity .Recherche des solutions de la forme f(u) avec u = x^alpha*t^beta . Quels sont les couples ( alpha, beta ) possibles. Insérer la fonction f(x^alpha*t^beta) dans l'équation différentielle, puis faire apparaître u par subs(x^alpha = u/(t^beta),equa_diff) . Choisir ( alpha, beta ) pour que le résultat ne dépende que de u

> f:='f':F:=(x,t)->f(x^alpha*t^beta);
eq(F);
expand(subs(x^alpha = u/(t^beta), eq(F)*x^2/(u*D(f)(u))));

F := proc (x, t) options operator, arrow; f(x^alpha...

lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...
lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...

u/D(f)(u)*lambda*`@@`(D,2)(f)(u)*alpha^2+lambda*alp...

F := proc (x, t) options operator, arrow; f(x^alpha...

lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...
lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...

u*lambda*`@@`(D,2)(f)(u)*alpha^2/D(f)(u)+lambda*alp...

Déterminer l'équation différentielle pour f(u), lorsque beta = -1/2 . Résoudre pour les conditions initiales données dans l'introduction

> assume(t>0,u>0,x>0,lambda>0):
eq(F);equ:=simplify(subs(x = u*sqrt(t),beta=-1/2,alpha=1, eq(F)*t));

lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...
lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...
lambda*(`@@`(D,2)(f)(x^alpha*t^beta)*(x^alpha)^2*al...

equ := lambda*`@@`(D,2)(f)(u) = -1/2*D(f)(u)*u

> f1:=subs(alpha=1,beta=-1/2,f);

>

f1 := f

> sol:=collect(dsolve({equ,f(0)=Theta[1],f1(infinity)=Theta[0]},f(u)),erf);

sol := f(u) = (-Theta[1]+Theta[0])*erf(1/2/lambda^(...

Tracez theta(x,t) pour Theta[0] = 0, Theta[1] = 100 , lambda = L^2/tau en posant x = X*L et t/tau = .1, 1, 10, 100, 1000

> assume(L>0,tau>0);f2:=simplify(subs(u=x/sqrt(t),Theta[0] = 0, Theta[1] = 100, x = X*L, lambda=L^2/tau , t=mu*tau,rhs(sol)));

f2 := -100*erf(1/2*X/mu^(1/2))+100

> plot([seq(f2,mu=( .1,1,10,100,1000))],X=0..10);

[Maple Plot]

Diffusion de la chaleur dans un mur d'épaisseur a compris entre x=0 et x=a .On recherche ici les solutions particulières de la forme theta(x,t) = X(x)*T(t)

Montrer que theta[n](x,t) = b[n]*sin(n*Pi*x/a)*exp(-t/tau[n]) est une solution particulière pour tau[n] = tau[1]/(n^2) , et donner tau[1] . Pourquoi n doit-il être entier?

> fn:=(x,t)->b[n]*sin(n*Pi*x/a)*exp(-t/tau[n]);
eqd:=eq(fn);eqtau:=rhs(eqd)/lhs(eqd);
Tau[1]:=solve(subs(tau[n]=tau[1]/n^2,eqtau)=1,tau[1]);

>

fn := proc (x, t) options operator, arrow; b[n]*sin...

eqd := -lambda*b[n]*sin(n*Pi*x/a)*n^2*Pi^2/a^2*exp(...

eqtau := 1/tau[n]/lambda/n^2/Pi^2*a^2

Tau[1] := a^2/lambda/Pi^2

On recherche maintenant la solution du problème sous la forme theta(x,t) = sum(theta[n](x,t),n = 0 .. infinity)

Exprimer theta(x,0) en fonction de theta[n](x,0)

> theta(x,0)=sum(fn(x,0),n=0..infinity);

theta(x,0) = sum(b[n]*sin(n*Pi*x/a),n = 0 .. infini...

Calculer les b[n] sachant que theta(x,0) = Theta[0] pour x compris entre 0 et a et Theta[1] en dehors de cet intervalle

> assume(m,integer):assume(n,integer);
eqb:=int((Theta[0]-Theta[1])*sin(n*Pi*x/a),x=0..a)=int(b[n]*sin(n*Pi*x/a)*sin(n*Pi*x/a),x=0..a);
B:=unapply(solve(eqb,b[n]),n);

eqb := a*(Theta[1]-Theta[0])*((-1)^n-1)/n/Pi = 1/2*...

B := proc (n) options operator, arrow; 2*(Theta[1]*...

> simplify([B(2*n),B(2*n+1)]);

[0, -4*(Theta[1]-Theta[0])/(2*n+1)/Pi]

Tracer la fonction theta(x,0) lorsque la somme est tronquée à 10 termes, 20 termes, 100 termes

> i:='i';theta[approx]:=unapply(sum(B(i)*sin(i*Pi*x/a),i=1..n),x,n);

i := 'i'

theta[approx] := proc (x, n) options operator, arro...

> plot([seq(simplify(theta[approx](a*x,k)/(Theta[0]-Theta[1])),k=(10,20,100))],x=0..1,color=[red,green,blue]);

[Maple Plot]

Définir la fonction theta(x,t) pour une somme tronquée à 10 termes. Utiliser les variables sans dimension X=ax et T = t/tau[1]

> theta[approxt]:=(X,T,n)->(sum(simplify(B(i)/(Theta[0]-Theta[1]))*sin(i*Pi*X)*exp(-T*i^2),i=1..n));

theta[approxt] := proc (X, T, n) options operator, ...

> theta[approxt](x,t,10);

4/Pi*sin(Pi*x)*exp(-t)+4/3/Pi*sin(3*Pi*x)*exp(-9*t)...
4/Pi*sin(Pi*x)*exp(-t)+4/3/Pi*sin(3*Pi*x)*exp(-9*t)...

Tracer cette fonction pour t/tau[1] = 10^k pour k entier de -4 à 0. Que se passe-t-il lorsque le nombre de termes passe 20 puis 100?

> plot([seq(theta[approxt](x,10^(-4),n),n = (10, 20, ...
plot([seq(theta[approxt](x,10^(-4),n),n = (10, 20, ...

[Maple Plot]

> plot([seq(theta[approxt](x,10^(-3),n),n = (10, 20, ...
plot([seq(theta[approxt](x,10^(-3),n),n = (10, 20, ...

[Maple Plot]

> plot([seq(theta[approxt](x,10^(-2),n),n = (10, 20, ...
plot([seq(theta[approxt](x,10^(-2),n),n = (10, 20, ...

[Maple Plot]

> plot([seq(theta[approxt](x,10^(-1),n),n = (10, 20, ...
plot([seq(theta[approxt](x,10^(-1),n),n = (10, 20, ...

[Maple Plot]

Comment varie la température en x = a/2 ? Tracer les variations.

> (theta[approxt](1/2,t,20));

4/Pi*exp(-t)-4/3*1/Pi*exp(-9*t)+4/5/Pi*exp(-25*t)-4...
4/Pi*exp(-t)-4/3*1/Pi*exp(-9*t)+4/5/Pi*exp(-25*t)-4...

> plot([seq(theta[approxt](1/2,t,n),n=(10,20,100))],t=0..0.1,color=[red,green,blue]);

[Maple Plot]

> plot([seq(theta[approxt](1/2,t,n),n=(10,20,100))],t=0..5,color=[red,green,blue]);

[Maple Plot]

>

>