Skip to content

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{eqnarray*} (\mathbf{u}\cdot\nabla) \mathbf{u}-\nu \Delta \mathbf{u}+\nabla p&=&0\\ \nabla\cdot \mathbf{u}&=&0 \end{eqnarray*}

where \nu is the viscosity of the fluid, \nabla = (\p_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

\begin{equation} \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 \end{equation}

The Newton Algorithm to solve nonlinear problem is

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

  1. choose u_0\in \R^n , ;
  2. 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{eqnarray*} 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{eqnarray*}

So the Newton algorithm become

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
// 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.

Fig. 1: Mesh and the velocity and pressure at Reynolds 200
NSNewtonTh
NSNewtonUP