Consider the spherical advection problem: describing the conservation of a property $u$ in a closed spherical domain.

$$ \frac{\partial u}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)=0 $$

In the domain $0<r<1,\ t\geq0$, with the following initial/boundary conditions:

$$ \left.\frac{\partial u}{\partial r}\right|_{r=a}=0, \\ u(r,t=0)=u_i $$ Where $\boldsymbol{v}=\boldsymbol{v}(r,t)$ and $u_i$ are known. ($\boldsymbol{v}\leq0$). Other posts (here, and here) suggested discretizing the conservative form. Forward explicit (downwind) discretization leads to: $$ r\to i\Delta r,\ t\to n\Delta t,\ u(r,t)\to u_i^n \\ \frac{\partial u}{\partial t} \to \frac{u^{n+1}_i-u^{n}_i}{\Delta t} \\ \frac{\partial(r^2u\ \boldsymbol{v})}{\partial r} \to \frac{ ((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n- (i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n }{\Delta r} $$ And the equation becomes: $$ \frac{u^{n+1}_i-u^{n}_i}{\Delta t}= -\frac{1}{(i\Delta r)^2} \frac{ ((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n- (i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n }{\Delta r} $$

My problems are:

  • I cannot use this equation at $r=1$ (last node) because i have no boundary condition to fix the ghost node at $r=1+\Delta r$.
  • I'm unsure how to write the equation at $r=0$. I can either use L'Hôpital to write $$ \frac{\partial u}{\partial t}+ \lim_{r\to0}\left(\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)\right) = \frac{\partial u}{\partial t}+ \lim_{r\to0}\left(\frac{1}{2r}\frac{\partial^2 }{\partial r^2}\left(r^2u\ \boldsymbol{v} \right)\right) =\\ \frac{\partial u}{\partial t}+ \frac{1}{2}\frac{\partial^3 }{\partial r^3}\left(r^2u\ \boldsymbol{v} \right) =0 $$ And discretize the above or I can use the forward derivative and write $$ \frac{u_1^{n}-u_0^{n}}{\Delta r}=0\therefore u_0^{n}=u_1^{n} $$

What is the best way to discretize the problem and how do I address the first and last nodes?


I found another post (here) with almost the same problem. One of the answers suggested an interpolation for the stencil at the last node, which i tried with unsuccessful results. The solution returns values greater than $1$ at $r=1$ and does not converge at all as the nodes increase (unstable). The image below shows the solutions for 21, 51 and 101 spatial nodes. The CFL number is below $1$. Solution should stabilize at $0.91$, and be less than $1$ for every instant.

results for 21/51/101 nodes

  • $\begingroup$ In case anyone wants to look further into the origin of this problem, this is the article I am trying to solve numerically. $\endgroup$
    – eg.Zeta
    Aug 24, 2016 at 11:28
  • 1
    $\begingroup$ I'm confused. With $\nu<0$, you have material flowing inward from the outside boundary, so your discretization is upwind-biased (not downwind), and you need a boundary condition at the outside of the sphere, which you haven't provided. Furthermore, it seems intuitively like the exact solution should develop a singularity at the origin, since all the mass arrives there eventually. $\endgroup$ Aug 25, 2016 at 10:58


Your Answer

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