# Domain and Pressure

**URL:** https://community.freefem.org/t/domain-and-pressure/2879
**Category:** General Discussion
**Created:** [January 13, 2024, 1:11pm UTC](https://community.freefem.org/t/domain-and-pressure/2879 "2024-01-13T13:11:29Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![Efe\_Taner](https://avatars.discourse-cdn.com/v4/letter/e/8e7dd6/32.png) [@Efe\_Taner](https://community.freefem.org/u/Efe_Taner)
#### Post date: [January 13, 2024, 1:11pm UTC](https://community.freefem.org/t/domain-and-pressure/2879/1 "2024-01-13T13:11:29Z")

</div>

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^2_psiT_(uy_dy(dT)+duy_dy(T)+ux_dx(dT)+dux_dx(T))  
+x^2\*(dy(psiT)_dy(dT)+dx(psiT)dx(dT))  
+xpsiT_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)
	       	          +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^2_psiT_(uy_dy(T)+ux_dx(T))  
+x^2\*(dy(psiT)_dy(T)+dx(psiT)dx(T))  
+xpsiT_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)
	       	             +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?

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [January 13, 2024, 1:33pm UTC](https://community.freefem.org/t/domain-and-pressure/2879/2 "2024-01-13T13:33:08Z")

</div>

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)

---

<div class="post-metadata">

### Author: ![Efe\_Taner](https://avatars.discourse-cdn.com/v4/letter/e/8e7dd6/32.png) [@Efe\_Taner](https://community.freefem.org/u/Efe_Taner)
#### Post date: [January 13, 2024, 1:45pm UTC](https://community.freefem.org/t/domain-and-pressure/2879/3 "2024-01-13T13:45:21Z")

</div>

I do not understand, Can you explain a bit?