Hi, I am new to FreeFem++, and I am trying to simulate the unsteady Navier-Stokes equations in a 2D channel [0,1]\times[0,1] with some periodic boundary conditions. I denote the tangential velocity by u_1, the normal velocity by u_2, and the pressure by Pr.

The boundary conditions are

u_1(0, y) = u_1 (1, y),

u_2(0, y) = u_2 (1, y),

\partial_xu_2(0, y) = \partial_xu_2(1, y),

Pr(0, y) = Pr(1, y) + a,

u_1(x, 0) = 0,

u_1(x, 1) = 0,

u_2(x, 0) = f(x,t),

u_2(x, 1) = 0,

where a is some nonnegative constant, and f(x,t) is a given function.

To deal with the periodic Dirichlet boundary conditions, I use the ‘periodic’ function in the definition of the element space. Since we cannot mix periodic and non-periodic elements, I define p := Pr - Pbar, where Pbar = -ax+1, so that p is periodic, and I use p in the variational formulation.

My big problem now is how to implement the periodic Neumann boundary condition? In my code, everything runs without problem but I did not specify anywhere the periodic Neumann boundary condition, so I guess the simulated solution is just wrong. How to specify periodic Neumann boundary conditions?

My code is below (with f=0 and a=0.5):

```
// Parameters
real dt = 0.5; // Time step
int meshSize = 50;
real L= 1;
real RE = 1000;
// Mesh boundary
border b1(s = 0, 1){x=s*L; y=0; label=1;}
border b2(s= 0, 1){x=L; y= s; label=2;}
border b3(s= 0, 1){x=L-s*L; y= 1; label=3;}
border b4(s= 0, 1){x=0; y= 1-s; label=4;}
func perio=[[4,y],[2,y]];
mesh Th = buildmesh(b1(meshSize) + b2(meshSize) + b3(meshSize) + b4(meshSize));
plot(Th, wait=true);
fespace Vh(Th, P1, periodic=perio);
Vh u1, v1, u2, v2, p, q, up1, up2, pp;
int i=0;
real alpha=1/dt;
real nu = 1/RE;
real epsi = 0.00001;
real a = 0.5;
problem NS (u1, u2, p, v1, v2, q, solver=Crout, init=i)
= int2d(Th)(
alpha*(u1*v1 + u2*v2)
+ nu * (
dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2)
)
- p*q*(epsi)
- p*dx(v1) - p*dy(v2)
- dx(u1)*q - dy(u2)*q
)
+ int2d(Th)(a*x*q*(epsi) + a*x*dx(v1) + a*x*dy(v2))
- int2d(Th)(1*q*(epsi)+1*dx(v1)+1*dy(v2))
+ int2d(Th)(
- alpha*convect([up1,up2],-dt,up1)*v1
- alpha*convect([up1,up2],-dt,up2)*v2)
+ on(1, 3, u1=0)
+ on(3, u2=0)
+ on(1, u2=0);
u1 = -cos(pi*x/4.0)*exp(-y*y/2.0);
u2 = 0;
p = -cos(pi*x/4.0)*exp(-y*y/2.0)+a*x-1;
for (i = 0; i <= 500000; i++){
// Update
pp = p;
up1 = u1;
up2 = u2;
//vort = dx(u2)-dy(u1);
// Solve
NS;
// Plot
plot([u1,u2], fill=true);
}
```

Thank you in advance for your answers,