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
2real R = 5.;
3real L = 15.;
4
5real nu = 1./50.;
6real nufinal = 1/200.;
7real cnu = 0.5;
8
9real eps = 1e-6;
10
11verbosity = 0;
12
13// Mesh
14border cc(t=0, 2*pi){x=cos(t)/2.; y=sin(t)/2.; label=1;}
15border ce(t=pi/2, 3*pi/2){x=cos(t)*R; y=sin(t)*R; label=1;}
16border beb(tt=0, 1){real t=tt^1.2; x=t*L; y=-R; label=1;}
17border beu(tt=1, 0){real t=tt^1.2; x=t*L; y=R; label=1;}
18border beo(t=-R, R){x=L; y=t; label=0;}
19border bei(t=-R/4, R/4){x=L/2; y=t; label=0;}
20mesh Th = buildmesh(cc(-50) + ce(30) + beb(20) + beu(20) + beo(10) + bei(10));
21plot(Th);
22
23//bounding box for the plot
24func bb = [[-1,-2],[4,2]];
25
26// Fespace
27fespace Xh(Th, P2);
28Xh u1, u2;
29Xh v1,v2;
30Xh du1,du2;
31Xh u1p,u2p;
32
33fespace Mh(Th,P1);
34Mh p;
35Mh q;
36Mh dp;
37Mh pp;
38
39// Macro
40macro Grad(u1,u2) [dx(u1), dy(u1), dx(u2),dy(u2)] //
42                        [u1,u2]'*[dx(v2),dy(v2)]] //
43macro div(u1,u2) (dx(u1) + dy(u2)) //
44
45// Initialization
46u1 = (x^2+y^2) > 2;
47u2 = 0;
48
49// Viscosity loop
50while(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)(
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) (
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.

Naver-Stokes newton

Table of content