# I am on trouble with convection (advection) - diffusion problem

Dear FreeFEM users,

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 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
// 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

1. why density rho at borders can be changed while it is defined as Dirichlet condtion?
2. why rho can be negative while the original values (at borders) are all positive?
3. what is wrong in my equations and codes?
4. 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

1 Like

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.

This scheme is instable, this is why you get bad value,
The Peclet (Pe) number of the test is to huge , and I do not understand

Pe = L V / nu , L= rr2-rr1 = 1.4e-3 , V = 5e-3 ,nu = 1.3e-9

Pe = 70 e-6/ 1.3 e-9 ~ 53 000

=> a boundary layer very thin.

1 Like

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

1. introducing “convect” function,
3. using finer mesh,
and
4. 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"

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) * 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;
}
``````
1 Like