I am solving the 1D diffusion equation by discretization using the method of lines. My problem is that I don't manage to ensure mass conservation. I have read many similar questions about the topic but it seems that the problems are related to advection terms (which I don't have here). I actually believe that my problem is related to the boundary condition implementation as changing the partition coefficient has an impact on how mass is conserved.

The equation I am trying to solve is represents the release/uptake of a component from a membrane or film into a liquid

$$ {\partial c\over\partial t}=D{\partial^2 c\over\partial x^2} \text{ for t>0 and 0<x<H} \quad[1] $$

$$ c(x,0)=1 \quad[2], \text{ initial condition} $$

$$ {dc\over dx} \mid_{(0,t)} =0\ \text{ for t>0,} \:[3], \text{ no flux on one side as it is insulated} $$

$$ D{dc\over dx} \mid_{(H,t)} =k_L(C^*_L(t)-C\mid_{(H,t)}) \quad[4],\text{ known flux on the other side} $$

$$ {dc_L\over dt} =-k_L a(C^*_L(t)-C\mid_{(H,t)}) \quad[5], \text{ mass conservation in the liquid} $$

$ \text{with } a=\frac{A}{V_L}=\frac{\text{exchange area}}{\text{liquid volume}} \quad[6] $

$ \text{ and } C^*_L = C_L ·K \quad[7], \text{ being K the partition coefficient defined as} \quad[7] $

$ K=\frac{ C^{\infty}_F }{C^{\infty}_L}=\frac{\text{concentration in the film after a long time}}{\text{concentration in the liquid after a long time}} \quad[8] $

I have discretized the previous equations as follows. Let the index $i$ the space coordinate $x_i$ so that $h=H/N$ and $i=0,1,\ldots,N$.

$$ \frac{dc_0}{dt}=D \frac{-2c_0+2c_1}{h^2} \quad [9] $$ $$ {dc_i\over dt} =D{ c_{i+1}-2c_i+c_{i-1}\over h^2} \text{ for i = 1...N-1} \quad[10] $$ and $$ \frac{dc_N}{dt}=D \frac{-2c_N+2c_{N-1}}{h^2}+\frac{2k_L}{h}(C^*_L-c_N) \quad [11] $$

I have evaluated the conservation of mass as the sum of the mass in the film and in the liquid. I have realised that several things did not affect the mass conservation:

  • Increasing or decreasing the number of nodes
  • Increasing or decreasing the diffusivity
  • Increasing or decreasing the mass transfer coefficient ($k_L$)
  • Increasing or decreasing the specific exchange area ($a$)

But the following did affect how mass was conserved:

  • Changing the partition coefficient ($K$)
  • Changing the film thickness ($H$)

What do you mean by "I have realised that several things did not affect the mass conservation". How did you measure conservation ? Are you asking if your method conserves total mass ?

My answer below is based on a conservation law I notice in your model.

You wrote an $r$ in your PDE but did not use it later, so I am assuming it is zero. Your model has this quantity conserved $$ m(t) = a \int_0^H c(x,t)dx + c_L(t) $$ wrt time, i.e., $$ m(t) = \textrm{constant} $$ To approximate the integral, we will use trapezoidal rule $$ m = a\left(c_0/2 + \sum_{i=1}^{N-1} c_i + c_N/2\right) h + c_L $$ Then this changes as $$ dm/dt = \frac{aD}{h} \left(c_1 - c_0 + \sum_{i=1}^{N-1} (c_{i-1} - 2 c_i + c_{i+1}) + c_{N-1} - c_N\right) + a k_L(C_L^* - c_N) - a k_L (C_L^* - c_N) $$ By summing up the terms inside the big bracket, you see it is zero, so that $$ dm/dt = 0 $$ Hence the total "mass" $m$ is conserved.

I have answered a similar question here Neumann boundary conditions diffusion equations methods of lines

  • $\begingroup$ It makes sense but I am not entirely sure of the role of the specific exchange area a . In my model it appears only in eq. [5] but not in equations [9]-[11]. Is it considered implicitely somehow? I am a bit confused. $\endgroup$
    – Toulousain
    Feb 26 '18 at 9:29
  • $\begingroup$ I found the mistake in my code (I miscalculated a). That's such an annoying little mistake for a complex problem, but I guess that's how it goes. Still, I believe your answer was very helpful in letting me discard discretization problems. $\endgroup$
    – Toulousain
    Feb 26 '18 at 10:36

The method of lines is not applicable to this equation (which is called advection-diffusion equation) in the same sense as it is applicable to the advection equation. There are no paths x(t) along which c is constant.

I would use a finite-difference approach throughout. You have already written down a centered finite difference in [10]; just use a forward-time difference, and then write down finite difference expressions for the spatial derivatives at the two boundaries. I could fill in details, but you seem to be well-familiar with how this works. For a forward time-difference, there will be a stability condition on the time step. As long as D and r are constant, this should automatically conserve mass.

One could also use a spectral solver, but implementing the boundary conditions in Fourier space would be some work.

  • $\begingroup$ I assume that you propose a Crank-Nicholson scheme or similar. I could certainly do that but I still do not understand why I cannot use the method of lines. Plus, I might want to implement cases (later) where D is not constant, but, of course, I need to solve this simpler case first. $\endgroup$
    – Toulousain
    Feb 26 '18 at 9:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.