Mass conservation in radially symmetric parabolic PDE problems

[ndalchau]

Dear all,
I'm trying to solve the 2d heat equation in a radially symmetric domain, numerically using the Crank-Nicolson method. i.e.

$$\dfrac{\partial u}{\partial t} = D\left( \dfrac{\partial^2u}{\partial r^2}+\dfrac{1}{r}\dfrac{\partial u}{\partial r}\right)$$

Applying the Crank-Nicolson method basically results in a recurrence relation:

$$(-q+\frac{z}{r})u_{i-1,j+1} + (1+2q)u_{i,j+1} - (q+\frac{z}{r})u_{i+1,j+1} = (q-\frac{z}{r})u_{i-1,j} + (1-2q)u_{i,j} + (q+\frac{z}{r})u_{i+1,j}$$

where $$q=\dfrac{\delta t}{2(\delta r)^2}$$, $$z=\dfrac{\delta t}{4\delta r}$$ and $$u_{i,j}$$ is the solution at $$r=i\delta r, t=j\delta t$$.

You can write this into a matrix equation of the form $$Au\{j+1\}=Bu\{j\}$$, which basically enables you to solve the problem. This all works very nicely for the 1d heat equation (the same differential equation but without the first spatial derivative). However, I now have the problem that the solution doesn't conserve mass. That is, given some initial condition across $$r\in (0,R)$$, the sum of the solution points decreases over time.

Clearly, if the column sums of A and B are the row vector of 1s, then I am guaranteed mass conservation. However, this is not the case for my Crank-Nicolson implementation here. This results from the $$\frac{z}{r}$$ changing values and signs between neighbouring rows, which means they don't cancel.

Anyone got any comments? Is there a way of solving these parabolic PDEs in radially symmetric domains that preserves the mass conservation law? Or is my implementation incorrect maybe?

Your assistance would be greatly appreciated,

Neil

 

[jvicens]

Dear Neil,

Thank you for sharing your question with the forum. It sounds like you are working on a very interesting and challenging problem.

One possible solution to your issue with mass conservation could be to modify the Crank-Nicolson method to include a correction term that takes into account the radial symmetry of your domain. This additional term could help ensure that the mass is conserved throughout the solution.

Another approach could be to use a different numerical method that is specifically designed for radially symmetric problems. For example, the Alternating Direction Implicit (ADI) method or the Fractional Step Method can both be used for solving parabolic PDEs in radially symmetric domains and may provide better mass conservation.

It is also important to make sure that your implementation is correct and that there are no errors in the code. You could consider double-checking your code or seeking feedback from other researchers in the field.

I hope these suggestions are helpful and I wish you success in solving your problem. Keep up the good work!