I've been trying unsuccessfully to solve the 1D diffusion equation in spherical co-ordinates using finite differences. For the life of me, I cannot get the thing to work. Is this a known phenomena? What is the usual method used for this seemingly simple task?
What do you mean by it "not working" ? Numerically unstable? Unphysical results? Does it work away from the origin?
Hi Mat Hunt
Long time, no see! Hope all is well.
I'm presuming that you're trying to solve T_t = T_{rr} + (2/r)T_r ? Perhaps with boundary conditions such as T_r=0 at r=0 and maybe T=const at r=1, say? Actually, that solution is T=const, and therefore I am presuming that you're doing something a little different from what I have suggested. Please let me know.
Andrew
Possibly, you should apply your FDM scheme in spherical coordinates. It is not like FDM will not give results. You can get some good research articles where researchers have solved the diffusion equation in spherical coordinates using FDM.
D. Andrew S. Rees, that's the equation I'm solving. I have zero flux at r=0, which is what people use and I have a Neumann boundary condition on the outer boundary. I'm essentially evaluating the equation at the point (t_{i+0.5},r_j) and obtaining a stencil that way via the Crank-Nicholson method.
Loïc O. Le Cunff By not working I mean that the results depend strongly on how fine the mesh is, giving completely different results depending upon how many points you use for time and distance meshes. The other thing is that when you compute the average quantity withing the domain, it's proportional to the flux on the boundary, this should be reproduced in the code, but it isn't when computed from the code.
Bharat Soni I have been doing this but it's remarkable sparse.
Possible...I saw some articles they were using solving diffusion equations in spherical polar coordinate systems.
Hi Mat Hunt
Theres been no significant comment since you opened this discussion. It just struck me that you didn't mention what the Neumann condition is at the outer boundary. Then I thought, pehaps, I need to specify something such as T_r=C, a constant, and solve for the steady-state solution. It is possible only if C=0. Therefore the long-term fate of a given r-dependent initial profile is that the temperature eventually evolves to being a constant.
Hi Mat Hunt if the attached file is yours I think the 10 eq. is wrong.
it has to be something like this:
3a * Ti+1,1 - (1+3a) * Ti+1,0 = -3a Ti,0 - (1-3a) Ti,0
- Badges
- Science topic
- Similar topics
- Mathematics