Recently, I am working on simulating diffusion - convection problem for a mechanical seal.
I defined 2d-donut like domain as shown in figure below.
Configurations is as follows;
Inner region: radius = 4.1e-3[m], filled with pure water (density = 1e-6 [mol/L])
Outer region: radius = 5.5e-3[m], filled with saline (density = 0.152 [mol/L])
Velocity field (ux, uy) is acquired by solving Reynolds equation in advance.
(In this post, I used velocity u = [x, y] just for simplification)
In this case, diffusion should lead some salt from outer region to inner one,
while convection takes some salt to outer region.
I want to know how much salt intrude to the inner region against higher pressure.
I derived a weak form of diffusion - convection equation, wrote script, and tested it
(so carefully).
When I executed my program, density of some region became below 0.
Also, density (rho) at inner / outer border is not 0 / 0.152 anymore but minus.
My script:
load "iovtk"
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 rho1 = 1.000e-6; // mol/L
real rho2 = 0.154e0; // mol/L
real nu = 1.342e-9;
real dt = 1e-1;
real tmax = 1e-0;
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;}
mesh Th = buildmesh(r1(-1 * 180) + r2(360), nbvx = 2^12, fixedborder = true);
fespace Uh(Th,P2);
Uh rho, rho0, ux = x, uy = y, v;
problem convecdiff(rho, v) =
// Time differentiation term
int2d(Th)(v * (rho) / dt) - int2d(Th)(v * (rho0) / dt)
// Convection term
+ int2d(Th)(v * [ux, uy]' * grad(rho))
// Diffusion term
+ int2d(Th)(nu * grad(v)' * grad(rho))
// Dirichlet conditions
+ 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;
}
I completely do not understand
why density rho at borders can be changed while it is defined as Dirichlet condtion?
why rho can be negative while the original values (at borders) are all positive?
what is wrong in my equations and codes?
Or, should I use Neumann conditions or source / sink terms instead of Dirichlet condition?
I spent several days for this problem but I am stuck now.
I appreciate if someone could give me advice.
Your diffusion term is very weak at this value of nu. You may have some numerical instability. Try increasing nu (e.g. nu=1.0) to get a stable condition and work from there
Thank you so much, cmd.
I had performed my simulation with nu = 1.0 and confirmed that
rho at borders are equal to rho1 and rho2 respectively
all through the simulation!
I had thought that diffusion term is not problematic since
setting ux = uy = 0 did not change rho at the outer border…
Unfortunately, I can not change value of nu so largely since it
came from diffusion coefficient of natrium ion (Na+) in pure water.
Now I am going to research how I should cope with the instability.
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!
