7

I want to discretize the wave equation

$$\frac{1}{c^2}\frac{\partial^2\psi\left(\vec{r},t\right)}{\partial t^2}=\triangle\psi\left(\vec{r},t\right)$$

in polar coordinates. I find the following relations:

$$\vec{r}=r\mathbf{e}_r$$ $$\triangle=\frac{1}{r}\partial_r r \partial_r+\frac{1}{r^2}\partial_\varphi^2$$

Define:

$$r=i\Delta r,\qquad\varphi=j\Delta\varphi,\qquad t=n\Delta t$$

Thus:

$$\psi\left(r,\varphi,t\right)=\psi\left(i\Delta r,j\Delta \varphi,n\Delta t\right)\equiv\psi_{i,j}^{\left(n\right)}$$

Using the forward difference in $r$ and central differences in $\varphi$ and $t$, I get:

$$\psi_{i,j}^{\left(n+1\right)}=2\psi_{i,j}^{\left(n\right)}-\psi_{i,j}^{\left(n-1\right)}+\frac{c^2 \Delta t^2}{\Delta r^2}\left[\psi_{i+2,j}^{\left(n\right)}+\frac{1}{i}\psi_{i+2,j}^{\left(n\right)}-2\psi_{i+1,j}^{\left(n\right)}-\frac{1}{i}\psi_{i+1,j}^{\left(n\right)}+\psi_{i,j}^{\left(n\right)}\right]+\frac{c^2\Delta t^2}{i^2\Delta r^2\Delta \varphi^2}\left[\psi_{i,j+1}^{\left(n\right)}+\psi_{i,j-1}^{\left(n\right)}-2\psi_{i,j}^{\left(n\right)}\right]$$

Now this is divergent at $r=0$, did I something wrong?

And how can I set restrictions to the steps $\Delta r$, $\Delta \varphi$ and $\Delta t$?

Andy
  • 393

1 Answers1

8

If you look at the Laplacian: $$ \nabla^2=\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right)+\frac1{r^2}\frac{\partial^2}{\partial\phi^2} $$ you can clearly see that this diverges at $r=0$ so discretization of this should also diverge.

There are three solutions to remedying the divergent feature that I can think of:

  1. Choose a problem such that you don't need to worry about $r=0$ (e.g., $r_{min}=0.1$)
  2. Use cell-centered values for $r$ (e.g., $r_{min}=\frac12dr$)
  3. Use two meshes a Cartesian one for $r<1$ and a polar one for $r\geq1$, merging them at $r\simeq1$ (requires at least 1 cell of overlap).

I've never used option #3, so I cannot speak anything about it. Option #1 works in certain cases, but not for all. Option #2 I have worked with and have seen in other hydrodynamics codes that I've used, so it might be your best bet.

As far as your discretization goes, I think it might be more clear to use $r_i$ and $\phi_j$ for the positions and $dr$ and $d\phi$ for the cell widths. Using this, $$ \frac{1}{r}\,\frac{\partial}{\partial r}\left(r\frac{\partial\psi}{\partial r}\right)=\frac{1}{r_i}\left(r_{i+1/2}\frac{\psi^n_{i+1,j}-\psi^n_{i,j}}{dr^2}-r_{i-1/2}\frac{\psi^n_{i,j}-\psi^n_{i-1,j}}{dr^2}\right) \\ \frac{1}{r^2}\frac{\partial^2\psi}{\partial\phi^2}=\frac{1}{r_i^2}\left(\frac{\psi^n_{i,j+1}-2\psi^n_{i,j}+\psi^n_{i,j-1}}{d\phi^2}\right) $$ where $r_{i+1/2}=\frac12\left(r_{i+1}+r_i\right)$.

Your choice for $dr$ and $d\phi$ are up to you. It might be easier to have $dr=d\phi$, but it might not be worth enforcing that behavior. The timestep, however, is limited by the Courant-Friedrichs-Levy condition. For the wave equation, your requirement is such that $$ c^2\frac{dt}{\min(dr^2,d\phi^2)}\leq\frac12 $$ Formally, the constant can be 1 instead of $\frac12$, but for multidimensional systems, it is appropriate to use $1/n_{dim}$. For time-stepping purposes, this is going to be the slow point of your code. It might be worth looking into implicit methods instead of the explicit method you are currently working with.

Kyle Kanos
  • 29,127