A two-column formalism for time-dependent modelling of stellar convection

# A two-column formalism for time-dependent modelling of stellar convection

I. Description of the method
Alexander Stökl CRAL, Université de Lyon, CNRS (UMR5574), École Normale Supérieure de Lyon, F-69007 Lyon, France
11email: alexander.stoekl@ens-lyon.fr
###### Key Words.:
hydrodynamics — Cepheids — convection — methods: numerical
###### Abstract

Context:In spite of all the advances in multi-dimensional hydrodynamics, investigations of stellar evolution and stellar pulsations still depend on one-dimensional computations. This paper devises an alternative to the mixing-length theory or turbulence models usually adopted in modelling convective transport in such studies.

Aims:The present work attempts to develop a time-dependent description of convection, which reflects the essential physics of convection and that is only moderately dependent on numerical parameters and far less time consuming than existing multi-dimensional hydrodynamics computations.

Methods:Assuming that the most extensive convective patterns generate the majority of convective transport, the convective velocity field is described using two parallel, radial columns to represent up- and downstream flows. Horizontal exchange, in the form of fluid flow and radiation, over their connecting interface couples the two columns and allows a simple circulating motion. The main parameters of this convective description have straightforward geometrical meanings, namely the diameter of the columns (corresponding to the size of the convective cells) and the ratio of the cross-section between up- and downdrafts. For this geometrical setup, the time-dependent solution of the equations of radiation hydrodynamics is computed from an implicit scheme that has the advantage of being unaffected by the Courant-Friedrichs-Lewy time-step limit. This implementation is part of the TAPIR-Code (short for The adaptive, implicit RHD-Code).

Results:To demonstrate the approach, results for convection zones in Cepheids are presented. The convective energy transport and convective velocities agree with expectations for Cepheids and the scheme reproduces both the kinetic energy flux and convective overshoot. A study of the parameter influence shows that the type of solution derived for these stars is in fact fairly robust with respect to the constitutive numerical parameters.

Conclusions:

## 1 Introduction

Convection is one of the persistent problems in stellar astrophysics. Almost all stars contain regions where convective transport is important; in the photosphere, in the envelope, or in the interior where nuclear burning occurs. A description of convection is therefore an essential ingredient to all types of investigations of stellar structure and evolution. Unfortunately, due to its nonlinear, nonlocal, and multi-length-scale nature, modelling convection turns out to be an intricate problem.

Conceptually, there are several different approaches to the numerical simulation of convective transport in stars. The most straightforward approach is the time-dependent solution of the equations of radiation hydrodynamics in a 3D (or at least 2D) domain to compute the convective flow patterns directly. However, this process is very expensive in computing time, prohibitively so for some applications. This particularly applies to problems with a large difference in relevant timescales, for instance between the thermal and acoustic timescale in the stellar interior, or the hydrodynamics and radiative timescale in the outer layers of luminous stars. An additional limit is imposed by the restricted spatial resolution. Even the most elaborate, high-resolution simulations of stellar convection are only capable of resolving the largest scales in the convective velocity field. The effect of turbulence on smaller length scales is effectively ignored, even though it is a possible interpretation to attribute the intrinsic numerical dissipation of the scheme to some (unknown) unresolved turbulence. In particular, hydrodynamics codes often include artificial viscosity for numerical stability, and ‘unresolved turbulence’ is the only physical mechanism that could be used to justify the inclusion of this additional dissipation.

Despite the poor description of turbulence in hydrodynamics computations, multi-dimensional simulations of solar granulation (e.g. Stein & Nordlund, 1998; Asplund et al., 2000; Wedemeyer et al., 2004) achieve remarkable quantitative agreement with observations.

A completely different and more subtle approach than trying to resolve turbulence on large numerical grids are convection models that use an equation, or a set of equations, to describe convective transport either with a heuristic parametrization or based on turbulence theory.

The most widely used of these 1D-descriptions is the well known mixing-length theory (MLT) (Böhm-Vitense, 1958; Cox & Giuli, 1968), which originates in ideas of Prandtl (1925). The MLT has been remarkably successful in stellar astrophysics in application to stars from white dwarfs to super giants, probably because of its simple yet flexible parametrization and its robust reference to the adiabatic temperature gradient.

Modern alternatives to the MLT are convection models in which the original single-eddy assumption of the MLT has been replaced with a full spectrum of turbulence (e.g. Canuto & Mazzitelli, 1991; Canuto, 1996; Canuto et al., 1996) by either assuming or computing a turbulent energy spectrum. For some types of stars, these models can also avoid the mixing-length scale as a free parameter. Convection models of this type are included in many stellar evolution codes as an alternative to MLT.

However, all of these models are, in a similar way to the original MLT, local theories that do not provide any information about overshoot. In stellar-evolution codes, this deficit is overcome by adding overshoot by means of a separate (typically diffusive) parametrization.

