Stokes in Cubic Pipe

General Discussion
1

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 Grad(u) [dx(u),dy(u),dz(u)]// EOM
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]) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3) //)';
		       - 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);
  }
2

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 Grad(u) [dx(u),dy(u),dz(u)]// EOM
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]) = 
 int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3) //)';
   	       - 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?

3

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).

4

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 Grad(u) [dx(u),dy(u),dz(u)]// EOM
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]) = 
 int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) +  Grad(u2)'*Grad(v2) +  Grad(u3)'*Grad(v3) //)';
   	       - 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);
 }