Hi respected Professor, and thank you so much for your valuable and very kind guidances. I get some errors when I modified my original Lid-driven cavity flow code for mpff specific and modified version of NSEs as following:
FF++ Error Prompt:
Modified NSEs:
Bugged Code:
[// Micropolar fluid flow in 2D of Lemma 1.
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 = 1.0, re = 1000.0, ee = 0.0001;
real number = 0;
Mh P, p, q;
Vh w, u, v, uT, vT, uTT, vTT, omega, omegaT, omegaTT;
real rho_nf = 1.0; // Density
real beta1 = 0.5, beta2 = 0.5;
real mu_nf = 1.0; // Viscosity
real sigma_nf = 1.0; // Conductivity
real B0 = 1.0; // Magnetic field strength
real[3] Omega = [0.0, 0.0, 1.0]; // Rotation vector
real nv = Th.nv;
real nt = Th.nt;
real[int] areaphi(nt), areahphi(nt), h(nt);
for (int m = 0; m <= T / dt; m++) {
Vh uold = u, vold = v, pold = p, omegaold = omega;
P = p;
// Intermediate velocity uT
solve Adual4unon(uT, w) =
int2d(Th)(uT * w / dt)
- int2d(Th)(convect([uold, vold], -dt, uold) * w / dt)
+ int2d(Th)(2 * Omega[2] * v * w) // Coriolis term
- int2d(Th)(sigma_nf * B0^2 * u * w) // Magnetic field term
+ int2d(Th)((mu_nf * (1 + beta1) / (1 + beta2)) * dx(dx(u)) * w) // Modified viscosity term
+ on(1, 2, 4, uT = 0)
+ on(3, uT = 1.0);
// Intermediate velocity vT
solve Adual4vnon(vT, w) =
int2d(Th)(vT * w / dt)
- int2d(Th)(convect([uold, vold], -dt, vold) * w / dt)
- int2d(Th)(2 * Omega[2] * u * w) // Coriolis term
- int2d(Th)(sigma_nf * B0^2 * v * w) // Magnetic field term
+ int2d(Th)((mu_nf * (1 + beta1) / (1 + beta2)) * dx(dx(v)) * w) // Modified viscosity 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);
// Solve for pressure p
solve pb4p(p, q) =
int2d(Th)((dx(p) * dx(q) + dy(p) * dy(q)))
+ int2d(Th)(((dx(uT) + dy(vT)) * q) / dt);
// Solve for uTT
solve pb4uTT(uTT, w) =
int2d(Th)(uTT * w / dt)
- int2d(Th)(uT * w / dt - dx(p) * w)
+ on(1, 2, 4, uTT = 0)
+ on(3, uTT = 1.0);
// Solve for vTT
solve pb4vTT(vTT, w) =
int2d(Th)(vTT * w / dt)
- int2d(Th)(vT * w / dt - dy(p) * w)
+ on(1, 2, 3, 4, vTT = 0);
// Solve for microrotation omegaTT
solve pb4omegaTT(omegaTT, w) =
int2d(Th)(omegaTT * w / dt)
- int2d(Th)(omegaT * w / dt - dx(p) * w)
+ on(1, 2, 3, 4, omegaTT = 0);
// Update velocity u
solve pb4u(u, w) =
int2d(Th)((dx(u) * dx(w) + dy(u) * dy(w)) / re + u * w / dt)
- int2d(Th)(uTT * w / dt)
+ on(1, 2, 4, u = 0)
+ on(3, u = 1.0);
// Update velocity v
solve pb4v(v, w) =
int2d(Th)((dx(v) * dx(w) + dy(v) * dy(w)) / re + v * w / dt)
- int2d(Th)(vTT * w / dt)
+ 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)(omegaTT * w / dt)
+ on(1, 2, 3, 4, omega = 0);
t = t + dt;
if (m % 100 == 0) {
ofstream os("MPFF" + real(t) + ".plt");
real m1, m2, m3, m4, m5, m6;
os << "TITLE = levelset " << endl;
os << "VARIABLES = X, Y, u, v, omega, p" << 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);
os << m1 << " " << m2 << " " << m3 << " " << m4 << " " << m5 << " " << m6 << 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;
}
}
}
]
