# Stokes in Cubic Pipe

Hello everyone,

I have been using the basic features of FreeFem for years, but now I am facing a problem that I can’t seem to make sense of. Specifically, I’m struggling to implement boundary conditions for a 3D Stokes solver.

For context, I have a cube with inflow and outflow conditions on the left and right faces, and I want to impose no-slip and no-penetration conditions on the other boundaries.

I’ve found that the documentation on this topic is possibly erroneous, and I’m having difficulty specifying the conditions I need, particularly on the left and right faces (Note that it NEEDS to be the left and right faces and I need to use `buildlayers` because I plan to make the problem more complex later on). I have no problem specify condition on the top and bottom face using labels 1 and 2. The issues seems to be with specifying labelmid (where the documentation says the other surfaces are labeled).

Does have any insights or suggestions on how to specify these boundary conditions? I would greatly appreciate any help.

-FluidMan

``````load "msh3" load "medit"
int nn=5;
int Wallt = 11;	//Pipe wall label
int Wallb = 22;	//Pipe wall label
int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Res = 5;   //Resolution of Solver

// Define mesh boundaries
border C01(t=0, 1){x=0; y=t; label=Inlet;}
border C02(t=0, 1){x=1-t; y=0; label=Wallb;}
border C03(t=0, 1){x=1; y=-t+1; label=Outlet;}
border C04(t=0, 1){x=t; y=1; label=Wallt;}

mesh Th2= buildmesh(C01(-Res) + C02(-Res) + C03(-Res) + C04(-Res));

fespace Vh2(Th2,P2);  Vh2 ux,uz,p2;
int[int] rup=[0,2],  rdown=[0,1], rmid=[11,3,22,4,33,5,44,6];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
zbound=[zmin,zmax],  labelmid=rmid,
reffaceup = rup,     reffacelow = rdown);

medit("c10x10x10",Th,wait=1);

// FFCS: testing 3d plots
plot(Th);

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (1-x)*(x)*y*(1-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
- div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,u1=0,u2=0,u3=0) + on(1,u1=0,u2=0,u3=0) + on(LEFT FACE,u1=fup,u2=0,u3=0)+ on(RIGHT FACE,u1=fup,u2=0,u3=0)+ on(BACK FACE,u1=0,u2=0,u3=0)+ on(FRONT FACE,u1=0,u2=0,u3=0);

for(int i=1;i<10;i++)
{
real yy=i/10.;
ux= u1(x,yy,y);
uz= u3(x,yy,y);
p2= p(x,yy,y);
plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
}

``````

hmm, I’m not sure it does what you want, but maybe this will help.

You can `buildlayers` and then rotate with `movemesh`

so

``````load "msh3" load "medit"

real Lx = 5.0, Ly=Lx, Lz=2*Lx;
int NN=5;
int[int] Labels=[0,1,2,3];
mesh Th2 = square(NN,NN,[Lx*x,Ly*y],label=Labels);

int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Faces = 66;	//Pipe outlet label

int[int] rup=[0,Outlet],  rdown=[0,Inlet], rmid=[0,Faces,1,Faces,2,Faces,3,Faces];
real zmin=-10,zmax=10;
mesh3 Th=buildlayers(Th2,int(NN*Lz/Lx),
zbound=[-0.5*Lz,0.5*Lz],  labelmid=rmid,
reffaceup = rup,     reffacelow = rdown);

plot(Th,wait=1);

real alpha=0.5*pi;
Th = movemesh(Th,[ x*cos(alpha)-z*sin(alpha)  ,y, x*sin(alpha)+z*cos(alpha)]);
plot(Th, wait=1);

fespace Vh2(Th2,P2);  Vh2 ux,uz,p2;

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (Lx-x)*(x)*y*(Ly-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
- div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(Faces,u1=0,u2=0,u3=0) + on(Inlet,u1=fup,u2=0,u3=0)+ on(Outlet ,u1=fup,u2=0,u3=0);

for(int i=1;i<10;i++)
{
real yy=i/10.;
ux= u1(x,yy,y);
uz= u3(x,yy,y);
p2= p(x,yy,y);
plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
}

``````

does this do as you want?

Thank you for the help. I am trying to work toward prescribing general Neumann conditions on all surfaces so I need to be able to specify velocity on all 6 surfaces of the cube. This code is closer but it still seems like conditions are initiable only on the bottom and top faces (then you rotated them to be the left and right faces).

You can separate them also

``````load "msh3" load "medit"

real Lx = 5.0, Ly=Lx, Lz=2*Lx;
int NN=5;
int[int] Labels=[0,1,2,3];
mesh Th2 = square(NN,NN,[Lx*x,Ly*y],label=Labels);

int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Face1 = 65;	//Pipe outlet label
int Face2 = 66;	//Pipe outlet label
int Face3 = 67;	//Pipe outlet label
int Face4 = 68;	//Pipe outlet label

int[int] rup=[0,Outlet],  rdown=[0,Inlet], rmid=[0,Face1,1,Face2,2,Face3,3,Face3];
real zmin=-10,zmax=10;
mesh3 Th=buildlayers(Th2,int(NN*Lz/Lx),
zbound=[-0.5*Lz,0.5*Lz],  labelmid=rmid,
reffaceup = rup,     reffacelow = rdown);

plot(Th,wait=1);

real alpha=0.5*pi;
Th = movemesh(Th,[ x*cos(alpha)-z*sin(alpha)  ,y, x*sin(alpha)+z*cos(alpha)]);
plot(Th, wait=1);

fespace Vh2(Th2,P2);  Vh2 ux,uz,p2;

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (Lx-x)*(x)*y*(Ly-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
- div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(Face1,u1=0,u2=0,u3=0)
+ on(Face2,u1=0,u2=0,u3=0)
+ on(Face3,u1=0,u2=0,u3=0)
+ on(Face4,u1=0,u2=0,u3=0)
+ on(Inlet,u1=fup,u2=0,u3=0)+ on(Outlet ,u1=fup,u2=0,u3=0);

for(int i=1;i<10;i++)
{
real yy=i/10.;
ux= u1(x,yy,y);
uz= u3(x,yy,y);
p2= p(x,yy,y);
plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
}

``````