Newton Method for the Steady Navier-Stokes equations

The problem is find the velocity field \(\mathbf{u}=(u_i)_{i=1}^d\) and the pressure \(p\) of a Flow satisfying in the domain \(\Omega \subset \mathbb{R}^d (d=2,3)\):

\[\begin{split}\begin{array}{rcl} (\mathbf{u}\cdot\nabla) \mathbf{u}-\nu \Delta \mathbf{u}+\nabla p&=&0\\ \nabla\cdot \mathbf{u}&=&0 \end{array}\end{split}\]

where \(\nu\) is the viscosity of the fluid, \(\nabla = (\partial_i )_{i=1}^d\), the dot product is \(\cdot\), and \(\Delta = \nabla\cdot\nabla\) with the same boundary conditions (\(\mathbf{u}\) is given on \(\Gamma\)).

The weak form is find \(\mathbf{u}, p\) such that for \(\forall \mathbf{v}\) (zero on \(\Gamma\)), and \(\forall q\):

\[\int_\Omega ((\mathbf{u}\cdot\nabla) \mathbf{u} ). \mathbf{v} + \nu \nabla \mathbf{u}:\nabla \mathbf{v} - p \nabla\cdot \mathbf{v} - q \nabla\cdot \mathbf{u} = 0\]

The Newton Algorithm to solve nonlinear problem is:

Find \(u\in V\) such that \(F(u)=0\) where \(F : V \mapsto V\).

choose \(u_0\in \mathbb{R}^n\) , ;
for ( \(i =0\); \(i\) < niter; \(i = i+1\))
1. solve \(DF(u_i) w_i = F(u_i)\);
2. \(u_{i+1} = u_i - w_i\);

break \(|| w_i|| < \varepsilon\).

Where \(DF(u)\) is the differential of \(F\) at point \(u\), this is a linear application such that:

\[F(u+\delta) = F(u) + DF(u) \delta + o(\delta)\]

For Navier Stokes, \(F\) and \(DF\) are:

\[\begin{split}\begin{array}{rcl} F(\mathbf{u},p) &=& \int_\Omega((\mathbf{u}\cdot\nabla) \mathbf{u} ). \mathbf{v} + \nu \nabla \mathbf{u}:\nabla \mathbf{v} - p \nabla\cdot \mathbf{v} - q \nabla\cdot \mathbf{u}\\ DF(\mathbf{u},p)(\mathbf{\delta u} ,\delta p) &=& \int_\Omega ((\mathbf{\delta u}\cdot\nabla) \mathbf{u} ). \mathbf v + ((\mathbf{u}\cdot\nabla) \mathbf{\delta u} ). \mathbf{v} \\ &+& \nu \nabla \mathbf{\delta u}:\nabla \mathbf{v} - \delta p \nabla\cdot \mathbf{v} - q \nabla\cdot \mathbf{\delta u} \end{array}\end{split}\]

So the Newton algorithm become:

// Parameters
real R = 5.;
real L = 15.;

real nu = 1./50.;
real nufinal = 1/200.;
real cnu = 0.5;

real eps = 1e-6;

verbosity = 0;

// Mesh
border cc(t=0, 2*pi){x=cos(t)/2.; y=sin(t)/2.; label=1;}
border ce(t=pi/2, 3*pi/2){x=cos(t)*R; y=sin(t)*R; label=1;}
border beb(tt=0, 1){real t=tt^1.2; x=t*L; y=-R; label=1;}
border beu(tt=1, 0){real t=tt^1.2; x=t*L; y=R; label=1;}
border beo(t=-R, R){x=L; y=t; label=0;}
border bei(t=-R/4, R/4){x=L/2; y=t; label=0;}
mesh Th = buildmesh(cc(-50) + ce(30) + beb(20) + beu(20) + beo(10) + bei(10));
plot(Th);

//bounding box for the plot
func bb = [[-1,-2],[4,2]];

// Fespace
fespace Xh(Th, P2);
Xh u1, u2;
Xh v1,v2;
Xh du1,du2;
Xh u1p,u2p;

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

// 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 = (x^2+y^2) > 2;
u2 = 0;

// Viscosity loop
while(1){
    int n;
    real err=0;
    // Newton loop
    for (n = 0; n < 15; n++){
        // Newton
        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 //stabilization term
            )
            - 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(1, du1=0, du2=0)
            ;

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

        real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
        err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;

        cout << n << " err = " << err << " " << eps << " rey = " << 1./nu << endl;
        if(err < eps) break; //converge
        if( n>3 && err > 10.) break; //blowup
    }

    if(err < eps){  //converge: decrease $\nu$ (more difficult)
        // Plot
        plot([u1, u2], p, wait=1, cmm=" rey = " + 1./nu , coef=0.3, bb=bb);

        // Change nu
        if( nu == nufinal) break;
        if( n < 4) cnu = cnu^1.5; //fast converge => change faster
        nu = max(nufinal, nu* cnu); //new viscosity

        // Update
        u1p = u1;
        u2p = u2;
        pp = p;
    }
    else{   //blowup: increase $\nu$ (more simple)
        assert(cnu< 0.95); //the method finally blowup

        // Recover nu
        nu = nu/cnu;
        cnu= cnu^(1./1.5); //no conv. => change lower
        nu = nu* cnu; //new viscosity
        cout << " restart nu = " << nu << " Rey = " << 1./nu << " (cnu = " << cnu << " ) \n";

        // Recover a correct solution
        u1 = u1p;
        u2 = u2p;
        p = pp;
    }
}

Note

We use a trick to make continuation on the viscosity \(\nu\), because the Newton method blowup owe start with the final viscosity \(\nu\).

\(\nu\) is gradually increased to the desired value.

NSNewtonUP — Fig. 26 Velocity and pressure

Naver-Stokes newton

Here you can use the search bar to search in all FreeFEM examples instead

The list of all matching examples on the right is updated for each new query

Clicking on a search hit or on an example on the right brings up the corresponding script from GitHub

You can also filter the list on the right by selecting one or more intersecting categories on the left tree

Newton Method for the Steady Navier-Stokes equations