Problème de Darcy couplée avec la chaleur — Les-mathematiques.net The most powerful custom community solution in the world

Problème de Darcy couplée avec la chaleur

j'ai essayé de programmer le code Freefem de la première formulation variationnelle qui est le suivante. $$

(V_{h,1})\left\{\begin{aligned}
\text { Chercher }\left(\mathbf{u}_{h}, p_{h}, T_{h}\right) \in \mathcal{W}_{h, 1} \times M_{h, 1} \times X_{h} \text { tel que : } \quad \quad \quad \quad \quad \quad \\
\qquad \forall \mathbf{v}_{h} \in \mathcal{W}_{h, 1}, \quad \int_{\Omega} \nu\left(T_{h}(\mathbf{x})\right) \mathbf{u}_{h}(\mathbf{x}) \cdot \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}-\int_{\Omega} p_{h}(\mathbf{x}) \operatorname{div} \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}= \\
\int_{\Omega} \mathbf{f}(\mathbf{x}) \cdot \mathbf{v}_{h}(\mathbf{x}) d \mathbf{x}, \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \\
\forall \mathbf{q}_{h} \in M_{h, 1}\text{, } \quad \int_{\Omega} q_{h}(\mathbf{x}) \operatorname{div} \mathbf{u}_{h}(\mathbf{x}) d \mathbf{x}=0,
\quad \quad \quad \quad \quad \quad \quad \quad \quad \\
\forall S_{h} \in X_{h}, \alpha \int_{\Omega} \nabla T_{h}(\mathbf{x}) \cdot \nabla S_{h}(\mathbf{x}) d \mathbf{x}+\int_{\Omega}\left(\mathbf{u}_{h} \cdot \nabla T_{h}\right)(\mathbf{x}) S_{h}(\mathbf{x}) d \mathbf{x}= \\
\int_{\Omega} g(\mathbf{x}) S_{h}(\mathbf{x}) d \mathbf{x}. \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad
\end{aligned}
\right.

$$ On discrétisé la vitesse $v$ par l'élément de Raviart-Thomas (RT0), $P0$ pour la pression $p$ et $P1$ pour la température $T$.

J’ai eu plusieurs problèmes.
Le premier est que la solution approchée n'est pas valable pour notre problème, c-à-d le code en n'a pas encore validé et je ne sais pas est ce que si le code que j'utilise est valable ou non ?

Le deuxième est de calculer l'erreur, j'ai traité un code mais le compilateur n'accepte pas ce code et je ne sais pas où le problème ? $$
erreur=||u-u_h||_{L^2(\Omega)^d}+||div( u-u_h)||_{L^2(\Omega)}+||p-p_h||_{L^2(\Omega)}+|T-T_h|_{H^1(\Omega)}.

$$ Voici le code que j'ai traité et ces joints des fichiers pdf
Merci cordialement.
Yassine OUZROUR.
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro Lablace(T) (dxx(T))+dyy(T)))//
macro GradientNorme(T) (dx(T)*dx(T)+dy(T)*dy(T))//
//donners
int  d=2;
real gm=1;
real alpha=3.;
real bitha=5;
//solution exacte
func phi=exp(-2*((x-1)^2+(y-2)^2));
func Pexct=gm*cos(pi*x/3)*cos(pi*y/3);
func Texct=gm*x^2*(x-3)^2*(y-3)^2*y^2;
func UXexct=gm*2*bitha*(1-x)*exp(-bitha*((x-1)^2+(y-1)^2));
func UYexct=gm*2*bitha*(1-y)*exp(-bitha*((x-1)^2+(y-1)^2));
//les fonctions donners
func fx =-(pi/3)*sin((pi/3)*x)*cos((pi/3)*y)+(Texct+1)*UXexct;
func fy =-(pi/3)*cos((pi/3)*x)*sin((pi/3)*y)+(Texct+1)*UYexct;
func g=-alpha*(2*(x-3)^2*(y-3)^2*y^2+4*x*(x-3)*(y-3)^2*y^2
+2*x^2*(y-3)^2*y^2+4*x*(x-3)*(y-3)^2*y^2
+2*(x-3)^2*(y-3)^2*x^2+4*(x-3)^2*y*(y-3)*x^2
+4*x^2*(y-3)*(x-3)^2*y+2*(y)^2*(x-3)^2*x^2
+UXexct*(2*x^2*(x-3)*(y-3)^2*y^2+2*x*(x-3)^2*(y-3)^2*y^2)
+UYexct*(2*x^2*(x-3)^2*(y-3)*y^2+2*x^2*(x-3)^2*(y-3)^2*y));
//////////////////
func Uexct=[UXexct,UYexct];
//Mesh
real L2erreur;
real[int] L2error(2);
//for(int n=0 ;n<2 ;n++){
//mesh Th=square(40*(n+1),40*(n+1));
mesh Th = square(100,100,[3*x,3*y]);
//Fespace
fespace Uh(Th,RT0);//Raviart-Thomas
Uh [ux, uy];
Uh [vx, vy];
Uh [FEUXexct,FEUYexct]=Uexct;
fespace Ph(Th, P0);
Ph p;
Ph q;
Ph FEPexct=Pexct;
fespace Zh(Th, P1);
Zh T;
Zh S;
Zh FETexct=Texct;
//Function
func nu=T+1;
func uT=Gradient(T)'*[ux,uy];
//func GN=dx(FETexct-T)*dx(FETexct-T)+dy(FETexct-T)*dy(FETexct-T);
//Problèm
problem Vh1([ux, uy,T,p],[vx, vy,S,q])
       =int2d(Th)(nu*[ux,uy]'*[vx,vy])
       +int2d(Th)(-p*(dx(vx)+dy(vy)))
       -int2d(Th)( fx*vx + fy*vy )
       -int2d(Th)(q*(dx(ux)+dy(uy)))
       -int2d(Th)(alpha*[dx(T),dy(T)]'*[dx(S),dy(S)])
       -int2d(Th)(uT*S)        
       +int2d(Th)(g*S) +on(1,2,3,4,T=0)
;
L2erreur= (int2d(Th)((FEPexct-p)^2))^(0.5)
        +int2d(Th)(GradientNorme(FETexct-T)^2)^(0.5)
        +(int2d(Th)((FEUXexct-ux)^2+(FEUYexct-uy)^2))^(0.5)
        +(int2d(Th)(Divergence((UXexct-ux),(UYexct-uy))^2))^0.5
        ;
Vh1;
//if(L2erreur<0.9){
plot(Th,wait=true, value=true ,cmm="Maillage de domaine");
plot(p,wait=true, value=true , cmm="Pression");
plot(FEPexct,wait=true, value=true , cmm="Pression exacte");
plot([ux, uy],wait=true, value=true , cmm="Vitesse");
plot([FEUXexct,FEUYexct],wait=true, value=true , cmm="Vitesse exacte");
plot(T,wait=true, value=true ,cmm="cheleur");
plot(FETexct,wait=true, value=true ,cmm="cheleur exacte");
                  // }
//cout << " L2error " << " = "<< L2erreur<<endl;
//}
105522
Connectez-vous ou Inscrivez-vous pour répondre.
Success message!