FreeFEM Documentation on GitHub

stars - forks

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\).

  1. choose \(u_0\in \mathbb{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{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:

  1 // Parameters
  2 real R = 5.;
  3 real L = 15.;
  4 
  5 real nu = 1./50.;
  6 real nufinal = 1/200.;
  7 real cnu = 0.5;
  8 
  9 real eps = 1e-6;
 10 
 11 verbosity = 0;
 12 
 13 // Mesh
 14 border cc(t=0, 2*pi){x=cos(t)/2.; y=sin(t)/2.; label=1;}
 15 border ce(t=pi/2, 3*pi/2){x=cos(t)*R; y=sin(t)*R; label=1;}
 16 border beb(tt=0, 1){real t=tt^1.2; x=t*L; y=-R; label=1;}
 17 border beu(tt=1, 0){real t=tt^1.2; x=t*L; y=R; label=1;}
 18 border beo(t=-R, R){x=L; y=t; label=0;}
 19 border bei(t=-R/4, R/4){x=L/2; y=t; label=0;}
 20 mesh Th = buildmesh(cc(-50) + ce(30) + beb(20) + beu(20) + beo(10) + bei(10));
 21 plot(Th);
 22 
 23 //bounding box for the plot
 24 func bb = [[-1,-2],[4,2]];
 25 
 26 // Fespace
 27 fespace Xh(Th, P2);
 28 Xh u1, u2;
 29 Xh v1,v2;
 30 Xh du1,du2;
 31 Xh u1p,u2p;
 32 
 33 fespace Mh(Th,P1);
 34 Mh p;
 35 Mh q;
 36 Mh dp;
 37 Mh pp;
 38 
 39 // Macro
 40 macro Grad(u1,u2) [dx(u1), dy(u1), dx(u2),dy(u2)] //
 41 macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],
 42                         [u1,u2]'*[dx(v2),dy(v2)]] //
 43 macro div(u1,u2) (dx(u1) + dy(u2)) //
 44 
 45 // Initialization
 46 u1 = (x^2+y^2) > 2;
 47 u2 = 0;
 48 
 49 // Viscosity loop
 50 while(1){
 51     int n;
 52     real err=0;
 53     // Newton loop
 54     for (n = 0; n < 15; n++){
 55         // Newton
 56         solve Oseen ([du1, du2, dp], [v1, v2, q])
 57             = int2d(Th)(
 58                     nu * (Grad(du1,du2)' * Grad(v1,v2))
 59                 + UgradV(du1,du2, u1, u2)' * [v1,v2]
 60                 + UgradV( u1, u2,du1,du2)' * [v1,v2]
 61                 - div(du1,du2) * q
 62                 - div(v1,v2) * dp
 63                 - 1e-8*dp*q //stabilization term
 64             )
 65             - int2d(Th) (
 66                     nu * (Grad(u1,u2)' * Grad(v1,v2))
 67                 + UgradV(u1,u2, u1, u2)' * [v1,v2]
 68                 - div(u1,u2) * q
 69                 - div(v1,v2) * p
 70             )
 71             + on(1, du1=0, du2=0)
 72             ;
 73 
 74         u1[] -= du1[];
 75         u2[] -= du2[];
 76         p[] -= dp[];
 77 
 78         real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
 79         err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
 80 
 81         cout << n << " err = " << err << " " << eps << " rey = " << 1./nu << endl;
 82         if(err < eps) break; //converge
 83         if( n>3 && err > 10.) break; //blowup
 84     }
 85 
 86     if(err < eps){  //converge: decrease $\nu$ (more difficult)
 87         // Plot
 88         plot([u1, u2], p, wait=1, cmm=" rey = " + 1./nu , coef=0.3, bb=bb);
 89 
 90         // Change nu
 91         if( nu == nufinal) break;
 92         if( n < 4) cnu = cnu^1.5; //fast converge => change faster
 93         nu = max(nufinal, nu* cnu); //new viscosity
 94 
 95         // Update
 96         u1p = u1;
 97         u2p = u2;
 98         pp = p;
 99     }
100     else{   //blowup: increase $\nu$ (more simple)
101         assert(cnu< 0.95); //the method finally blowup
102 
103         // Recover nu
104         nu = nu/cnu;
105         cnu= cnu^(1./1.5); //no conv. => change lower
106         nu = nu* cnu; //new viscosity
107         cout << " restart nu = " << nu << " Rey = " << 1./nu << " (cnu = " << cnu << " ) \n";
108 
109         // Recover a correct solution
110         u1 = u1p;
111         u2 = u2p;
112         p = pp;
113     }
114 }

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.

NSNewtonTh

Fig. 25 Mesh

NSNewtonUP

Fig. 26 Velocity and pressure

Naver-Stokes newton

Table of content