Equation de Chaleur sur FreeFem++

Réalisée par Charaf Rezrazi et Ayoub Outaghat

Problème traité

Dans ce qui suit, nous allons formuler une résolution d'une équation mise sous la forme suivante : $$\begin{cases} \frac{\partial u(x,y,t)}{\partial t} - \alpha\triangle u(x,y,t) = f(x,y,t) & ((x,y),t) \in \omega\times ]0,T[ \\u(x,y,t)=0 & ((x,y),t)\in\delta\omega\times]0,T[ \\u(x,y,0)=u_{0}(x,y) & ((x,y),t)\in\omega \end{cases}$$

Ceci revient à une formulation variationnelle : $$\int_{\omega}\frac{u^{m-1}v}{\triangle t} -\int_{\omega}\frac{u^{m}v}{\triangle t} -\alpha\int_{\omega}\triangledown u \triangledown v = \int_{\omega}fv$$

Maillage

Solution exacte et second membre

A partir d'un script Maple, on reproduit l'équation pour une fonction donnée : $$u_{exacte}(x,y,t)=(x ^ 2 - 4) * (y ^ 2 - 4) * (t - 9)$$ Qui par la suite nous retourne la fonction second membre qu'on va utiliser dans notre équation: $$f(x,y,t)= x ^ 2 * y ^ 2 - 2 * t * x ^ 2 - 2 * t * y ^ 2 + 14 * x ^ 2 + 14 * y ^ 2 + 16 * t - 128$$ Pour la valeur initiale de la solution, nous allons visualiser la variation pour chaque cas.

Visualisation

Explication du code

//-------------------------
// Domaine et maillage
//-------------------------

real[int] edges = [-2,2,-2,2]; // Coins du domaine
int n = 50; // Nombre d'elements
mesh pave = square(n,n,[edges[0]+(edges[1]-edges[0])*x,edges[2]+(edges[3]-edges[2])*y]); // Maillage rectangulaire
int BORDURE = 10;
int[int] labels = [1,BORDURE,2,BORDURE,3,BORDURE,4,BORDURE] ;
pave = change(pave,label=labels);
savemesh(pave,"pave.msh");
plot(pave,ps="pave.png"); // Graphique et exportation
            

Dans cette partie, le domaine rectangulaire est défini, ainsi que une généralisation de la bordure, et le nombre d'éléments du maillage. On exporte par la suite le maillage.

//-------------------------
// Constantes et espace elements finis
//-------------------------
real t=0,T = 1.; // [0,T] Temps
int M = 50; // Nombre d'elements temps
real dt = T/M; // Pas de discretisation
real alpha = 1; // Conductivite thermique
int i = 0; // Index

fespace VH(pave,P2);
VH u,v,old,f,uex;

//-------------------------
// Condition initiales
//-------------------------
VH[int] uh(M+2); // Vecteur regroupant les valeurs de uh
uh[0] = exp(x)*exp(y); // Valeur initiale de uh a t=0
old = uh[0];
            

Ici on definit les constantes ainsi que la valeur de la solution à t=0.

//-------------------------
// Calcul et exportation
//-------------------------

// Formulation variationnelle du probleme
problem thermic(u,v)=
int2d(pave)(dt*alpha*(dx(u)*dx(v)+dy(u)*dy(v)))
+int2d(pave)(u*v)
-int2d(pave)(old*v)
+on(BORDURE,u=0);
            

La formulation variationnelle traduite en syntaxe FreeFem++

