Navier Stokes with Newton method

General Discussion
1

Dear all,

I am trying to solve the Navier Stokes problem (a channel flow with a parabolic inflow) using the Newton method. I have the following code (inspired from what I have found on the FreeFEM tutorial Newton Method for the Steady Navier-Stokes equations).

verbosity = 0;

{
string Createrep="mkdir -p ./SOL_VTU ";                         
exec(Createrep);
}


mesh Th = square(30,30); 

func int labelbord(){
  return (1*(x==0)+2*(x==1)+3*(y==0)+4*(y==1)+10);} //function to put a label on each boundary edge of the mesh;

Th = change(Th, flabel=labelbord());

fespace Xh(Th, P2);
Xh u1, u2;
Xh v1,v2;
Xh du1,du2;

fespace Mh(Th,P1);
Mh p;
Mh q;
Mh dp;


// Macro
macro Grad(u1,u2) [dx(u1), dy(u1), dx(u2),dy(u2)] //
macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)], [u1,u2]'*[dx(v2),dy(v2)]] //
macro div(u1,u2) (dx(u1) + dy(u2)) //

// Initialization
u1 = 0;
u2 = 0;
real nu=1; 
real err=1; 
real eps = 1e-6; 

for (int n = 0; n < 100; n++){
   
    solve Oseen ([du1, du2, dp], [v1, v2, q])  
    = int2d(Th)(nu * (Grad(du1,du2)' * Grad(v1,v2))+ UgradV(du1,du2, u1, u2)' * [v1,v2]
                    + UgradV(u1, u2,du1,du2)' * [v1,v2] - div(du1,du2) * q - div(v1,v2) * dp 
                    - 1e-8*dp*q)
    - int2d(Th) ( nu * (Grad(u1,u2)' * Grad(v1,v2)) + UgradV(u1,u2, u1, u2)' * [v1,v2] - div(u1,u2) * q - div(v1,v2) * p)
      +on(11, du1=y*(1-y), du2=0) + on(12,du1=0, du2=0) + on(13,du1=0, du2=0) + on(14, du1=0 , du2=0) ; //+ on(10, du1=0, du2=0);
    
    u1[]-=du1[];
    u2[]-=du2[];
    p[]-=dp[];

    if (n>1){
    real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
    err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
    }

    if(err < eps) break; //converge

    if( n>3 && err > 10.) break; //blowup

    cout << err << endl; 

    load "iovtk"
    int[int] fforder=[1, 1];
    savevtk("./SOL_VTU/SOL_REF.vtu", Th, [u1, u2], p, dataname = "velocity pressure", order=fforder);
  }

However, the results I obtain are false (the velocity is too high and I do not obtain a poiseuille as expected).

Could someone help me to understand what is wrong with this code ?

Thank you in advance,

Best regards,

Loïc,

2

This is obviously wrong.

    u1[]-=du1[];
    u2[]-=du2[];
    p[]-=dp[];
3

Why ?
and how to update the solution otherwise ?

Loïc,

4

I think this is the boundary condition that is wrong, it should be
du1=u1-y*(1-y)

5

Because u1[], u2[], and p[] are all the same. So you are basically updating the solution with the Newton step three times. My bad, you are not using a vectorial fespace, so I was wrong.
There is a steady-state NS solver available out-of-the-box at FreeFem-sources/examples/hpddm/navier-stokes-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub.

6

But yes, see the above answer: at convergence, your linear form will never evaluate to 0 because of the on, so it is not consistent.

7

Thank you ! Change the boundary conditions solves the problem.

I will also see FreeFem-sources/examples/hpddm/navier-stokes-2d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub

Best regards,

Loïc,