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