Numerical Solution of transient heat equation with source in cylindrical co-ordinates

39 Views Asked by At

I am attempting to use the finite difference scheme on a uniform grid in $r$. I have a no flux boundary condition at $r=0$ which is proving to be a pain to implement. I'm having difficulty in doing this, I have checked and rechecked my code and derivation of the code but I cannot see an error. I am using the Crank-Nicholson method as I'm told that nice and stable and although I get something nice, physically it's wrong. The stencil I get is: $$2\alpha u_{i+1,1}-(2\alpha+1)u_{i+1,0}=-2\alpha u_{i,1}+(2\alpha-1)u_{i,0}-\frac{1}{2}(Q_{i+1,0}+Q_{i,0})\delta t$$ This doesn't work and I don't know what's going on.

Do I need to abandon this completely or is there something I can do to rescue this?