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;
}
}
}
}