Another type of one-equation models consists of a time-dependent equation for the turbulent kinetic energy (e.g. Stellingwerf, 1982; Kuhfuß, 1986; Gehmeyr & Winkler, 1992) using heuristic approximations for individual terms. Containing a diffusion term for the turbulent kinetic energy, they also allow a simple form of non-locality. These convection models are mainly geared towards computations of stellar pulsations (e.g. Bono & Stellingwerf, 1994; Feuchtinger, 1999b; Kolláth et al., 2002) where one is interested in the time-dependence of the convective transport.

Finally, the most complete way to model turbulence is the ‘Reynolds stress approach’ by solving a set of moment equations (Canuto, 1992, 1993, 1997; Xiong, 1989; Xiong et al., 1997) that are terminated by a closure model at the third or fourth order. These closures are based on either the quasi-normal approximation for the fourth-order moments or use a parametrization with reference to measured data, hydrodynamics (large-eddy) simulations, or concepts such as the ‘plume model’ (see below). Turbulence models of this type are able to describe convective transport in a time-dependent and non-local way. For applications of the Reynolds stress model to stellar surface convection zones, see Kupka (1999); Kupka & Montgomery (2002); Montgomery & Kupka (2002).

The two-column scheme, presented in this paper, is in-between the above categories, combining a hydrodynamics simulation of the convective fluid flow with a parametrized, predefined geometry of the flow patterns. This setup is almost as simple as a 1D description; it describes up- and downstream with two parallel radial columns, and fluid flow over an interface between those two columns allows a basic circulating convective motion. The two-column model could therefore be regarded as a simplistic 2D hydrodynamics scheme, limiting the horizontal range of the grid to just two cells. The very coarse description of the convective velocity field effectively implies that the most extensive convective flow patterns generate the majority of the convective transport. However, this assumption does not differ significantly from what is assumed in multi-dimensional hydrodynamics computations that are also unable to resolve the full spectrum of convective turbulence. Since multi-dimensional hydrodynamics achieve, in spite of this limitation, a good agreement with observations, the scenario of macroscopic convective patterns with distinct up- and downstream regions and little sub-structure (as also observed in the solar granulation) appears to be a sufficient description of the actual physics (Nordlund et al., 1997). The two-column approach may therefore also give reasonable results.

In the two-column scheme, the convective flux is computed directly from hydrodynamics and not from a heuristic model. Although it is not without numerical parameters, there is no ‘mixing-length parameter’, nor anything equivalent. The method is also intrinsically non-local and the thickness of convective regions and the amount of overshoot are obtained consistently.

The basic idea of modelling convection using separate radial stratifications for up- and downstream regions is in fact an old one. In the 1960’s to 1970’s, predating advances in computing power that enabled 2D and 3D hydrodynamics computations to become possible, several similar two- or multi-stream models were devised. From observations, the solar granulation pattern appeared to be separable to almost distinct hot and cool areas; it was therefore a logical first step to place two stratifications next to each other in order to construct more realistic models of the solar photosphere (Margrave & Swihart, 1969; Nordlund, 1976). More recently a two-stream model was applied by Lesaffre et al. (2005) to investigate the convective Urca process in supernova-progenitor white dwarfs. In geophysics, a similar concept, known as ‘plume model’, was introduced by Morton et al. (1956). In its basic form, this model considers only plumes that are immersed in a static, surrounding medium, although there are also models that consider both up- and downdrafts and their interaction (e.g. Telford, 1970; Wang & Albrecht, 1986; Chatfield & Brost, 1987; Randall et al., 1992). The idea of separated up- and downwards streams was also used to construct closures for turbulence models (Abdella & McFarlane, 1997; Zilitinkevich et al., 1999; Lappen & Randall, 2001; Gryanik & Hartmann, 2002; Gryanik et al., 2005; Canuto et al., 2007)

However, all existing multi-stream models differ from the present attempt in that the ‘two-column-scheme’ introduced below is based on fully implicit time-dependent radiation hydrodynamics in both radial and horizontal directions without any ad-hoc assumptions or parametrizations for the physical coupling of the two columns. The term ‘two-column’ (in contrast to ‘two-stream’) was chosen intentionally because the discretization scheme resembles that of two 1D discretizations placed beside each other in two parallel columns. In analogy to ‘2D’ for ‘two-dimensional’, we will occasionally refer to ‘two-column’ as ‘2C’ in the following sections.

