Lid-driven cavity problem for MPFF code gives error for a modified NSEs model

General Discussion
1

Hi,

Dear respected Professors @fb77 @julienG , thank you for your kind guidance and help.

I am struggling to fix this bug in my code, please help to make the code run smoothly.

NSEs modified momentum model:

image
image775×50 2.37 KB

error:

image
image962×490 9.39 KB

code snippet:

Blockquote
// Micropolar fluid flow with viscoelastic tensor: Lid-driven cavity flow test using DEVSS stabilization

real n = 100;
mesh Th = square(n, n);
plot(Th);
fespace Vh(Th, P2);
fespace Mh(Th, P2);
real al = 0.5, dt = 0.002, dT = dt / 2.0, t = 0, T = 15.0;
real re = 1000.0, ee = 0.0001;
real rhonf = 1.73;
real munf = 1.0;
real beta1 = 0.5;
real beta2 = 0.5;
real sigmanf = 1.05;
real B0 = 3.0;
real Omega = 1.025;
real lambda = 1.05;

real number = 0;

Mh P, p, q;
Vh w, u, v, uT, vT, uTT, vTT, omega, omegaT, omegaTT;

real nv = Th.nv;
real nt = Th.nt;

real[int] areaphi(nt), areahphi(nt), h(nt);

u = 0; // Initial condition for velocity u
v = 0; // Initial condition for velocity v
omega = 1.0; // Initial condition for microrotation omega

for (int m = 0; m <= T / dt; m++) {
Vh uold = u, vold = v, pold = p, omegaold = omega;
P = p;

// Intermediate velocity prediction (uT, vT)
solve Adual4unon(uT, w) =
    int2d(Th)(rhonf * uT * w / dt)
    - int2d(Th)(rhonf * convect([uold, vold], -dt, uold) * w / dt)
    + int2d(Th)(2 * rhonf * Omega * (dy(vold) * w - dx(uold) * w))
    + int2d(Th)(munf * ((1 + beta1) / (1 + beta2)) * (dx(uT) * dx(w) + dy(uT) * dy(w)))
    - int2d(Th)(sigmanf * B0^2 * uT * w)
    + int2d(Th)(lambda * (dx(uold) * dx(w) + dy(uold) * dy(w)))  // Extra stress tensor term
    + int2d(Th)((dx(uT) - dx(uold)) * dx(w) + (dy(uT) - dy(uold)) * dy(w))  // DEVSS stabilization term
    + on(1, 2, 4, uT = 0)
    + on(3, uT = 1.0);

solve Adual4vnon(vT, w) =
    int2d(Th)(rhonf * vT * w / dt)
    - int2d(Th)(rhonf * convect([uold, vold], -dt, vold) * w / dt)
    + int2d(Th)(2 * rhonf * Omega * (dy(uold) * w - dx(vold) * w))
    + int2d(Th)(munf * ((1 + beta1) / (1 + beta2)) * (dx(vT) * dx(w) + dy(vT) * dy(w)))
    - int2d(Th)(sigmanf * B0^2 * vT * w)
    + int2d(Th)(lambda * (dx(vold) * dx(w) + dy(vold) * dy(w)))  // Extra stress tensor term
    + int2d(Th)((dx(vT) - dx(vold)) * dx(w) + (dy(vT) - dy(vold)) * dy(w))  // DEVSS stabilization term
    + on(1, 2, 3, 4, vT = 0);

// Solve for microrotation omegaT
solve Adual4omeganon(omegaT, w) =
    int2d(Th)(omegaT * w / dt)
    - int2d(Th)(convect([uold, vold], -dt, omegaold) * w / dt)
    + on(1, 2, 3, 4, omegaT = 0);

// Pressure correction step
solve pb4p(p, q) =
    int2d(Th)((dx(p) * dx(q) + dy(p) * dy(q)))
    + int2d(Th)(((dx(uT) + dy(vT)) * q) / dt);

// Final velocity update (u, v)
solve pb4u(u, w) =
    int2d(Th)((dx(u) * dx(w) + dy(u) * dy(w)) / re + u * w / dt)
    - int2d(Th)(uT * w / dt - dx(p) * w)
    + int2d(Th)((dx(u) - dx(uold)) * dx(w) + (dy(u) - dy(uold)) * dy(w))  // DEVSS stabilization term
    + on(1, 2, 4, u = 0)
    + on(3, u = 1.0);

solve pb4v(v, w) =
    int2d(Th)((dx(v) * dx(w) + dy(v) * dy(w)) / re + v * w / dt)
    - int2d(Th)(vT * w / dt - dy(p) * w)
    + int2d(Th)((dx(v) - dx(vold)) * dx(w) + (dy(v) - dy(vold)) * dy(w))  // DEVSS stabilization term
    + on(1, 2, 3, 4, v = 0);

// Solve for final microrotation omega
solve pb4omega(omega, w) =
    int2d(Th)((dx(omega) * dx(w) + dy(omega) * dy(w)) / re + omega * w / dt)
    - int2d(Th)(omegaT * w / dt)
    + on(1, 2, 3, 4, omega = 0);

// Compute stress tensor components
Vh tauxx = lambda * 2 * dx(u);
Vh tauxy = lambda * (dy(u) + dx(v));
Vh tauyy = lambda * 2 * dy(v);

t = t + dt;

if (m % 100 == 0) {
    ofstream os("lemma1_tau" + real(t) + ".plt");
    real m1, m2, m3, m4, m5, m6, m7, m8, m9;
    os << "TITLE = levelset " << endl;
    os << "VARIABLES = X, Y, u, v, omega, p, tauxx, tauxy, tauyy" << endl;
    os << "ZONE, N=" << nv << ", E=" << nt << ", F=FEPOINT, ET=triangle" << endl;
    for (int i = 0; i < nv; i++) {
        m1 = Th(i).x;
        m2 = Th(i).y;
        m3 = u(m1, m2);
        m4 = v(m1, m2);
        m5 = omega(m1, m2);
        m6 = p(m1, m2);
        m7 = tauxx(m1, m2);
        m8 = tauxy(m1, m2);
        m9 = tauyy(m1, m2);
        os << m1 << " " << m2 << " " << m3 << "  " << m4 << "  " << m5 << "  " << m6 << " " << m7 << " " << m8 << " " << m9 << endl;
    }
    os << "\n\n\n\n" << endl;
    int a2, bt2, c2;
    for (int i = 0; i < nt; i++) {
        a2 = Th[i][0]
        bt2 = Th[i][1] + 1;
        c2 = Th[i][2] + 1;
        os << a2 << " " << bt2 << " " << c2 << endl;
    }
}

}

 //**************************output results************************************

Blockquote

2

In your variational formulations you must separate the terms containing the unknowns and the others. This means fror Adual4unon(uT, w)

+int2d(Th)(...)// terms containing uT
+int2d(Th)(...)// terms not containing uT

In particular the term dx(uT)*dx(w) must be separated from the term -dx(uold)*dx(w)

1 Like