I'm trying to solve 1D Euler equations with gravity in spherical coordinates using a finite-difference TVD MacCormack method on a non-uniform grid of $N$ components, following the method provided in Appendix A in this paper. My aim is to evaluate stationary solutions.

The discretization for the source term in the second equation should guarantee that hydrostatic equilibrium, if reached, is maintained. I am using fixed BC at the lower boundary for all phisical variables and free-flow BC at the upper boundary; since the TVD correction requires 2 cells on both side of each cell, BC have to be imposed using two ghost cells:

\begin{equation} \begin{cases} \rho_1=\rho_2=1 \\\ u_1=u_2=0 \\\ p_1=p_2=1 \end{cases} ,\qquad \begin{cases} \rho_{N}=\rho_{N-1}=\rho_{N-2} \\\ u_{N}=u_{N-1}=u_{N-2} \\\ p_{N}=p_{N-1}=p_{N-2} \end{cases} \end{equation}

I was testing the code by using as IC a particular equilibrium solution, i.e.:

\begin{equation} \begin{cases} \rho_{st}(r)=e^{-r}\\\ v_{st}(r)=0\\\ p_{st}(r)=(1+r)e^{-r}\\\ \end{cases}, \end{equation}

with $g=-\frac{\partial \phi}{\partial r}$, $\phi=\frac{r^2}{2}$ and $Q_\text{net}=0$. The code works well for all middle-grid points except near the boundary; in fact, if we consider the predictor step for the $(N-2)$-th component of the second equation at the first time step we have

\begin{equation} \begin{aligned} \big(\rho u & r^2\big)_{N-2}^{1/2} =\big(\rho u r^2\big)_{N-2}-\lambda \Big[\big(\rho u^2 r^2\big)_{N-1}-\big(\rho u^2 r^2\big)_{N-2}\Big]\\\ &-\frac{\lambda}{2}\Big[ \Big(\big(\rho r^2\big)_{N-1}+\big(\rho r^2\big)_{N-2}\Big)\Big(\phi _{N-1}-\phi _{N-2}\Big)+\big(r^2_{N-1}+r^2_{N-2}\big)\big(p_{N-1}-p_{N-2}\big)\Big] \end{aligned} \end{equation}

where $\lambda=\frac{\Delta t}{r_{N-1}-r_{N-2}}$.

In the first time step, since $v(r)=0$ and given the BC for pressure at the outer boundary, the equation reduces to

$$ \big(\rho u r^2\big)_{N-2}^{1/2}=\Big(\big(\rho r^2\big)_{N-1}+\big(\rho r^2\big)_{N-2}\Big)\Big(\phi _{N-1}-\phi _{N-2}\Big) $$

This clearly shows that the equilibrium at this grid point is not preserved. The same thing happens near the lower boundary in the corrector step for index $3$. If I run the code for more time steps, the error propagates throughout the interior and the solution deviates from equilibrium overall points.

Is this due to my choice of BC? Or it's a more general problem of the reconstruction of MacCormack method? Do you think that using WENO schemes should be better to achieve better stationary solution (I'll have to modify the equations in the future like in the previously mentioned paper)?

If you have any reference for useful books or papers, it would be appreciated!

  • 1
    $\begingroup$ Hi, Welcome to SC SE! When you say upper and lower boundary, do you mean r = 0, and r = R? Also, the bc you described isn't necessarily great, I've never seen someone use the B.C. you specified, and I'm not sure of the possible impacts. $\endgroup$ – EMP Aug 8 '19 at 19:36
  • $\begingroup$ Yes of course I mean boundaries at $r=1$ and $r=R$. This kind of boundary condition are the ones used in the article I mentioned (see par. 2.3) and they are suitable for the physical phenomenon that is being modeled $\endgroup$ – Andrea Caldiroli Aug 8 '19 at 20:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.