The remaining paper is structured as follows. The next section, Sect. 2, introduces the two-column discretization scheme and its geometric derivation. The equations of radiation hydrodynamics are given in Sect. 3 in analytical form, while Sect. 4 describes their discretization and considers details such as artificial viscosity and radiative transport. Section 5 presents the deployed solution algorithm, followed by a demonstrating example and some parameter studies in Sect. 6. Finally, Sect. 7 draws the conclusions and summarizes the paper. A verification of the method by comparison with detailed 2D hydrodynamics computations as well as applications in time-dependent calculations of Cepheids’ pulsations will be given in the forthcoming part II paper of this series.

## 2 The two-column discretization scheme

Figure 1 shows the setup of the two-column discretization scheme and the localization of the primary variables (for a complete listing of the primary variables, see Table 1). The columns do not correspond directly to an individual convective cell but should be considered as a representation of all up- and downstream flows, respectively. Correspondingly, the interface between the two columns represents the sum of all contact surfaces between up- and downdrafts.

In the 2C-scheme, the horizontal components of fluid flow and radiation, which actually occur in both the - and -direction of spherical geometry, are described each by just one ‘horizontal’ variable. For these variables, and respectively (see Table 1), the subscript ‘’ does not refer to spherical geometry components but is used more generally to denote ‘horizontal’ variables.

The geometrical configuration of the discretization scheme is specified by two parameters. The first parameter is the typical horizontal length scale , which can be interpreted as the diameter of the convective cells or the typical horizontal distance between up- and downdrafts. In contrast to the sketch in Fig. 1, the two columns, in general, do not have the same size. The second parameter ( for column fraction) specifies the fraction of the sphere that is associated with column 1, and correspondingly, with , is that allocated to column 2.

These two parameters, and , with their straightforward geometrical meaning are the main free parameters of the convection model. In Sect. 2.4, we will introduce a third constitutive parameter correlated with horizontal advection. Other numerical parameters, such as solution accuracy, radial advection scheme, artificial viscosity, boundary conditions, and grid resolution have only a minor effect on the solution.

An equivalent yet more concrete quantity than is the number of convective cells on a sphere . In principle, it is possible to assign different values of to individual shells, or to specify an analytical relation that defines , for instance as a function of the radius. However, in the absence of a robust physical indication for the behavior of , the present implementation uses the same for every shell independent of the radius. That way, the up- and downdrafts are assumed to retain their identity throughout the convective region, which appears reasonable for photospheric convection zones of moderate depth. The relative cross-sections and must remain constant in all cases since changing their values would tilt the interface between the columns from the radial (coordinate) direction.

By definition, and are related by

 D=2r√2/N, (1)

which is based on the assumption of circular convective cells as sketched in Fig. 2. Since the formalism of the two columns is symmetric apart from in the relative cross-sections, it makes no difference whether these circular cells are considered as up- or as downdrafts. The right hand side of Fig. 2 illustrates the computation of the interface area between the two columns. Using Eq. 1 and summing for columns, we obtain

 Aintf=π2√N/2(r2i+1−r2i). (2)

Despite the geometric motivation given in Fig. 2, this picture should not be interpreted literally. By combining both the - and -directions of spherical geometry to one generic horizontal variability, the direct correlation with the three-dimensional configuration disappears. It is therefore not sensible to interpret expressions from the two-column formalism using, for example, a specific slice through the setup shown in Fig. 2.

However, the only points where the geometrical configuration enters the scheme are the definitions of the interface area in Eq. 2 and of the horizontal derivatives in Eqs. 3 & 4. Since these equations are coupled with each other by the requirement to ensure that Gauß’s theorem applies also in the discrete case, a different geometrical picture would change Eqs. 2 – 4 by only a constant factor.

#### Horizontal derivatives

Using the typical horizontal length scale , we can approximate derivatives in the horizontal direction by

 1r∂∂θX≃1DΔθ(X)=12r√N/2Δθ(X), (3)

where ‘’ is again used to denote the generic horizontal extension of the 2C-scheme.

For second-order derivatives with respect to , we adopt the estimate

 (4)

which is based on the assumption of a basically periodic variation of in -direction (and consequently of as well) due to the alternating succession of up- and downdraft columns. This already requires over-stretching of the geometrical picture but since these derivatives exist at a less important point (in the term of the viscous forces, Eqs. 40 & 44), the approximative evaluation of is acceptable.

### 2.1 Scalar discretization

The discretization of scalar physical variables and equations uses discretization volumes similar to that highlighted in gray in Fig. 1. The volume of the scalar cells (distinguished by the prefix ‘S’ for scalar) is computed to be the appropriate fraction of the shell between the radii and

 S_vol1=cf14π3(r3i+1−r3i)S_vol2=cf24π3(r3i+1−r3i). (5)

The advective fluxes (transported ‘volume’ during a time step) for the scalar discretization are indicated by arrows in Fig. 1. Radial advection consists of two contributions; one from the proper motion of the fluid, and one due to movement of the adaptive grid

 S_flux1=cf1[4πr2iu1,iδt−4π3(rnewi3−roldi3)]S_flux2=cf2[4πr2iu2,iδt−4π3(rnewi3−roldi3)] (6)

