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.

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./dt*C*v + K*grad©’ grad(v))*v)

- int2d(Th)(1./dtconvect([uxh, uyh], -dt, Cold)v)

- int2d(Th)(Colddivvel

- int2d(Th)(B

*Cv*v)

+ int2d(Th)(B

*C*v)

;

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)(Cold

*divvel*v), which is not correct in my opinion, since I do not have the condition \nabla\cdot u =0.

If I consider the - int2d(Th)(Cold

*divvel*v) in the problem, the solution explodes in few iterations in time.

Do you have any suggestion?

Thank you very much.

Best regards,

Gabriele