Domain and Pressure

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 div(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))
+x
psiT
dx(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)

+ 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))
+x
psiT
dx(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)

+ 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?

I haven’t looked at everything, but part of the issue must be that your operators are wrong. Cylindrical coordinates changes things (eg weak form of laplacian)

I do not understand, Can you explain a bit?