Average condition for pressur in stokes equation

General Discussion
1

Hello everyone,i have a small question is it possible to add a condition that says the average of pressure that we want to calculate should equal to 0 ,here is my code that calculate the stokes equation in porous medium
"load “msh3”
load “iovtk”

real Mu = 1.;
real pEps = 1.e-10;
real r=0.1;

border C01(t=-r, r){x=-1; y=t; label=10;}
border C02(t=-1, -r){x=t; y=r; label=1;}
border C03(t=-1, -r){x=t; y=-r; label=2;}
border C04(t=r, 1){x=-r; y=t; label=3;}
border C05(t=-1, -r){x=-r; y=t; label=4;}
border C06(t=-r, r){x=t; y=1; label=11;}
border C07(t=-r, r){x=t; y=-1; label=22;}
border C08(t=-1, -r){x=r; y=t; label=5;}
border C09(t=r, 1){x=r; y=t; label=6;}
border C10(t=r, 1){x=t; y=r; label=7;}
border C11(t=r, 1){x=t; y=-r; label=8;}
border C12(t=-r, r){x=1; y=t; label=36;}

int n = 100;

mesh Th = buildmesh(C01(-n)+C02(-n)+C04(-n)+ C06(-n)
+C09(n)+C10(-n)+C12(n)+C11(n)+C08(n)+ C07(n)+C05(-n) + C03(n));

border C1(t=-1, 1){x=-1; y=t; label=44;}
border C2(t=-1, 1){x=t; y=1; label=45;}
border C3(t=1, -1){x=1; y=t; label=46;}
border C4(t=1, -1){x=t; y=-1; label=47;}

mesh Th1 = buildmesh(C1(-n) + C2(-n) + C3(-n) + C4(-n));

plot(Th, wait=true);

fespace All(Th, [P2, P2, P1], periodic=[[10,y],[36,y],[11,x],[22,x]]);
All [Ux, Uy, p];
All [Uhx, Uhy, ph];

// Macros remain identical
macro grad(u) [dx(u), dy(u)] //
macro Grad(U) [grad(U#x), grad(U#y)] //
macro div(ux, uy) (dx(ux) + dy(uy)) //
macro Div(U) div(U#x, U#y) //

// Modified problem definition
problem Stokes ([Ux, Uy, p], [Uhx, Uhy, ph])
= int2d(Th)(
Mu * (Grad(U) : Grad(Uh))
- p * Div(Uh)
- ph * Div(U)
- pEps * p * ph
)
+ int2d(Th)( 1 * Uhx + 0* Uhy)
+ on(1,2,3,4,5,6,7,8, Ux=0, Uy=0)
;

// Solve
Stokes;

// Plot
plot(Ux, fill=true, value=true, wait=true, cmm=“vitesse Ux”);
plot(Uy, fill=true, value=true, wait=true, cmm=“Vitesse Uy”);
plot(p, fill=true, value=true, wait=true, cmm=“Pression p”);

real area = int2d(Th1)(1.0);

real Kxx = -int2d(Th)(Ux) / area;
real Kyy = -int2d(Th)(Uy) / area;

cout << "Kxx = " << Kxx << endl;
cout << "Kyy = " << Kyy << endl;
real ara = int2d(Th)(1.0);
real pmoy = int2d(Th)(p) / ara;
cout << "Moyenne de p = " << pmoy << endl;

fespace Vh1(Th, P1);
Vh1 Ux1 = Ux, Uy1 = Uy, p1 = p;

int[int] order = [1, 1, 1]; 
savevtk("StokesSolution.vtu", Th, Ux1, Uy1, p1, 
        dataname="Velocity_X Velocity_Y Pressure", order=order);

" and thank you so much for your help

2

Hello,
Your formulation indeed sets the average pressure to zero.
To see that, take in the variational formulation Uhx=0, Uhy=0, ph=1.
It gives \int_\Omega \mathop{\rm div} u+\epsilon\int_\Omega p=0.
You have \int_\Omega \mathop{\rm div} u=\int_{\partial \Omega}u\cdot N=0 since u=0 on \partial\Omega, thus \int_\Omega p=0.

3

oh ok thank you so much, but when i calculate the average in the end it always give something li 10^-5 or something like but never 0,is it normal to have value near to 0 but never 0 ??