Poiseuille Flow Test with Extra Stress Tensor: Fatal error in Convect (R2) operator: loop => velocity too high ? or NaN F. Hecht

Dear respected Prof @fb77 thank you so much for the enormous help you give on this platform for all of us. Here again I face some issue with my code for Poiseuille Flow Test as following:
I am solving micropolarpolar fluid flow standard NSEs model with viscoelastic stress tensor in this code but I always run into this issue:

{ – Solve :
min -19.3012 max 13.117
– Solve :
min 0 max 0
– Solve :
min -2.99727e+13 max -2.99727e+13
Fatal error in Convect (R2) operator: loop => velocity too high ??? or NaN F. Hecht
current line = 110
Assertion fail : (0)
line :4939, in file lgfem.cpp
Assertion fail : (0)
line :4939, in file lgfem.cpp
err code 6 , mpirank 0}

I fixed the initial and BCs as in the standard Poiseuille Flow Test, but still the issue continues.
code:

{real uMax = 1.0;
real omegaMax = 1.0;
real Mu = 0.01;
real Kappa = 0.005;
real Gamma = 0.001;
real J = 0.01;
real dpdx = -1.01;
real beta = 1.0 / 9.0, Wi = 1.0; // Weissenberg number

real al = 0.5, dt = 0.001, dT = dt / 2.0, t = 0, T = 20.0, re = 100.0, ee = 0.0001;

int nn = 16;
real L = 8.0;
real D = 1.0;
int Inlet = 1;
int Wall = 2;
int Outlet = 3;

// Define velocity profile
func uIn = 0.4 * y * (1 - y); // Updated velocity profile
func omegaIn = omegaMax * 0.1 * D * (1 - D);

border b1(t = 0., 1.) { x = L * t; y = 0.; label = Wall; };
border b2(t = 0., 1.) { x = L; y = D * t; label = Outlet; };
border b3(t = 0., 1.) { x = L - L * t; y = D; label = Wall; };
border b4(t = 0., 1.) { x = 0.; y = D - D * t; label = Inlet; };

// FEM Mesh
int nnL = max(2., L * nn * 2);
int nnD = max(2., D * nn * 2);

mesh Th = buildmesh(b1(nnL) + b2(nnD) + b3(nnL) + b4(nnD));
plot(Th, cmm = “mesh”, ps = “PoiseuilleFlowmesh.eps”, fill = 0, value = 1, wait = 1);

fespace Vh(Th, P2);
fespace Mh(Th, P2);

Mh P, p, q;
Vh w, u, v, omega, tau11, tau12, tau21, tau22, G11, G12, G21, G22;

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

cout << "Nodes: " << nv << endl;
cout << "Triangles: " << nt << endl;

for (int m = 0; m <= T / dt; m++) {

Vh uOld = u, vOld = v, pOld = p, omegaOld = omega;
Mh P = p;

// Solve for intermediate velocity
solve Adual4unon(u, w) =
    int2d(Th)(u * w / dt)
    - int2d(Th)(convect([uOld, vOld], -dt, uOld) * w / dt)
    + int2d(Th)((dx(tau11) + dy(tau12)) * w)  // Viscoelastic contribution
    + on(Wall, u = 0)
    + on(Inlet, u = uIn)
    + on(Outlet, u = 0);

solve Adual4vnon(v, w) =
    int2d(Th)(v * w / dt)
    - int2d(Th)(convect([uOld, vOld], -dt, vOld) * w / dt)
    + int2d(Th)((dx(tau21) + dy(tau22)) * w)  // Viscoelastic contribution
    + on(Wall, v = 0)
    + on(Inlet, v = 0)
    + on(Outlet, v = 0);

solve Adual4omeganon(omega, w) =
    int2d(Th)(omega * w / dt)
    - int2d(Th)(convect([uOld, vOld], -dt, omegaOld) * w / dt)
    + on(Wall, omega = 0)
    + on(Inlet, omega = omegaIn)
    + on(Outlet, omega = 0);

solve pb4p(p, q) =
    int2d(Th)((dx(p) * dx(q) + dy(p) * dy(q)))
    + int2d(Th)(((dx(u) + dy(v)) * q) / dt);

Vh tau11old = tau11, tau12old = tau12, tau21old = tau21, tau22old = tau22;
Vh G11old = G11, G12old = G12, G21old = G21, G22old = G22;

// ********************************** Solve the constitutive model equation ********************* //


Vh dudx = dx(u); // Compute derivative inside the domain
solve Oldroyd11 (tau11, w) =
int2d(Th)(tau11 * w / dt)
+ int2d(Th)(-1 / dt * convect([u, v], -dt, tau11old) * w)
- int2d(Th)(((1 - beta) * 2 / Wi * dudx + 2 * (dudx * tau11old + dy(u) * tau12old) - tau11old / Wi) * w)
+ on(Wall, tau11 = 0)
+ on(Inlet, tau11 = 2.0 * Wi * (1 - beta) * dudx) // Now using precomputed dudx
+ on(Outlet, tau11 = 0);


solve Oldroyd12 (tau12, w) =
    int2d(Th)(tau12 * w / dt)
    + int2d(Th)(-1 / dt * convect([u, v], -dt, tau12old) * w)
    - int2d(Th)(((1 - beta) / Wi * (dy(u) + dx(v)) + (dx(u) * tau12old + dy(u) * tau22old + dy(v) * tau12old + dx(v) * tau11old)) * w - w * tau12old / Wi)
    + on(Wall, tau12 = 0)
    + on(Inlet, tau12 = (1.0 - beta) * dx(v))
    + on(Outlet, tau12 = 0);

tau21 = tau12;

solve Oldroyd22 (tau22, w) =
    int2d(Th)(tau22 * w / dt)
    + int2d(Th)(-1 / dt * convect([u, v], -dt, tau22old) * w)
    - int2d(Th)(((1 - beta) * 2 / Wi * dy(v) + 2 * (dx(v) * tau12old + dy(v) * tau22old) - tau22old / Wi) * w)
    + on(Wall, tau22 = 0)
    + on(Inlet, tau22 = 2.0 * Wi * (1 - beta) * dy(v))
    + on(Outlet, tau22 = 0);

t = t + dt;

if (m % 100 == 0) {
ofstream os("Bmmpff_bearlub_new" + 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 = tau11(m1, m2); // This corresponds to tauxx
    m8 = tau12(m1, m2); // This corresponds to tauxy
    m9 = tau22(m1, m2); // This corresponds to tauyy
    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] + 1;
    bt2 = Th[i][1] + 1;
    c2 = Th[i][2] + 1;
    os << a2 << " " << bt2 << " " << c2 << endl;
	}

}
}
}

Dear Hussein,
Thank you for your return.
There are several mistakes in your code: sign error in div(tau), missing projection step for the pressure. And the variables are not updated at the right time.
The scheme is unstable with separate resolution of velocity and stress. It is more stable to solve all variables simultaneously.
I added a boundary term on the outlet in order to be compatible with the Poiseuille flow.
micropl.edp (8.6 KB)

1 Like

Thank you prof for this valuable help