where is the time step during which the grid adaptivity alters the radius of the grid point from to . Note that the individual radial velocities and were used in the two columns.

The horizontal advective flux between the two columns is computed to be

 S_fluxθ=Aintfuθ,iδt=π2√N/2(r2i+1−r2i)uθ,iδt. (7)

For both the horizontal velocity and the horizontal advective flux , a positive sign corresponds to a fluid flow from column 1 to column 2 by convention. The analogous convention is also used for the horizontal radiative flux .

### 2.2 Vector discretization – radial

Vector-type variables and equations are discretized on a staggered mesh where the radial part closely resembles that given by Dorfi et al. (2006). To define discretization volumes for the radial vector components, we start by defining ‘averaged’ radii located in between the grid point positions and denoted by

 ¯¯¯r3i≡r3i+12=12(r3i+r3i+1). (8)

Figure 3 shows two of these averaged radii, and , which establish the discretization volumes centered around the grid point .

Using the definition of , we can now compute the volumes (with the prefix ‘V’ for vector) of these cells

 V_vol1=cf14π312(r3i+1−r3i−1)V_vol2=cf24π312(r3i+1−r3i−1). (9)

For the advective flux in the radial direction, we interpolate the velocity assuming flux conservation, i.e.

 r2i+1/2ui+1/2=12(r2iui+r2i+1ui+1), (10)

and therefore in analogy with Eq. 6, we obtain

 V_flux1=cf1[4π12(r2iu1,i+r2i+1u1,i+1)δt−4π3(¯¯¯rnewi3−¯¯¯roldi3)]V_flux2=cf2[4π12(r2iu2,i+r2i+1u2,i+1)δt−4π3(¯¯¯rnewi3−¯¯¯roldi3)]. (11)

As indicated by the arrows in Fig. 3, the horizontal flux between the two radial vector volumes and is composed of two parts correlated with and

 V_fluxθ=π2√N/2[(¯¯¯r2i−r2i)uθ,i+(r2i−¯¯¯r2i−1)uθ,i−1]δt. (12)

### 2.3 Vector discretization – horizontal

The discretization of the horizontal components of vector variables and equations (namely of and ) uses discretization volumes as illustrated in Fig. 4. The corresponding volumes and fluxes are labeled with the prefix ‘H’ for horizontal. The discretization cell with is centered on the interface between the columns and considers half the volume of the sphere

 H_vol=124π3(r3i+1−r3i). (13)

The second half of the volume is assumed to mirror the physics of the first one. Even though the discretization volume represents only one half of the shell, it therefore describes the horizontal components of the entire sphere.

From comparison with Eq. 5, we observe that ; the flux in the radial direction is assembled in the same way from and (see Eq. 6)

 H_flux=12[4πr2i(cf1u1,i+cf2u2,i)δt−4π3(rnewi3−roldi3)]. (14)

Since the discretization scheme represents a large number of convective cells distributed over the sphere, the two columns can be considered as part of a sequence of alternating up- and downdrafts. Along this sequence, the direction of horizontal fluid flow and radiation switches its sign repeatedly. This implies that a right-hand orientated flow (as illustrated in the lower part of Fig. 4) is confronted with an equal flow in the opposing direction when reaching the (in this case) right hand cell boundary. To allow for this effect, a dissipation term was included in the equation of motion that could be interpreted as annihilation of the momentum of the two opposing flows

 Fanhl=−|S_fluxH|(cf1ρ1+cf2ρ2)uθ (15)

where is the momentum in the cell and describes the volume fraction swept against the horizontal boundary. In the equation of internal energy, this dissipated energy enters as

 Eanhl=|S_fluxH|(cf1ρ1+cf2ρ2)u2θ. (16)

The location at which this energy is deposited depends on the direction of horizontal flow. In the case sketched in Fig. 4, the dissipated energy would be deposited into column 2.

As we will make use of it in the set of discrete equations (Table 3), we finally define , the averaged density appropriate for the ‘horizontal’ discretization volume

 ¯¯¯ρθ=(cf1ρ1+cf2ρ2). (17)

Despite the similar notation, this horizontally averaged density should not be mistaken for the radially averaged densities and that are introduced in Sect. 4.5.

### 2.4 Horizontal advection and energy conservation

