I wish to solve the following using the finite volume method: $$\frac{\partial u}{\partial t}=\frac{D}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u}{\partial r}\right)+Q(t,r)$$ with the following boundary conditions: $$\frac{\partial u}{\partial r}\Bigg|_{r=0}=0,\quad -D\frac{\partial u}{\partial r}\Bigg|_{r=R}=hT(t,R)$$ Integrating the PDE (multiplied with $r$) from $r_{j-\frac{1}{2}}$ to $r_{j+\frac{1}{2}}$ to get: $$\frac{d}{dt}\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rudr=D\left[ r\frac{\partial u}{\partial r}\right]_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}+\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rQ(t,r)dr$$ I denote for general $f(t, r) $: $$\tilde{f}_{j}(t)=\frac{1}{\delta r}\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rf(t,r)dr$$ which makes my equation: $$\frac{d\tilde{u}_{j}}{dt}=\frac{D}{\delta r}\left[ r\frac{\partial u}{\partial r}\right]_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}+\tilde{Q}_{j}$$ Now is this a valid approximation? $$r\frac{\partial u}{\partial r}\Bigg|_{r=r_{j+\frac{1}{2}}}=r_{j+\frac{1}{2}}\cdot\frac{\tilde{u}_{j+1}-\tilde{u}_{j}}{\delta r}$$ and likewise for the other flux. I'm also having issues with implementing the inner boundary condition at $r=0$ which is confusing me a lot.

  • $\begingroup$ Apparently the formula for $r \partial_r{u}$ is meant to be for $j+1/2$, otherwise it is fine, what is confusing there? For the BC at $r$=0, if it is Dirichlet type, e.g., $u$=1, then that should be no problem. For Neumann $\partial_r u =0$ should also be easy. Probably you would not want a finite (nonzero) radial derivative at $r$=0, that would make your function non-differentiable there, which would be physically possible only for a $\delta$-function source at the origin. For a more complex BC at $r$=0, it would be fine to set it at some finite (but small) distance from the origin. $\endgroup$ Apr 1, 2022 at 21:29
  • $\begingroup$ Agreed with @MaximUmansky this all seems fine (besides the typo of $j+1/2$). $\endgroup$
    – user20857
    Apr 2, 2022 at 19:31
  • $\begingroup$ I'm confused about how to do the outer boundary condition. How do I link the value of u on the outer boundary to $\tilde{u}$ on the centre of the last cell? $\endgroup$ Apr 20, 2022 at 11:33
  • $\begingroup$ Store your unknowns at cell centers. This means you will not have $u$ at $r=0$ or $r=R$. You will implement the boundary conditions as fluxes. If your first unknown is $u_1$, then it will be located at $\delta r/2$. $\endgroup$
    – cfdlab
    Nov 28, 2023 at 13:37
  • $\begingroup$ Your definition of $\tilde f_j$ is not correct, even dimension-wise. You need to divide by $r_j$. $\endgroup$
    – cfdlab
    Nov 28, 2023 at 13:51

1 Answer 1


The simple answer is, as I now realise, that you are free to compute the fluxes however you please, and the fluxes are evaluated on the cell boundaries, so I can just insert the actual values for the flux and both the inner and outer boundaries.

  • $\begingroup$ Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Nov 29, 2023 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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