Hello,
I try to solve the base flow of a hated cylinder (2D)(natural convection) and try to understand the outer domain and pressure relationship. However, the pressure values sometimes increase and decrease when the domain increases.
I use the Newton method for iteration,
the code is;
fespace Mh(Th,[P2,P2,P1,P2]);
fespace Nh(Th,[P2,P2,P1]);
fespace XXh(Th,P2);
fespace Xh(Th,P1);
Mh [ux,uy,p,T]=[0,0,0,0];
Mh [dux,duy, dp, dT];
Mh [psix,psiy,psip, psiT];
macro grad1(v) [dx(v),dy(v)] //
macro div(ux,uy) (dx(ux)+dy(uy)) //
macro grad2(ux,uy) [dx(ux),dy(uy)] //
//mathematical solution with weak form (psi functions)
varf dNS([dux,duy,dp,dT],[psix,psiy,psip,psiT])=
int2d(Th)(x^2psiT(uydy(dT)+duydy(T)+uxdx(dT)+duxdx(T))
+x^2*(dy(psiT)dy(dT)+dx(psiT)dx(dT))
+xpsiTdx(dT)) // heat equation (cylindrical and weak form)
+int2d(Th)( x^2*psip*div(dux,duy) + x*dux*psip) // continuity (cylindrical and weak form)
+int2d(Th)( x^2*psiy*(uy*dy(duy)+duy*dy(uy)
+ux*dx(duy)+dux*dx(uy))
+ x^2*psix*(uy*dy(dux)+duy*dy(ux)
+ux*dx(dux)+dux*dx(ux))) //u.grad(u)(convective term of N-S)
+int2d(Th)( x^2*Pr*(2*grad2(dux,duy)' *grad2(psix,psiy)+div(duy,dux)*div(psiy,psix))
+ x*Pr*(2*psix*dx(dux)+psiy*div(duy,dux))
+ 2*Pr*dux*psix
- dp*(x^2*div(psix,psiy) + 2*x*psix)
- x^2*1e-8*dp*psip ) // diffusion and pressure term of N-S
-int2d(Th)(x^2*Ra*Pr*psiy*dT) // bousinessq approximation for buoyancy
//boundary conditions
+ on(axis,dux=0-ux, duy=0-uy, dT=0-T)
+ on(outer,dux=0-ux, duy=0-uy, dT=0-T)
+ on(sphere1, dux=0-ux, duy=0-uy, dT=0-T)
+ on(sphere2, dux=0-ux, duy=0-uy, dT=1-T);
//Solution for main variables
varf NS([dux,duy,dp,dT],[psix,psiy,psip,psiT])=
int2d(Th)(x^2psiT(uydy(T)+uxdx(T))
+x^2*(dy(psiT)dy(T)+dx(psiT)dx(T))
+xpsiTdx(T)) // heat equation (cylindrical and weak form)
+ int2d(Th)( x^2*psip*div(ux,uy) + x*ux*psip) // continuity (cylindrical and weak form)
+int2d(Th)( x^2*psix*(uy*dy(ux)
+ux*dx(ux))
+ x^2*psiy*(uy*dy(uy)
+ux*dx(uy))) //u.grad(u)(convective term of N-S)
+int2d(Th)( Pr*(x^2*(2*grad2(ux,uy)' *grad2(psix,psiy)+div(uy,ux)*div(psiy,psix))
+ x*(2*psix*dx(ux)+psiy*div(uy,ux))
+ 2*ux*psix)
- p*(x^2*div(psix,psiy) + 2*x*psix )
- x^2*1e-8*p*psip ) // diffusion and pressure term of N-S
-int2d(Th)(x^2*Ra*Pr*psiy*T); //bousinessq approximation for buoyancy
//iteration part
real iterNUM=1; //number of iterations
real errNS=1; //Navier-Stokes equ error
real tolST=1e-8; //target value for iteration
real normux0 = 1; //initial velocity
real normuy0 = 1; //initial velocity
real normp0=1; //initial pressure
real normT0=1; //initial temperature
real normNS=1; //initial N-S
real deltaerr=1; //delta of error
real normp=1; //target for pressure
real normT=1; //target for temperature
real normux; //target for velocity
real normuy; //target for velocity
while (normNS > tolST)
{
cout << " Iteration number =" << iterNUM <<endl;
cout << " Solving Navier-Stokes Equations" <<endl;
cout << “======================================================” <<endl;
matrix A = dNS(Mh, Mh,solver=sparsesolver);//matrix solution
real[int] b1 = NS(0,Mh); //constant part coming from initial guess
real[int] b2 = dNS(0,Mh); //constant part coming from Dirichlet BC's
real[int] b = b1 - b2;
real[int] q = (Mh.ndof);
q = A^-1*b ;
dux[] =q;
ux[] = ux[]- dux[];
errNS= normNS;
//plot([ux,uy], wait= false, value=true, fill=true);
// updating the L2 norms
normux = int2d(Th)(x*((dux)*(dux) )); // L2 norm
normuy = int2d(Th)(x*( (duy)*(duy) )); // L2 norm
normp = int2d(Th)(x*(dp)*(dp)); // L2 norm
normT = int2d(Th)(x*(dT)*(dT)); // L2 norm
normNS = sqrt(normux+normp+normT+normuy);
deltaerr= normNS-errNS;
cout << " NS residual = " << normNS <<endl;
cout << " velocity error X= " << sqrt(normux) <<endl;
cout << " velocity error Y= " << sqrt(normuy) <<endl;
cout << " pressure error= " << sqrt(normp) <<endl;
cout << " temperature error= " << sqrt(normT) <<endl;
cout << " NS total error= " << errNS <<endl;
cout << "======================================================" <<endl;
cout << "======================================================" <<endl;
cout << " Iteration number= " << iterNUM <<endl;
cout << " Stokes residual= " << normNS <<endl;
cout << "======================================================" <<endl;
iterNUM=iterNUM+1;
}
Can somebody help so that I can understand the reason?