Solving non linear system of equations with Newton Raphson method

General Discussion
1

Hi , I am trying to solve steady boussinesq approximation of governing equations of flow compromising conservation of mass, momentum, energy. I weak formulated the equation and I need to use Newton raphson method to solve the non-linear residuals. I am requesting some guidance in how to solve problem of multiple coupled non linear residuals.

2

You have an example with comments in the FreeFem documentation p.75 :
Newton Method for the Steady Navier-Stokes equations

3

Thanks for your reply. I tried it. But I am implemented the following code and I made an error here. How to correct and identify it?

// Parameters
real Pr = 0.7;
real Ra = 1000;
real gamma = 100; // penalty parameter
real eps = 1e-5;

// Mesh
mesh Th = square(33, 33); // Labels: left=1, right=2, bottom=3, top=4
plot(Th, wait=true);

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

fespace Mh(Th, P1);
Mh T, q, dT, pp;

// Stream Function Space
fespace Ph(Th, P2); 
Ph psi, phi;

// Macros
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)) //
macro DivInner(u1, u2, v1, v2) ((dx(u1) + dy(u2)) * (dx(v1) + dy(v2))) //
macro Graddotgrad(T, q) (dx(q) * dx(T) + dy(q) * dy(T)) //
macro Advect(u1, u2, T) (u1 * dx(T) + u2 * dy(T)) //

// Initialization
u1 = 0;
u2 = 0;
T = 0;

// Viscosity loop
while (1) {
    int n;
    real err = 0;

    // Newton loop
    for (n = 0; n < 15; n++) {
        // Solve velocity equation
        solve velocity ([du1, du2], [v1, v2])
            = int2d(Th)(
                Pr * (Grad(du1, du2)' * Grad(v1, v2))
                + UgradV(du1, du2, u1, u2)' * [v1, v2]
                + UgradV(u1, u2, du1, du2)' * [v1, v2]
                + gamma * DivInner(du1, du2, v1, v2) // divergence penalty
                - Ra * Pr * dT * v2
		)
		-int2d(Th)( Pr * (Grad(u1, u2)' * Grad(v1, v2))
                + UgradV(u1, u2, u1, u2)' * [v1, v2]
                + gamma * DivInner(u1, u2, v1, v2) // divergence penalty
                - Ra * Pr * T * v2
            )
            + on(1, u1=0, u2=0)
            + on(2, u1=0, u2=0)
            + on(3, u1=0, u2=0)
            + on(4, u1=0, u2=0);

        // Solve temperature equation
        solve Temperature ([dT], q)
            = int2d(Th)(
                Advect(du1, du2, T) * q
                + Advect(u1, u2, dT) * q
                + Graddotgrad(dT, q)
            )
            - int2d(Th)(
                Advect(u1, u2, T) * q
                + Graddotgrad(T, q)
            )
            + on(3, T = 1) // Hot bottom wall
            + on(4, T = 0); // Cold top wall

        // Update fields
        u1[] -= du1[];
        u2[] -= du2[];
        T[] -= dT[];

        // Compute error with safeguard
        real Lu1 = max(1e-10, u1[].linfty);
        real Lu2 = max(1e-10, u2[].linfty);
        real LT = max(1e-10, T[].linfty);
        err = du1[].linfty / Lu1 + du2[].linfty / Lu2 + dT[].linfty / LT;

        // Stop if convergence is reached
        if (err < eps) break;
    }

    // Check convergence
    if (err < eps) {
        cout << "Converged after " << n << " Newton iterations.\n";
        break;
    }
}
4

There are several mistakes in your code:

– default labels bottom=1, right=2, top=3, left=4

– in a variational formulation you have to separate (ie in different int2d()) the terms linear with respect to the unknown and the terms which do not depend on the unknown

– boundary conditions on() : only the unknowns are allowed (du1,du2) in “velocity”, dT in “Temperature”

– it is better to resolve both velocity and temperature together in order to have fast convergence of the Newton iteration

– You have to set the boundary conditions on the updated values u1-du1, u2-du2, T-dT. Thus to set these values to 0,0,1 respectively you need to set

+ on(1,2,3,4, du1=u1, du2=u2)  //(meaning u1-du1=0, u2-du2=0)
+ on(1, dT = T-1) //(meaning T-dT=1)

boussinesq.edp (2.6 KB)

I am doubtful of the penalty method for div u=0. Probably it gives wrong results because of no inf-sup stability. You have to use rather a variational formulation in velocity-pressure (and temperature in your case).

5

Very thanks sir! Yes, the results are wrong. I am not aware of the stability condition. Without penalty method, I am not getting convergence. What should I look into?

6

The inf-sup condition of Babuska-Brezzi is very classical, see for example

see also “the system of Stokes for fluids” in FF documentation p.69

You have to use a couple of spaces for velocity and pressure satisfying the inf-sup condition. You can take for example P2 for velocity and P1 for pressure (Taylor-Hood choice). This gives
boussinesq.edp (3.0 KB)