// Erreurs
real[int] errL2(M+1);
real[int] errL2P(M+1);
real[int] errH1(M+1);
// Calcul en temps discretise
for(i;i<=M;i++){
	t=i*dt; // Incrementation du temps
	f = x ^ 2 * y ^ 2 - 2 * t * x ^ 2 - 2 * t * y ^ 2 + 14 * x ^ 2 + 14 * y ^ 2 + 16 * t - 128; // Fonction second membre qui change en fonction de t
	uex = (x ^ 2 - 4) * (y ^ 2 - 4) * (t - 9); // Solution exacte evaluee avec Matlab auparavant
	thermic; // Résolution du probleme
	//plot(u,wait=0,fill=1,cmm="t="+t,ps="snapshots/"+fname+"_"+i+"_frame_t"+t+".png"); // Graphique et exportation
	old = u; // On remplace la valeur de U(m-1) par U(m)
	uh[i+1] = old; // (Facultatif)

	// Evaluation des erreurs
	errL2[i] = sqrt(int2d(pave)((u-uex)^2));
	errH1[i] = sqrt(int2d(pave)((dx(u)-dx(uex))^2+(dy(u)-dy(uex))^2));
	errL2P[i] = errL2[i]/sqrt(int2d(pave)((dx(u)+dy(u))^2));
}
            

Dans cette partie, on définit des tableaux qui vont contenir les valeurs des distances. Ainsi que la fonction second membre calculée avec Maple à partir d'une solution exacte aussi définie. Les deux fonctions ont été définies à l'intérieur de la boucle car elles sont dépendentes de t.

// Exportation des erreurs
// Entete de fichier en Markdown
{
	ofstream f("evaluated_log.txt");
	f << "| T Value | L2 Evaluated Difference | H1 Evaluated Difference | L2 Evaluated Percentage   |" << endl;
	f << "|---------|-------------------------|-------------------------|---------------------------|" << endl;
};
for(int j=0;j<=M;j++){
	{
   		ofstream f("evaluated_log.txt",append);
   		f << "|"+j*dt+"|"+errL2[j]+"|"+errH1[j]+"|"+errL2[j]/errL2P[j]+"|" <<endl;
 	};
}
            

On exporte les valeurs dans finchier texte sous forme de tableau en Markdown pour la facilité de sa mise en forme.

Tableau des erreurs

T Value	L2 Evaluated Difference	H1 Evaluated Difference	L2 Evaluated Percentage
0	317.952	356.388	29.0086
0.02	316.983	354.851	22.8965
0.04	316.029	353.615	20.0583
0.06	315.086	352.465	18.1924
0.08	314.153	351.36	16.7931
0.1	313.228	350.281	15.6692
0.12	312.31	349.223	14.7279
0.14	311.399	348.179	13.9172
0.16	310.495	347.149	13.2051
0.18	309.597	346.129	12.5704
0.2	308.705	345.119	11.9982
0.22	307.818	344.117	11.4776
0.24	306.936	343.124	11.0005
0.26	306.06	342.138	10.5605
0.28	305.189	341.159	10.1526
0.3	304.323	340.186	9.77279
0.32	303.461	339.219	9.41772
0.34	302.604	338.258	9.08462
0.36	301.751	337.302	8.77118
0.38	300.902	336.352	8.47541
0.4	300.058	335.407	8.19564
0.42	299.218	334.467	7.93042
0.44	298.381	333.531	7.67846
0.46	297.548	332.6	7.43868
0.48	296.719	331.673	7.2101
0.5	295.894	330.75	6.99185
0.52	295.072	329.832	6.78317
0.54	294.254	328.917	6.58339
0.56	293.438	328.007	6.39188
0.58	292.626	327.1	6.2081
0.6	291.818	326.196	6.03155
0.62	291.012	325.296	5.86178
0.64	290.209	324.4	5.69837
0.66	289.41	323.507	5.54096
0.68	288.613	322.617	5.38921
0.7	287.818	321.73	5.24279
0.72	287.027	320.846	5.10143
0.74	286.238	319.965	4.96486
0.76	285.452	319.087	4.83283
0.78	284.668	318.212	4.70511
0.8	283.887	317.34	4.58151
0.82	283.108	316.47	4.46182
0.84	282.331	315.603	4.34586
0.86	281.557	314.738	4.23347
0.88	280.785	313.876	4.12448
0.9	280.015	313.016	4.01876
0.92	279.247	312.159	3.91617
0.94	278.481	311.304	3.81658
0.96	277.718	310.451	3.71986
0.98	276.956	309.6	3.62591
1	276.196	308.752	3.53461