Periodic Neumann Boundary Conditions (2D Channel)

General Discussion
1

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,

2

  1. just a first remark, the characteristics method in not correctly implemented with periodic probem ( no way to day).

  2. no problem with Periodic Neumann become the function u will be peridic so dx(u)
    is also periodic.

  3. the characteristic method need a CFL close too 1 , a here the velocity is proportional to RE.

  4. p is defined through a constant so int2d(Th)(p) == 0 to set this constante.

Some correction but you have to change the method for do the Non linear part.

perio-NS.edp (1.4 KB)

1 Like
3

Thank you very much for your answer! It is very clear now. I will see how to implement the convective part.

4

you can just do for the term u.\nabla u the previous u like up.\nabla u ,
in this cas you mush rebuild the matrix in the loop so init = 1

1 Like