The time step constraint in radiation hydrodynamics

ABSTRACT

Explicit radiation hydrodynamic simulations of the atmospheres of massive stars and of convection in accretion discs around white dwarfs suffer from prohibitively short time steps due to radiation. This constraint is related to the cooling time rather than the radiative pressure, which also becomes important in hot stars and discs. We show that the radiative time step constraint is governed by the minimum of the sum of the optically thick and thin contributions rather than the smaller one of the two. In simulations with the Pencil Code, their weighting fractions are found empirically. In three-dimensional convective accretion disc simulations, the Deardorff term is found to be the main contributor to the enthalpy flux rather than the superadiabatic gradient. We conclude with a discussion of how the radiative time step problem could be mitigated in certain types of investigations.

KEYWORDS:

• Numerical stability
• radiative cooling
• convection
• accretion discs
• hot stars

Acknowledgments

This paper is dedicated to Ed Spiegel. If it was not for the radiative time step problem, the Brandenburg and Spiegel (2006) paper would have been published by now! We thank Matthias Steffen for sharing with us his experience and knowledge regarding the time step constraint in stellar convection simulations. We are also grateful to the three referees for their thoughtful comments. AB also acknowledges Fazeleh (Sepideh) Khajenabi for her work on the radiative time step problem while visiting Nordita in the spring of 2010. UD thanks Akshay Bhatnagar for useful discussions. Simulations presented in this work have been performed with computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm. The source code used for the simulations of this study, the Pencil Code, is freely available on https://github.com/pencil-code/. The DOI of the code is http://doi.org/10.5281/zenodo.2315093. The setups of runs and corresponding data are freely available on https://www.nordita.org/∼brandenb/projects/tstep/.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Axel Brandenburg http://orcid.org/0000-0002-7304-021X

Upasana Das http://orcid.org/0000-0003-2302-2280

Correction Statement

This article has been republished with a minor change. This change does not impact the academic content of the article.

Notes

1 See http://www.astro.uu.se/∼bf/co5bold_main.html.

2 For a second order discretisation, for example, we have ${\mathrm{\nabla }}^2{f}_i=(f_{i+1}-2f_i+f_{i-1})/\delta x^2$ in one dimension, but ${\mathrm{\nabla }}^2{f}_{ij}=(f_{i+1j}+f_{ij+1}-4f_{ij}+f_{i-1j}+f_{ij-1})/\delta x^2$ in two dimensions, so at the centre point, $f_{ij}$, the coefficient increases from $2/\delta x^2$ to $4/\delta x^2$ in two dimensions, and to $6/\delta x^2$ in three dimensions. This applies analogously also to the sixth order discretisation used in the Pencil Code, except that the coefficient now increases to $49/18\approx 2.72$ instead of 2 per direction as for the second order case. Also, for a given function f, the value of $-\left({\mathrm{\nabla }}^{2}f\right)/f$ depends on the numerical scheme. For a checkerboard pattern of f (e.g. an alternating sequence with $-1$, $+1$, etc, in one spatial dimension), using a second order scheme, the value is $-4/\delta z^2$ per direction, while the analytic value is $-{\pi }^2/\delta z^2\equiv k_{\mathrm{N}\mathrm{y}}^2$, which is more than twice as much.

Additional information

Funding

This research was supported in part by the National Science Foundation under the Astronomy and Astrophysics (Division of Astronomical Sciences) Grants Program (grant 1615100), and the University of Colorado through its support of the George Ellery Hale visiting faculty appointment.