Convection-diffusion problem with pre-computed velocity field

General Discussion
1

Dear FreeFem users,

I’m trying to solve a convection-diffusion problem of this form:

\partial_t c - K\triangle c + \nabla\cdot(u c) - B(cv - c)=0

where K is the diffusion coefficient and it is constant, the term B(cv - c) is a distributed source in a region of the domain and u is the velocity field which has been previously computed (from Darcy’s law). In figure you can see u.

Screenshot 2020-11-01 at 21.28.55
Screenshot 2020-11-01 at 21.28.551382×992 19.8 KB

I have tried to solve the equation in this form:
\partial_t c - K\triangle c + u\cdot\nabla c + (\nabla\cdot u )c - B(cv - c)=0

by using the Freefem convect operator on the time derivative and the term u\cdot\nabla c.
Here it is the piece of code:

macro grad(fun) [dx(fun), dy(fun)] //eom
Vh divvel = dx(uxh) + dy(uyh);
problem Cevol(C, v, solver=GMRES) =
int2d(Th)(1./dtCv + Kgrad©’grad(v))
- int2d(Th)(1./dtconvect([uxh, uyh], -dt, Cold)v)
- int2d(Th)(Colddivvelv)
- int2d(Th)(BCvv)
+ int2d(Th)(BCv)
;
I would like to know if this is the correct way to treat the convection term of such an equation. The thing is that in my results the problem is solved with an acceptable solution when I discard the term - int2d(Th)(Colddivvelv), which is not correct in my opinion, since I do not have the condition \nabla\cdot u =0.
If I consider the - int2d(Th)(Colddivvelv) in the problem, the solution explodes in few iterations in time.

Do you have any suggestion?

Thank you very much.

Best regards,
Gabriele

2

In this schema explicite in the part c(\nabla.u) if (\nabla.u) is large the the schema
become instable son the first thing is to implicite this term and warning to with the signe of
int2d(Th)(Colddivvel v)
it is a + not a -.

problem Cevol(C, v, solver=sparsesolver) =
int2d(Th)(1./dt*C* v + K*grad(c)'*grad(v) + C*divvel* v)
- int2d(Th)(1./dt*convect([uxh, uyh], -dt, Cold)*v)
- int2d(Th)(B*Cv* v)
+ int2d(Th)(B*C* v)
;

Best Regards, F. Hecht