# Mass conservation in 1d diffusion by method of lines

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

• 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. – Toulousain Feb 26 '18 at 9:29
• 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. – 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.

• 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. – Toulousain Feb 26 '18 at 9:33