Thank you for your insightful comment, frederichecht.
I had known the word “Peclet number” but never knew that larger Peclet number causes instability and alters boundary values.
With a help of your comments, I have tried to avoid numerical instability by
- introducing “convect” function,
- adding additional inner & outer regions,
- using finer mesh,
and - increasing dt value.
Finally, I could get somewhat reasonable results as below. Now I can go forward and will investigate further. I appreciate your comments so much!
load "iovtk"
load "ff-Ipopt"
load "UMFPACK64"
macro grad(A) [dx(A), dy(A)] //
macro div(A) (dx(A) + dy(A)) //
real rr1 = 4.1e-3; // inner diameter = 8.2mm.
real rr2 = 5.5e-3; // outer diameter = 11mm.
real rr1A = rr1 - 0.1e-3;
real rr2A = rr2 + 0.1e-3;
real rho1 = 1.000e-6; // mol/L
real rho2 = 0.154e0; // mol/L
real nu = 1.342e-9;
real dt = 1e+1;
real tmax = 2e+2;
int n = 4;
int tmesh1 = ltime();
border r1(theta = 0, 2*pi) {x = rr1 * cos(theta); y = rr1 * sin(theta); label = 1;}
border r2(theta = 0, 2*pi) {x = rr2 * cos(theta); y = rr2 * sin(theta); label = 2;}
border r1A(theta = 0, 2*pi) {x = rr1A * cos(theta); y = rr1A * sin(theta); label = 1;}
border r2A(theta = 0, 2*pi) {x = rr2A * cos(theta); y = rr2A * sin(theta); label = 2;}
mesh Th0 = buildmesh(r1(-1 * 180 * n) + r2(360 * n), nbvx = 2^16, fixedborder = true);
mesh Th1 = buildmesh(r1A(-1 * 180 * n) + r1(180 * n), nbvx = 2^16, fixedborder = true);
mesh Th2 = buildmesh(r2(-1 * 360 * n) + r2A(360 * n), nbvx = 2^16, fixedborder = true);
mesh Th = Th0 + Th1 + Th2;
fespace Uh(Th,P2);
Uh rho, rho0, ux = x, uy = y, v;
problem convecdiff(rho, v) =
int2d(Th)(((1.0/dt) * rho * v) + (nu * grad(rho)' * grad(v)))
- int2d(Th)((1.0/dt) * convect([ux, uy], -dt, rho0) * v)
+ on(r1, rho=rho1) + on(r2, rho=rho2) ;
int k = 1;
for (real t = 0.0; t <= tmax; t = t + dt) {
rho0 = rho;
convecdiff;
savevtk("results/data" + k + ".vtk", Th,
rho, [ux, uy, 0],
dataname = "rho u",
bin = true);
k = k + 1;
}