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)$$:

 $$(\mathbf{u}\cdot\nabla) \mathbf{u}-\nu \Delta \mathbf{u}+\nabla p$$ $$=$$ $$0,\\ \nabla\cdot \mathbf{u}$$ $$=$$ $$0$$



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 some boundary conditions ( $$\mathbf{u}$$ is given on $$\Gamma$$)

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

 $$\int_\Omega ((\mathbf{u}\cdot\nabla) \mathbf{u} ). \bm 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$$)
• solve $$DF(u_i) w_i = F(u_i)$$;
• $$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 : [-

 $$F(\mathbf{u},p) = \int_\Omega$$ $$((\mathbf{u}\cdot\nabla) \mathbf{u} ). \bm 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} ). \bm v + ((\mathbf{u}\cdot\nabla) \mathbf{\delta u} ). \bm v \\$$ $$+$$ $$\nu \nabla \mathbf{\delta u}:\nabla \mathbf{v} - \delta p \nabla\cdot \mathbf{v} - q \nabla\cdot \mathbf{\delta u}$$



So the Newton algorithm become Example[NSNewton.edp]

 ...
for( n=0;n< 15;n++)
{ solve Oseen([du1,du2,dp],[v1,v2,q]) =
- div(du1,du2)*q - div(v1,v2)*dp
- 1e-8*dp*q // stabilization term
)
- div(u1,u2)*q - div(v1,v2)*p
)
+ on(1,du1=0,du2=0) ;
u1[] -= du1[];  u2[] -= du2[]; p[]  -= dp[];
err= du1[].linfty + du2[].linfty + dp[].linfty;
if(err < eps) break;
if( n>3 \)||

     }


With the operator:

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


We build a computation mesh the exterior of a 2d cylinder.

 real R = 5,L=15;
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?;

/  FE Space Taylor Hood
fespace Xh(Th,P2);// for volicity
fespace Mh(Th,P1);// for pressure
Xh u1,u2,v1,v2,du1,du2,u1p,u2p;
Mh p,q,dp,pp;

// intial guess with B.C.
u1 = ( x^2+y^2) > 2;
u2=0;


Finally we use trick to make continuation on the viscosity $$\nu$$, because the Newton method blowup owe start with the final viscosity $$\nu$$

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

// stop test for Newton
real eps=1e-6;

verbosity=0;
while(1)  //  Loop on viscosity
{   int n;
real err=0; // err on Newton algo ...

... put the new the Newton  algo here

if(err < eps)
{ // converge decrease $\nu$ (more difficult)
plot([u1,u2],p,wait=1,cmm=" rey = " + 1./nu , coef=0.3,bb=bb);
if( nu == nufinal) break;
if( n < 4) cnu=cnu^1.5; // fast converge => change faster
nu = max(nufinal, nu* cnu); // new vicosity
u1p=u1;  u2p=u2;  pp=p; //  save correct solution ...
}
else
{   //  blowup increase $\nu$  (more simple)
assert(cnu< 0.95); //  the method  finally  blowup
nu = nu/cnu; //  get previous value of viscosity
cnu= cnu^(1./1.5); // no conv. => change lower
nu = nu* cnu;  // new viscosity
cout << " restart nu = " << nu << " Rey= "<< 1./nu << "  (cnu = " << cnu << " ) \nabla";
// restore a correct solution ..
u1=u1p;
u2=u2p;
p=pp;
}
}


[htbp]

[NSNewton] Mesh and the velocity and pressure at Reynolds $$200$$