Advection in the radial direction is considered using a second-order van Leer advection scheme (van Leer, 1974, 1977), but for horizontal advection, i.e. fluid flow from one column to the other, a higher order scheme is obviously not applicable. The most straightforward approach is donor cell advection, but, in general, we can allow some variation within the columns, as illustrated in Fig. 5 for the example of density. This produces the following advection scheme

 ˜ρ={(1−λ)ρ1+λρ2for Col.~{}1→Col.~{}2(1−λ)ρ2+λρ1for Col.~{}1←Col.~{}2, (18)

Despite this potentially unphysical behavior for larger values of , the formalism in Eq. 18 provides an additional free parameter for adjusting the convection zones obtained from the 2C-scheme with reference to established results.

#### Energy conservation

The equations of radiation hydrodynamics (Eqs. 26 – 31) are discretized conservatively. Hence, the scheme conserves mass, momentum, internal energy, as well as the moments of radiation. Although analytically equivalent, this does not translate into conservation of the total energy in the discrete case. Usually, this is not crucial for 1D computations. Moreover, the adaptive grid provides a fine grid resolution for all gradients and accordingly minimizes spatial discretization errors. In case of the horizontal components in the 2C-scheme, we now have to consider a very coarse spatial representation where, in particular, advection from one column to another requires some attention.

In the present discretization, advection of total energy consists of three components: internal, radiative, and kinetic energy. The former two are treated accurately by the advection terms in the corresponding equations of internal energy and radiation energy. In contrast, the advection of kinetic energy is modeled only indirectly by density and momentum transport. Analytically, advection of kinetic energy, momentum, and density are related by

 (19)

Integration over a cell volume provides the discrete (approximate) equivalent

 ∑i˜12ρu2iFluxi≃u∑i˜ρuiFluxi−12u2∑i˜ρiFluxi (20)

where the quantities with a tilde are advected over cell boundaries with the transported volumes . Using this formula, we can now compute the kinetic energy effectively transported by the advection of mass and momentum.

In application to horizontal advection from one column to the other, we obtain for column 1

 ˜12ρu2V_fluxθ≃u1˜ρuV_fluxθ−12u21˜ρV_fluxθ (21)

and column 2

 ˜12ρu2V_fluxθ≃u2˜ρuV_fluxθ−12u22˜ρV_fluxθ. (22)

Comparing these two lines, it becomes apparent that the advected kinetic energy is not identical in both columns: a certain flux of density and momentum over the interface between the two columns induces a change in kinetic energy in column 1 as given by Eq. 21, while column 2 experiences a change according to Eq. 22. The advection process therefore creates an error in the kinetic energy balance, and consequently also in the conservation of total energy.

To allow for this deficit in the total energy balance, we compute the difference and place it as a source term into the equation of internal energy

In the simplest case, horizontal transport uses donor cell advection or, more generally, a formalism as in Eq. 18. Assuming that we compute and analogously, i.e. with the same , we can further simplify

where is given by

 ρ∗={(1−λ)¯¯¯ρ1−λ¯¯¯ρ2forV_fluxθ>0(1−λ)¯¯¯ρ2−λ¯¯¯ρ1forV_fluxθ<0 (25)

and is the radially averaged density (see Sect. 4.5). For (donor cell), becomes the upstream value, in which case we could write .

From Eq. 24, we have , i.e.  always acts as a source term for the internal energy. therefore effectively describes the dissipation of kinetic energy in the course of advection from one column to another. This dissipation increases with the radial velocity difference between the two columns. In a convection zone, the two columns hold opposing up- and downdraft motions and is quite large; . In these cases, the dissipation term becomes indispensable to the total energy balance; in the examples presented in Sect. 6, it can account for more than 30% of the energy throughput (i.e. luminosity).

Depending on the direction of the horizontal flow, the dissipated energy is deposited in the receiving column. In doing so, the contribution from Eq. 24 must be divided radially to be consistent with the scalar discretization of the equation of internal energy.

### 2.5 Radial distribution of grid points

In the preceding paragraphs, we constructed the two columns discretization scheme with reference to a given radial distribution of the grid points . We now have to adopt a method to determine these grid point positions.

For obvious reasons, the convective fluid flow prohibits a Lagrangian grid customarily used in stellar models. A spatially fixed Eulerian grid is also poorly suited to our needs for two reasons. First, advection alters the stellar structure. Starting from an initial, purely radiative model, the star shrinks significantly with the onset of advective transport. Secondly, this scheme is intended to be used in computing stellar pulsations, i.e. to follow the convective circulating motion while the entire envelope moves in- and outward in the course of stellar pulsation.

To meet these requirements, the code uses an adaptive grid equation (Dorfi & Drury, 1987) that redistributes the grid points continuously according to the evolving physical structures and therefore provides high resolution as needed, e.g. at photospheric gradients, while following radial movements of these features due to structural changes or stellar pulsation. This adaptive grid equation is solved implicitly together with the physical equations. Since the grid equation is an elliptic differential equation, this approach is only possible with an implicit solving method.

In the application to the two-column scheme, the same grid point distribution is used for both columns. Otherwise horizontal advection would become far more complicated and – a serious issue for the implicit solving algorithm – non-local with respect to the grid index .

Variables from both columns are used as ‘grid-weights’, in particular the grid adapts according to gradients in density, internal energy, and in each column. The grid resolution is therefore increased in both columns in identical ways, even though, in general, only the physical structure in one of them actually demands this high grid resolution.

## 3 Physical equations

The physics within the two-column geometrical setup is computed from the equations of radiation hydrodynamics (e.g. Mihalas & Mihalas, 1984). The radiation field is thereby described using the first three gray moments of the intensity , , and , which correspond to the radiative energy density, radiative flux, and radiative pressure, respectively. An Eddington factor closes the moment equations. Neglecting scattering, the source function of radiation, , is given by the Stefan-Boltzmann law , and and are the Rosseland and Planck mean opacities. The gas pressure and gas temperature are given by the equation of state. Self gravity is described by the gravitational potential , which is assumed to be spherically symmetric, i.e. we do not allow for the (negligible) gravitational interaction between up- and downstreams. is Newton’s gravitational constant. Artificial viscosity, discussed in detail in Sect. 4.1, enters in the form of the viscous pressure tensor .

The system of analytical equations is given by:

Equation of continuity

 ∂∂tρ+\@vec∇⋅(\@vecuρ)=0 (26)

Equation of motion

 ∂∂t(ρ\@vecu)+\@vec∇⋅(\@vecuρ\@vecu)+\@vec∇P+ρ\@vec∇ϕ−4πcκRρ\@vecH+\@vec∇⋅\@tensQ=0 (27)

Equation of internal energy

 ∂∂t(ρe)+\@vec∇⋅(\@vecuρe)+P\@vec∇⋅\@vecu−4πκPρ(J−S)+\@tensQ:\@vec∇\@vecu=0 (28)

Poisson equation

 Δϕ=4πGρ (29)

 ∂∂tJ+\@vec∇⋅(\@vecuJ)+c\@vec∇⋅\@vecH+\@tensK:\@vec∇\@vecu+cκPρ(J−S)=0 (30)

 ∂∂t\@vecH+\@vec∇⋅(\@vecu\@vecH)+c\@vec∇⋅\@tensK+\@vecH⋅\@vec∇\@vecu+cκRρ\@vecH=0 (31)

#### Radiation equations for high optical depths

The difference gradually vanishes with increasing optical depth, i.e. towards the interior of a star.

Therefore, the coupling term between radiative energy and gas energy becomes numerically unresolvable for high optical depths. To derive still the correct contribution from this coupling for the equation of internal energy, the corresponding term is expressed by the radiation energy equation and inserted into the equation of internal energy (Feuchtinger, 1999a). This corresponds to evaluating the sum ‘Equation of energy’ + ‘Radiation energy equation’, where terms with cancel out each other. The conversion factor relates the zeroth moment of the intensity to the radiation energy density

 ∂∂t(ρe+4πcJ)+\@vec∇⋅[\@vecu(ρe+4πcJ)]+ +P\@vec∇⋅\@vecu+4πc\@tensK:\@vec∇\@vecu+4π\@vec∇⋅\@vecH+\@tensQ:\@vec∇\@vecu=0. (32)

Inwards of a predefined stellar depth, the equation of internal energy is substituted with this sum, i.e. Eq. 3 is solved instead of Eq. 28.

## 4 Discrete set of equations

After introducing the analytical form of the equations of radiation hydrodynamics, we now develop their discrete version. Table 1 summarizes the primary variables, the corresponding discrete equations, and the closures of the system. As illustrated in Fig. 1, and , as well as the ‘horizontal’ variables and , are integral quantities for both columns. All other variables, , , , , and exist in duplicates, assigned individually to the two columns.

### 4.1 Artificial viscosity

In the continuum description of fluids, shock fronts – and, in the present case, horizontal shear flows – may become indefinitely sharp. In hydrodynamics codes, the smallest physical length scale is given by the mesh size of the numerical grid; on this length scale, numerical dissipation intrinsic to the spatial discretization becomes effective. In the present implementation, the adaptive grid continuously refines to resolve all gradients properly on the grid. This reduces the intrinsic numerical dissipation and can thus lead to a runaway effect of successively steepening gradients and subsequent grid refinement. It is therefore necessary to include an artificial viscosity as a measure of broadening narrow physical features on a predefined length scale. Consequently, this also limits the maximum grid resolution to which the adaptive grid will be refined to.

In this way, artificial viscosity, by specifying the minimum length scale in the computation, plays a more important role than in usual Lagrangian or Eulerian hydrodynamics codes.

Due to the small overall dissipation of the numerical scheme, it is also sometimes necessary to include some extra viscosity to limit amplitudes and velocities, e.g. of stellar pulsations. However, in the results presented in Sect. 6, the influence of the artificial viscosity always remains negligible and is apparent only in a minor smoothing of velocity spikes.

For the artificial viscosity, the geometry-independent description provided by Tscharnuter & Winkler (1979) was adopted. In this description, modeled by analogy with the ordinary (molecular) fluid viscosity, the viscous pressure tensor reads

 \@tensQ=−μQ([\@vec∇\@vecu]sym−113\@vec∇⋅\@vecu) (33)

where the viscosity coefficient contains parameters for ‘linear’ (pseudo-molecular) viscosity and ‘quadratic’ viscosity (where ‘quadratic’ refers to the quadratic dependency on the velocity field, which causes it to act in a way similar to a turbulent viscosity)

The use of the maximum implies that expanding flows are unaffected by viscosity. The viscous length scale is set to the characteristic extension of the problem (and of the numerical grid), e.g. the radius in spherical geometry. is the local speed of sound. For the symmetric velocity gradient, the notation was introduced. The symmetric description ensures that rotation, which does not affect the physical structure, remains unaffected by viscosity

 [\@vec∇\@vecu]sym=12(\@vec∇\@vecu+(\@vec∇\@vecu)T). (35)

The contributions of artificial viscosity to the equations of motion and internal energy follow directly from the viscous pressure tensor. The viscous force is computed to be the divergence of the viscous pressure

 \@vecfQ=\@vec∇⋅\@tensQ. (36)

The viscous energy dissipation is obtained by contraction of the viscous pressure tensor with the gradient of the velocity field. Since is symmetric, there is no difference between using the velocity gradient or the symmetric velocity gradient

 ϵQ=\@tensQ:\@vec∇\@vecu. (37)

To apply this recipe in the present case, Eqs. 33 – 37 must be evaluated by assuming spherical geometry. Since the 2C-scheme describes all types of horizontal variability and dynamics with only one interface between the two radial columns, the - and -directions of spherical coordinates are not considered separately and we adopt the identities and . To allow for that, all derivatives in the -direction are assumed to be taken on the great circle, i.e. for . Also note that is symmetric by definition and we therefore finally have four independent entries for the viscous pressure tensor in the two-column geometry: , , , and .

In principle, viscosity couples fluid flows in different coordinate directions. In the 2C-scheme, due to the combined discretization of - and -components, the corresponding terms in the spherical symmetric description become ambiguous in interpretation; the two-column representation of the three-dimensional flow is too simplistic to enable a proper modelling of this effect. The viscous interaction between the two directions of fluid flow was therefore neglected by assuming for the viscosity in the -direction, and in the radial direction. For the viscosity in the radial direction, we then obtain

 Qrr = −μQ23(∂ur∂r−urr) (38) Qrθ = −μQ121r∂ur∂θ (39)
 fQr=3r∂∂r3(r3Qrr)+2r∂Qrθ∂θ (40)
 ϵQr=−μQ23(∂ur∂r−urr)2−μQ(1r∂ur∂θ)2. (41)

For the viscosity in the -direction, we arrive at

 Qrθ = −μQ12(∂uθ∂r−uθr) (42) Qθθ = −μQ131r∂uθ∂θ (43)
 fQθ=23r∂∂r3(r3Qrθ)+21r∂∂θ(4Qθθ) (44)
 ϵQθ=−μQ(∂uθ∂r−uθr)2−μQ423(1r∂uθ∂θ)2. (45)

In Eq.44, an additional factor 2 was included for because we are considering forces in both the - and -directions, even though they are combined in the discretization process.

In the discrete case, derivatives with respect to radius transform into differences between radial indices (), and derivatives in the -direction are discretized using Eqs. 3 & 4.

The various terms of the artificial viscosity (radial – horizontal, shear – non-shear) include separate coefficients to allow for their individual adjustment. Table 2 presents the coefficients with the parameters on which they depend. The computation of the coefficients is similar to that described by Eq. 34, except for details related to the staggered-mesh location of the involved variables; the turbulent (‘quadratic’) viscosity parameter is only used for the radial, non-shear part ( and ). Where appropriate, the ’s are evaluated separately in each column, although the viscosity parameters are the same in both columns. In total, there are 5 viscosity parameters, although until now only three (except for testing purposes) were actually used in the computations. The default values adopted in the examples presented in Sect. 6 are , , and .

The discretization of the viscous force and energy dissipation uses the same discretization volumes as the corresponding equations, i.e. the equation of motion and equation of internal energy. To emphasize the volume-integrated variables, the discrete forces and energies are denoted by capital letters; and .

For the viscous force in the radial direction, we obtain for column 1

 FQ1 = −cf18π3rΔr{μQ1¯¯¯r3(Δru1Δrr−¯¯¯u1¯¯¯r)} (46) +μQshearN4r2(u1−u2)12V_vol

and for column 2

 FQ2 = −cf28π3rΔr{μQ2¯¯¯r3(Δru2Δrr−¯¯¯u1¯¯¯r)} (47) −μQshearN4r2(u1−u2)12V_vol.

Note that both shear forces are discretized with instead of and to allow them to cancel out each other for the two columns.

The corresponding viscous energy dissipation reads

 EQ1 = −μQ123(Δru1Δrr−¯¯¯u1¯¯¯r)2S_vol1 (48) −μQshearN8(¯¯¯u1−¯¯¯u2¯¯¯r)2S_vol1 EQ2 = −μQ223(Δru2Δrr−¯¯¯u1¯¯¯r)2S_vol2 (49) −μQshearN8(¯¯¯u1−¯¯¯u2¯¯¯r)2S_vol2.

This formalism for the viscosity in the radial direction closely resembles – except of course for the shear part – the customary 1D viscosity description given, e.g., by Dorfi (1998) or Feuchtinger (1999a).

In the horizontal direction, discretization yields a viscous force

 FQθ = −2π¯¯¯rΔr{μQθshearr3(ΔruθΔr¯¯¯r−¯¯¯uθr)} (50) +μQθ4N3uθ¯¯¯r2H_vol,

and a viscous energy dissipation

 EQθ1 = −μQθshear(Δr¯¯¯uθΔrr−uθ¯¯¯r)2S_vol1 (51) −μQθ4N3u2θ¯¯¯r2S_vol1 EQθ2 = −μQθshear(Δr¯¯¯uθΔrr−uθ¯¯¯r)2S_vol2 (52) −μQθ4N3u2θ¯¯¯r2S_vol2.

In the moment description of radiation, the second moment of the intensity – which corresponds to the radiation pressure – is assumed to be of the following form

 \@tensK=⎛⎜⎝KrrKθθKϕϕ⎞⎟⎠=⎛⎜ ⎜ ⎜⎝KrrJ−Krr2J−Krr2⎞⎟ ⎟ ⎟⎠ (53)

with the radial component given by a scalar Eddington factor

 Krr=feddJ. (54)

The same Eddington factor is taken for both columns

 Krr,1=feddJ1Krr,2=feddJ2, (55)

and, for the examples presented in this paper, it has been set to a constant value of for simplicity.

The discrete equation of radiative flux in the horizontal direction (i.e. for ), does not use the full time-dependent equation Eq. 31 but only its stationary part

 \@vec∇⋅\@tensK+κRρ\@vecH=0. (56)

In this way, we did not have to discretize the horizontal component of , which is, in a similar way to artificial viscosity, ambiguous in interpretation in the context of the 2C-scheme. Considering the simplistic discretization of horizontal exchange between the two columns, this stationary, diffusion-like description remains sufficient.

Note that the assumption for the radiative pressure in Eq. 53 will in general not be consistent with the horizontal radiative flux computed from Eq. 56. However, a more consistent description is not reasonably possible given the limited resolution in horizontal direction of the 2C-scheme. Adopting a more elaborate description would also require solving the detailed 2D radiative transport to obtain the required Eddington factors (e.g.  and ). And after all, there is no point in improving the radiative transport beyond the level of approximation of the hydrodynamics part.

### 4.3 The discrete equations of radiation hydrodynamics

Using the discretization scheme presented in Sect. 2 and the results from Sect. 4.1 and Sect. 4.2, we obtain the discrete version of Eqs. 26 – 31 & Eq. 56. Table 3 provides the full discrete set of equations of radiation hydrodynamics. These physical equations are completed by the equations for the radially averaged densities, Eqs. 60 & 61, and by the adaptive grid equation.

As an example of the discrete form of conservative equations and to illustrate the notation of the advective contributions, the discrete equations of continuity are given here for both columns

 δ[ρ1S_vol1]+[˜ρ1S_flux1]i+1−[˜ρ1S_flux1]i+˜ρ1S_fluxθ=0 (57)
 δ[ρ2S_vol2]+[˜ρ2S_flux2]i+1−[˜ρ2S_flux2]i−˜ρ2S_fluxθ=0. (58)

The notation indicates a difference of between the new and old time level separated by the time step . Advection represented by the terms with S_flux, occurs in the radial direction both at the radii and as well as in the horizontal direction over the interface between the two columns. Note that the horizontal transport terms with correspond to each other.

In Table 3, an abbreviated notation was adopted for the advective terms by summarizing all three contributions. In this form, advective terms, e.g. those from Eq. 57, are written as

 ∑˜ρ1S_flux1,θ≡[˜ρ1S_flux1]i+1−[˜ρ1S_flux1]i+˜ρ1S_fluxθ. (59)

Spatial differences are denoted as and in the radial and horizontal direction, respectively. Averaged quantities – where the precise definition depends on the context – are written with overhead dashes.