Moving and reactive boundary conditions in moving-mesh hydrodynamics

ABSTRACT

We outline the methodology of implementing moving boundary conditions into the moving-mesh code manga. The motion of our boundaries is reactive to hydrodynamic and gravitational forces. We discuss the hydrodynamics of a moving boundary as well as the modifications to our hydrodynamic and gravity solvers. Appropriate initial conditions to accurately produce a boundary of arbitrary shape are also discussed. Our code is applied to several test cases, including a Sod shock tube, a Sedov–Taylor blast wave, and a supersonic wind on a sphere. We show the convergence of conserved quantities in our simulations. We demonstrate the use of moving boundaries in astrophysical settings by simulating a common envelope phase in a binary system, in which the companion object is modelled by a spherical boundary. We conclude that our methodology is suitable to simulate astrophysical systems using moving and reactive boundary conditions.

hydrodynamics, methods: numerical

1 INTRODUCTION

Numerical simulations have played a crucial role in understanding the hydrodynamics of many astrophysical phenomena, including gas, stars, discs, galaxies, and large-scale structure. Two methodologies are commonly employed to solve the hydrodynamic equations: smooth particle hydrodynamics (SPH) and grid-based solvers. SPH codes have excellent conservation properties but resolve shocks poorly due to their smoothing nature. On the other hand, grid-based solvers are superior at modelling shocks due to the use of Godunov schemes. However, they suffer from a lack of Galilean invariance and errors in angular momentum conservation. Furthermore, even with adaptive mesh refinement, grid-based solvers encounter difficulty with highly refined regions that move with large velocities relative to the mesh. This is due to the challenges associated with adapting the mesh to the anticipated motion of the system (Springel 2010).

In recent years, a hybrid of these methods has been developed in an attempt to capture the best characteristics of both. This is the arbitrary Lagrangian–Eulerian (ALE) scheme, and software that use ALE schemes are known as moving-mesh codes (for a review, see Donea et al. 2004). In an ALE scheme, the mesh moves along with the fluid, combining the superior shock capturing of grid-based solvers with the conservation properties of SPH. Springel (2010, hereinafter S10) described one such scheme that has proven successful, which is implemented in the code arepo. The scheme constructs an unstructured mesh from an arbitrary distribution of points using a Voronoi tessellation (Okabe, Boots & Sugihara 1992). This guarantees that the mesh will be well defined, unique, and continuously deformable; thus, finite-volume methods can be applied in a manner similar to that of an Eulerian code. In addition, the lack of Galilean invariance in Eulerian codes is rectified if the mesh cells move along with the local flow. It has also been argued that ALE schemes are superior at capturing boundary layer instabilities such as Kelvin–Helmholtz instabilities (S10), though Lecoanet et al. (2016) cautioned that numerical noise can masquerade as solutions of such instabilities. Although arepo was originally developed with cosmological simulations in mind (see for instance Vogelsberger et al. 2014), it has also been used in a number of problems including stellar mergers (Zhu et al. 2015), tidal disruption events (Goicovic et al. 2019), and common envelope evolution (CEE) (Ohlmann et al. 2016).

A moving-mesh hydrodynamic solver has also been developed for the N-body simulation code changa (Charm N-body GrAvity solver) (Jetley et al. 2008, 2010; Menon et al. 2015). This moving-mesh solver is known as manga, which is described in detail by Chang, Wadsley & Quinn (2017, hereinafter CWQ17). In Prust & Chang (2019, hereinafter PC19), we incorporate individual time-steps into manga as well as the equations of state implemented in the open-source stellar evolution code mesa (Paxton et al. 2011, 2013, 2015, 2018, 2019). We use these tools to simulate CEE in binary star systems.

In fluid dynamics, it is often desirable to set boundary conditions on a surface. For example, in modelling the flow around an airfoil, the edges of the simulation box can be treated as inflow or outflow boundaries and the airfoil itself as a reflective boundary, which allows no matter to pass through it. Indeed, Wan & Zhong (2019) use a 2D ALE scheme with moving boundaries to model a pitching NACA 0012 airfoil. Immersed boundary methods can also be useful, especially when the relative displacement of the boundaries is large (Hedayat & Borazjani 2019).

Such schemes are common in fields such as aerodynamics, but they also have interesting and useful applications in astrophysics. For example, Chen, Ivanova & Carroll-Nellenback (2020) use two separate boundary conditions in their simulations of asymptotic-giant-branch (AGB) binaries using the adaptive mesh refinement code astrobear. On a cylindrical surface surrounding the binary, they set the fluid velocity to be supersonic in the radially outward direction, ensuring that no information propagates into the cylinder from the outside. They also set an inner boundary condition beneath the photosphere of the AGB star requiring that the radial fluid velocity oscillate in time, which approximates stellar pulsations (Bowen 1988). Additionally, outflow conditions can be useful in astrophysical problems where accretion is present, such as in the growth of intermediate-mass black holes in the early universe (for example, Toyouchi et al. 2019). However, due to the difficulty of implementing moving boundaries, such phenomena are typically modelled using static boundary conditions even in cases where moving boundaries are a more natural fit to the system in question.

The implementation of boundary conditions into Eulerian grid-based solvers, which we discuss in Section 3, is straightforward. However, because Eulerian grid codes have a static mesh, it is at best very difficult to have a boundary that moves with respect to the mesh. On the other hand, moving-mesh codes are ideal for moving boundaries since the mesh-generating points can move with respect to one another. In ALE schemes, the velocities of the mesh-generating points are typically taken to be equal or similar to the fluid velocity to gain the advantages offered by a Lagrangian scheme. However, the mesh-generating points can in principle move at any velocity relative to the fluid.

S10 carried out a simple 2D simulation using moving boundaries in arepo in which a curved boundary (thought of as a ‘spoon’) stirs a fluid. They generate their boundary using two sets of points on either side of and equidistant from the desired boundary surface (see their fig. 38), on which they set reflective boundary conditions. They also require that both sets of points move together as a rigid body and move their spoon at a constant speed along a pre-determined circular path.

We adopt a similar strategy for our implementation of moving boundaries into manga, requiring the motion of the cells adjacent to the boundary to be identical to the motion of the boundary itself. However, we also set the velocity of the boundary such that its motion is reactive to external forces. The aerodynamic drag force can be computed by integrating the gas pressure over the surface of the boundary, and the gravitational force by integrating over the volume of the boundary.

We have organized this paper as follows. In Section 2, we summarize the algorithm of manga prior to the implementation of moving boundary conditions. We discuss the changes to the hydrodynamic and gravity solvers to implement moving boundaries in Section 3. In Section 4, we run simulations of four test cases: a supersonic wind, a Sod shock tube, a Sedov–Taylor explosion, and a common envelope phase in a binary star system. We discuss our results and possible future directions in Section 5.

2 PREVIOUS ALGORITHM

The ALE algorithm that is implemented in manga is briefly summarized as follows. We refer the reader to CWQ17 and PC19 for detailed discussions. manga solves the Euler equations, which written in conservative form are

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla }\cdot \rho \boldsymbol{v}= 0, \end{eqnarray}$$

(1)

$$\begin{eqnarray} \frac{\partial \rho \boldsymbol{v}}{\partial t} + \boldsymbol{\nabla }\cdot \rho \boldsymbol{v}\boldsymbol{v}+ \boldsymbol{\nabla }P = -\rho \boldsymbol{\nabla }\Phi , \end{eqnarray}$$

(2)

and

$$\begin{eqnarray} \frac{\partial \rho e}{\partial t} + \boldsymbol{\nabla }\cdot \left(\rho e + P\right)\boldsymbol{v}&=& -\rho \boldsymbol{v}\cdot \boldsymbol{\nabla }\Phi , \end{eqnarray}$$

(3)

where ρ is the density, |$\boldsymbol{v}$| is the fluid velocity, Φ is the gravitational potential, e = ϵ + v²/2 is the specific energy, ϵ is the specific internal energy, and P(ρ, ϵ) is the pressure. Equations (1) to (3) can be written in a compact form by introducing a state vector |$\boldsymbol{\mathcal {U}}=(\rho , \rho \boldsymbol{v}, \rho e)$|⁠:

$$\begin{eqnarray} \frac{\partial \boldsymbol{\mathcal {U}}}{\partial t} + \boldsymbol{\nabla }\cdot \boldsymbol{\mathcal {F}}= \boldsymbol{\mathcal {S}}, \end{eqnarray}$$

(4)

where |$\boldsymbol{\mathcal {F}}=(\rho \boldsymbol{v}, \rho \boldsymbol{v}\boldsymbol{v}+ P, (\rho e + P)\boldsymbol{v})$| is the flux function and |$\boldsymbol{\mathcal {S}}= (0, -\rho \boldsymbol{\nabla }\Phi , -\rho \boldsymbol{v}\cdot \boldsymbol{\nabla }\Phi$|⁠) is the source function. To solve equation (4), we adopt the same finite-volume strategy as S10. For each cell, the integral over the volume of the ith cell V_i defines the charge of the ith cell to be

$$\begin{eqnarray} \boldsymbol{U}_i = \int _{V_{i}} \boldsymbol{\mathcal {U}}dV = \boldsymbol{\mathcal {U}}_i V_i. \end{eqnarray}$$

(5)

We then use Gauss’ theorem to convert the volume integral over the divergence of the flux in equation (4) to a surface integral

$$\begin{eqnarray} \int _{V_{i}} \boldsymbol{\nabla }\cdot \boldsymbol{\mathcal {F}}dV = \int _{\partial V_{i}} \boldsymbol{\mathcal {F}}\cdot \hat{\boldsymbol{n}}dA, \end{eqnarray}$$

(6)

where ∂V_i is the boundary of the cell. We now take advantage of the fact that the volumes are Voronoi cells with a finite number of neighbours to define an integrated flux

$$\begin{eqnarray} \sum _{j \in {\rm neighbors}} \boldsymbol{F}_{ij} A_{ij} = \int _{\partial V_{i}} \boldsymbol{\mathcal {F}}\cdot \hat{\boldsymbol{n}}dA, \end{eqnarray}$$

(7)

where |$\boldsymbol{F}_{ij}$| and A_ij are the average flux and area of the common face between cells i and j. The discrete time evolution of the charges in the system is given by

$$\begin{eqnarray} \boldsymbol{U}_i^{n+1} = \boldsymbol{U}_i^n + \delta t \sum _j \hat{\boldsymbol{F}}_{ij} A_{ij} + \Delta t\boldsymbol{S}_i, \end{eqnarray}$$

(8)

where |$\hat{\boldsymbol{F}}_{ij}$| is an estimate of the half time-step flux between the initial |$\boldsymbol{U}_i^n$| and final states |$\boldsymbol{U}^{n+1}_i$| and |$\boldsymbol{S}_i = \int _{V_{i}} \boldsymbol{\mathcal {S}}dV$| is the time-averaged integrated source function.

We estimate the flux |$\hat{\boldsymbol{F}}_{ij}$| across each face as follows.

• Use the gradient estimates at the initial time-step to predict the half time-step cell-centred values.

• Drift the cells a half time-step and rebuild the Voronoi tessellation at the half time-step.

• Estimate the half time-step state vector (in the rest frame of the moving face) at the face centre |$\boldsymbol{\tilde{r}}_{ij}$| between the neighbouring i and j cells by linear reconstruction.

• Estimate the (half time-step) velocity |$\boldsymbol{\tilde{w}}_{ij}$| of the face, following S10, and boost the state vector from the ‘lab’ frame (the rest frame of the simulation box) to the rest frame of the face to find the flux along the normal of the face, i.e. in the direction from i to j.

• Estimate the flux |$\hat{\boldsymbol{F}}_{ij}$| across the face using an Harten-Lax-van Leer-Contact (HLLC) or Harten-Lax-van Leer (HLL) approximate Riemann solver implemented following Toro (2009).

• Boost the solved flux back into the ‘lab’ frame.

We can then use the estimated fluxes to time-evolve the charges |$\boldsymbol{U}_i$| following equation (8) using the full time-step δt and apply changes owing to the source terms.

To conclude our summary of manga, we outline the main computation loop as a whole.

• We determine a valid Voronoi tessellation of the mesh-generating points using the voro++ library.

• Using the volume of the Voronoi cell and integral quantities, |$\boldsymbol{U}_i^n$|⁠, the conserved and primitive variables are determined. The local gradients and half time-step conserved variables are calculated.

• The flux across the faces is calculated as discussed above.

• The new state |$\boldsymbol{U}^{n+1}_i$| is determined from the fluxes. New velocities for the mesh-generating points are determined. The mesh-generating points drift by a full time-step.

3 HYDRODYNAMICS OF A MOVING BOUNDARY

We now discuss our changes to the manga algorithm to implement moving boundary conditions. The treatment is restricted to reflective boundaries, though the extension to inflow or outflow conditions is straightforward.

3.1 Flux across a reflective boundary

Let us consider a face on which we would like to set a boundary condition and boost our state vectors and flux functions into the rest frame of the face. If we rotate the face such that its normal lies along the x-direction, as is done in manga, then we can write the flux normal to the face as

$$\begin{eqnarray} \boldsymbol{\mathcal {F}}_{x} = (\rho v_{x},\,\rho v_{x}^{2} + P,\,\rho v_{x} v_{y},\,\rho v_{x} v_{z},\,(\rho e + P) v_{x}). \end{eqnarray}$$

(9)

Here, we have expanded the momentum flux tensor |$\rho \boldsymbol{v}\boldsymbol{v}$| into its three components. The flux of mass, momentum, and energy that pass through the face is determined by the values of the state vector |$\boldsymbol{\mathcal {U}}_{\rm {L}}$|⁠, |$\boldsymbol{\mathcal {U}}_{\rm {R}}$|⁠, and the flux function |$\boldsymbol{\mathcal {F}}_{x,\rm {L}}$|⁠, |$\boldsymbol{\mathcal {F}}_{x,\rm {R}}$| on either side of the face. We see that under the transformation v_x → −v_x, all terms in |$\boldsymbol{\mathcal {F}}_{x}$| change sign with the exception of the x-momentum flux |$\rho v_{x}^{2} + P$|⁠. Conversely, when this transformation is applied to the state vector |$\boldsymbol{\mathcal {U}}=(\rho , \rho v_{x}, \rho v_{y}, \rho v_{z}, \rho e)$|⁠, the only term that changes sign is the x-momentum ρv_x. Consider now the HLL approximate Riemann solver (Einfeldt 1988), which computes the flux across the face as

$$\begin{eqnarray} \hat{\boldsymbol{F}}_{\rm {HLL}} = \frac{ \alpha _{+}\boldsymbol{\mathcal {F}}_{x,\rm {L}} + \alpha _{-}\boldsymbol{\mathcal {F}}_{x,\rm {R}} - \alpha _{+}\alpha _{-}(\boldsymbol{\mathcal {U}}_{\rm {R}}-\boldsymbol{\mathcal {U}}_{\mathrm{ L}}) }{ \alpha _{+}+\alpha _{-} }, \end{eqnarray}$$

(10)

where |$\alpha _{\pm } = \rm {max}[0, \pm \lambda _{\pm }(\boldsymbol{\mathcal {U}}_{\rm {L}}), \pm \lambda _{\pm }(\boldsymbol{\mathcal {U}}_{\rm {R}})]$|⁠, the maximum and minimum eigenvalues of the Jacobians of each state are λ_± = v_x ± c_s and c_s is the sound speed. We now impose the condition that the state vector and flux on either side are equal except for the fluid velocities normal to the face, which are opposite: |$v_{x} = v_{x,\rm {L}} = -v_{x,\rm {R}}$|⁠. That is,

$$\begin{eqnarray} \boldsymbol{\mathcal {F}}_{x,\rm {R}}(v_{x}) = \boldsymbol{\mathcal {F}}_{x,\rm {L}}(-v_{x}) \quad \text{and}\quad \boldsymbol{\mathcal {U}}_{\rm {R}}(v_{x}) = \boldsymbol{\mathcal {U}}_{\rm {L}}(-v_{x}). \end{eqnarray}$$

(11)

This implies α₊ = α₋ = c_s + |v_x| and simplifies the HLL flux to the form

$$\begin{eqnarray} \hat{\boldsymbol{F}}_{\rm {HLL}} = \frac{ \boldsymbol{\mathcal {F}}_{x,\rm {L}}(v_{x}) + \boldsymbol{\mathcal {F}}_{x,\rm {L}}(-v_{x}) - \alpha _{+}[\boldsymbol{\mathcal {U}}_{\rm {L}}(-v_{x})-\boldsymbol{\mathcal {U}}_{\rm {L}}(v_{x})] }{2}. \end{eqnarray}$$

(12)

If we now plug in |$\boldsymbol{\mathcal {F}}$| and |$\boldsymbol{\mathcal {U}}$| in component form as in (9), we obtain

$$\begin{eqnarray} \hat{\boldsymbol{F}}_{\rm {HLL}} = \left(0,\,\rho v_{x}^{2}+P+\rho v_{x}(c_{\rm {s}}+|v_{x}|),\,0,\,0,\,0\right). \end{eqnarray}$$

(13)

We see that all terms in |$\hat{\boldsymbol{F}}_{\rm {HLL}}$| vanish except for the term corresponding to the flux of x-momentum – that is, the component of momentum normal to the face. Thus, the reflective boundary condition (equation 11) prevents the transport of any quantity besides momentum normal to the face for the HLL Riemann solver; mass, energy, and entropy cannot be transported. This momentum flux corresponds to the contact pressure exerted on the face. The derivation for the HLLC solver (Toro, Spruce & Speares 1994) proceeds in a similar fashion. The implementation of reflective boundary conditions requires modifications to the Riemann solver and – for schemes that are accurate to second order or higher in space – the gradient estimation, both of which are discussed below. The alternate condition |$\boldsymbol{\mathcal {F}}_{x,\rm {L}}(v_{x}) = \boldsymbol{\mathcal {F}}_{x,\rm {R}}(v_{x})$| and |$\boldsymbol{\mathcal {U}}_{x,\rm {L}}(v_{x}) = \boldsymbol{\mathcal {U}}_{x,\rm {R}}(v_{x})$| corresponds to an outflow boundary condition. However, we focus on reflective conditions in this paper for the reason that inflow or outflow conditions in a moving-mesh code would necessitate the creation or destruction of mesh-generating points, which is not yet implemented in manga. A separate publication (Chang, Davis & Jiang 2020) details the development of inflow and outflow conditions for radiation, which is applicable to fluids on a static mesh.

3.2 Initialization

In the initialization step at the beginning of the simulation, we declare some user-specified region as the initial boundary. Cells with centres within this region are defined to be a part of the boundary. We will refer to these cells as ‘boundary cells’ and to all other cells as ‘gas cells’. The boundary cells will comprise the boundary for the entirety of the simulation, as boundary cells are never added or removed. We define the boundary particles to have equal mass so that the density of the boundary is constant, provided that the density of mesh-generating points is constant within the boundary. The total mass of the boundary is an input parameter in our simulations.

To accurately produce the desired boundary shape, it is necessary to have a sufficiently high density of mesh-generating points near the edges of the boundary. To this end, it is sometimes beneficial to initialize the simulation with a larger resolution in the region containing the boundary than in the rest of the simulation box. We demonstrate this in Section 4.4. In principle, a boundary of arbitrary shape can be constructed, given a sufficiently high density of mesh-generating points.

3.3 Edge cells

We define all gas cells that share a face with a boundary cell as ‘edge cells’. These cells function as normal gas cells for the most part; however, the motion of their mesh-generating points always matches that of the boundary. The result of this is a layer of gas cells that do not move with respect to the boundary, as in S10. This is important because the surface of the boundary comprises the shared faces between the edge cells and the boundary cells, which we refer to as ‘boundary faces’. The boundary faces are redrawn with every Voronoi computation, which happens several times during each time-step. Therefore, if the mesh-generating points were allowed to flow past the boundary, the surface of the boundary would deform, which is undesirable. Although the edge cells are the only gas cells that have direct contact with the boundary, they communicate its effects to the rest of the gas through their properties and motion.

3.4 Gradient estimation and linear reconstruction

A second-order accurate code in space demands an appropriate estimate for the state vector at the face centres of each cell. The state vector on the face between the ith cell and its jth neighbour is taken to be

$$\begin{eqnarray} \tilde{\boldsymbol{\mathcal {U}}}_{ij} = \boldsymbol{\mathcal {U}}_{i} + (\tilde{\boldsymbol{r}}_{ij} - \boldsymbol{c}_{i}) \cdot \nabla \boldsymbol{\mathcal {U}}_{i}, \end{eqnarray}$$

(14)

where |$\tilde{\boldsymbol{r}}$| and |$\boldsymbol{c}_{i}$| are the centroids of the face and ith cell, respectively. Equation (14) requires an estimate of the gradient of the state vector in the ith cell. To calculate the gradients, we follow the procedure of Steinberg, Yalinewich & Sari (2016), who improved upon the prescription of S10 in using the Gauss–Green theorem to estimate these gradients. We refer the reader to CWQ17 for a more detailed discussion.

We modify the gradient estimation algorithm for edge cells by modifying the properties of their neighbouring boundary cells. Consider an edge cell with state vector |$\boldsymbol{\mathcal {U}}_{i}$|⁠, fluid velocity |$\boldsymbol{v}_{i}$|⁠, and neighbours with states |$\boldsymbol{\mathcal {U}}_{j}$|⁠, where the normal vectors of the faces shared with its neighbours are |$\hat{\boldsymbol{n}}_{ij}$|⁠. For each boundary cell neighbour, we first copy the state vector of the edge cell into the boundary cell |$(\boldsymbol{\mathcal {U}}_{j} = \boldsymbol{\mathcal {U}}_{i})$|⁠. We then modify |$\boldsymbol{\mathcal {U}}_{j}$| such that the component of its fluid velocity |$\boldsymbol{v}_{j}$| normal to the shared face is equal and opposite to |$\boldsymbol{v}_{i}$|⁠. That is,

$$\begin{eqnarray} \boldsymbol{v}_{j} = \boldsymbol{v}_{i} - 2 (\hat{\boldsymbol{n}}_{ij} \cdot \boldsymbol{v}_{i}) \hat{\boldsymbol{n}}_{ij}, \end{eqnarray}$$

(15)

which enforces the condition (11) for reflective boundary conditions.

We also modify the slope limiter in the presence of a boundary. manga uses a slope limiter that reduces numerical oscillations near strong gradients

$$\begin{eqnarray} \langle \nabla \boldsymbol{\mathcal {U}}\rangle _{i}^{S10} = \alpha _{i}^{S10} \langle \nabla \boldsymbol{\mathcal {U}}\rangle _{i}, \end{eqnarray}$$

(16)

where |$\alpha _{i}^{S10}$| is defined in S10. As noted by S10, this limiter is not total variation diminishing (TVD), so spurious oscillations can still occur. For this reason, we follow the suggestion by Duffell & MacFadyen (2011) and apply an additional slope limiter

$$\begin{eqnarray} \langle \nabla \boldsymbol{\mathcal {U}}\rangle _{i}^{\mathrm{ DM}} = \alpha _{i}^{\mathrm{ DM}} \langle \nabla \boldsymbol{\mathcal {U}}\rangle _{i}^{S10}, \end{eqnarray}$$

(17)

where

$$\begin{eqnarray} \alpha _{i}^{\mathrm{ DM}} = \min \left(1, \psi ^{\prime }_{ij}\right). \end{eqnarray}$$

(18)

Here,

$$\begin{eqnarray} \psi ^{\prime }_{ij} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\max [\theta (\boldsymbol{\mathcal {U}}_{j}-\boldsymbol{\mathcal {U}}_{i})/\Delta \boldsymbol{\mathcal {U}}_{ij}, 1] & |\Delta \boldsymbol{\mathcal {U}}_{ij}| \gt 0 \\ 1 & \mbox{otherwise} \end{array}\right. \end{eqnarray}$$

(19)

and

$$\begin{eqnarray} \Delta \boldsymbol{\mathcal {U}}_{ij} = \langle \nabla \boldsymbol{\mathcal {U}}\rangle _{i} \cdot (\boldsymbol{c}_{j} - \boldsymbol{c}_{i}). \end{eqnarray}$$

(20)

This limiter is TVD if θ < 0.5, and setting θ = 1 reduces it to the S10 limiter. For all simulations presented in this paper, we use θ = 0.49. However, when performing a reconstruction on a face between an edge cell and a boundary cell we set θ = 0, effectively turning off this additional slope limiter. We find that this marginally improves the conservation properties of our simulations and the stability of the code.

3.5 Time-stepping

In PC19, we discuss the implementation of individual time-steps into manga, which we summarize here. The basic (universal) time-stepping algorithm for a second-order accurate (in time) integrator can be broken down into three stages, the initial time (a), the half time-step (b), and the full time-step (c).

• Initial time t = 0: Determine Voronoi cells using current positions of mesh-generating points. Estimate gradients and construct half time-step predictions. Zero out the changes to the charges, e.g. |$\delta \boldsymbol{U}=0$|⁠.

• Half time-step t = 0.5δt: Drift Voronoi cells to half time-step positions. Construct Voronoi cells at half time-steps. Compute the gradients and perform linear reconstruction to the half time-step cell faces to compute fluxes. Perform a Riemann solution and incorporate source terms with the full time-step. Place these changes into |$\delta \boldsymbol{U}$|⁠.

• Full time-step t = δt: Drift Voronoi cells to full time-step positions. Advance charges to be |$\boldsymbol{U}^{n+1} = \boldsymbol{U}^{n} + \delta \boldsymbol{U}$|⁠. Reset the state to be the new initial time.

We note that the full time-step for any cell involves actions at a full time-step (a,c) and the half time-step (b). For individual time-steps, we assign to each particle a time-step level that we refer to as a ‘rung’, where each rung has a smaller time-step than the rung below it by a factor of 2. Therefore, at the half time-step for some rung i, full time-steps are executed for all rungs j > i. Cells are allowed to change rungs only at their full time-steps. Additionally, we limit the time-step of any cell to be less than |$\sqrt{2}$| of the minimum time-step of its neighbours. This ensures that the time-step changes by no more than a factor of 2 (and the rung by no more than 1) over a distance of two cells. This is similar to the time-step smoothing over an SPH kernel used by Saitoh & Makino (2009), and we find that it allows for stable integrations in manga. PC19 finds that individual time-stepping decreases the computation time by a factor of 4 to 5 for a CEE simulation.

When a moving boundary is present, we place the boundary and edge cells on the highest rung (smallest time-step) of all cells present. Besides improving the stability of the code, this mitigates the risk of gas cells penetrating into the boundary or other pathological behaviours. Because the number of edge cells is typically small relative to the total number of cells, we find that this does not substantially slow down the simulation code.

3.6 Riemann solver

manga uses either an HLLC or HLL approximate Riemann solver to estimate the flux through each face of a cell at the half time-step. We have made several modifications to the Riemann solver in the presence of a boundary.

• A Riemann solution is not performed for boundary cells. Because the state vector |$\boldsymbol{\mathcal {U}}$| of a boundary cell is determined by those of its neighbours, as described in Section 3.4, there is no need to compute the time rate of change of the boundary cell charges.

• For boundary faces, we modify the state vector |$\boldsymbol{\mathcal {U}}$| as well as the flux function |$\boldsymbol{\mathcal {F}}$| inside of the boundary cell according to equation (11). The reconstruction gives the states on either side of the boundary face in the edge cell and boundary cell |$\boldsymbol{\mathcal {U}}_{\rm {e}}$|⁠, |$\boldsymbol{\mathcal {U}}_{\rm {b}}$| as well as the flux functions |$\boldsymbol{\mathcal {F}}_{\rm {e}}$|⁠, |$\boldsymbol{\mathcal {F}}_{\rm {b}}$|⁠. We rotate these such that the normal direction to the face lies along the x-direction. To enforce reflective boundary conditions, we ignore the |$\boldsymbol{\mathcal {U}}_{\rm {b}}$| and |$\boldsymbol{\mathcal {F}}_{\rm {b}}$| given by the reconstruction and instead set them using condition (11). Denoting the (rotated) x-velocity in the edge cell as |$v_{x,\rm {e}}$|⁠, we require

  $$\begin{eqnarray} \boldsymbol{\mathcal {U}}_{\rm {b}} = \boldsymbol{\mathcal {U}}_{\rm {e}}(-v_{x,\rm {e}}) \quad \text{and}\quad \boldsymbol{\mathcal {F}}_{\rm {b}} = \boldsymbol{\mathcal {F}}_{\rm {e}}(-v_{x,\rm {e}}). \end{eqnarray}$$

  (21)

  The new |$\boldsymbol{\mathcal {U}}_{\rm {b}}$| and |$\boldsymbol{\mathcal {F}}_{\rm {b}}$| along with |$\boldsymbol{\mathcal {U}}_{\rm {e}}$| and |$\boldsymbol{\mathcal {F}}_{\rm {e}}$| are then used by the Riemann solver to estimate the flux |$\hat{\boldsymbol{F}}$| across the face, and the result is used in equation (8) to evolve the charges in time.
• The contact pressure P_con at each boundary face, computed by the Riemann solver as in (13), is given by

  $$\begin{eqnarray} P_{\rm {con}}= \bar{P} - \frac{1}{2} \hat{\boldsymbol{n}}\cdot (\boldsymbol{v}_{\rm {b}}-\boldsymbol{v}_{\rm {e}}) \bar{\rho } \bar{c_{\rm {s}}} + \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0 & \hat{\boldsymbol{n}}\cdot \boldsymbol{v}_{\rm {e}} \le 0 \\ 2 \bar{\rho } (\hat{\boldsymbol{n}}\cdot \boldsymbol{v}_{\rm {e}})^{2} & \hat{\boldsymbol{n}}\cdot \boldsymbol{v}_{\rm {e}} \gt 0 \end{array}\right. \end{eqnarray}$$

  (22)

  with the condition that the contact pressure must be non-negative. Here P, |$\hat{\boldsymbol{n}}$|⁠, ρ, c_s, are the pressure, normal to the face in the direction of the boundary cell, density, and sound speed. The bar denotes an average of the values on each side of the face, which are obtained from the state vectors |$\boldsymbol{\mathcal {U}}_{\rm {e}}$| and |$\boldsymbol{\mathcal {U}}_{\rm {b}}$|⁠. For a face with area A, this results in a hydrodynamic force

  $$\begin{eqnarray} \boldsymbol{F}_{\rm {hy}} = \hat{\boldsymbol{n}}A P_{\rm {con}} \end{eqnarray}$$

  (23)

  at the face. We sum all such forces at boundary faces to find the total hydrodynamic force on the boundary.

3.7 Gravity solver

manga uses a tree-based gravity solver originally implemented in gasoline, which uses multipole moments up to fourth order (hexadecapole) to represent the mass distribution within cells at each level of the tree (Wadsley, Stadel & Quinn 2004). After a gravity solution is performed, we intercept the gravitational accelerations |$\boldsymbol{a}_{i}$| on all boundary particles and average them to find the gravitational acceleration of the boundary |$\boldsymbol{a}_{\rm {gr}}$|⁠.

For all particles present regardless of type, manga uses spline softening to soften the gravitational forces within radius h of the particle. All gas cells in manga have equal softening lengths h_gas, defined at the start of the simulation as the mean separation between cell centres in the region of highest density

$$\begin{eqnarray} h_{\rm gas} = 2(m_{\rm {max}} / \rho _{\rm {max}})^{1/3}. \end{eqnarray}$$

(24)

Here, ρ_max is the maximum density and m_max is the cell mass at the location of ρ_max. We require boundary cells to have this same softening length, which guarantees that the self-gravity of the boundary does not impact its motion.

3.8 Motion of boundary

The motion of our boundaries is reactive to both hydrodynamic and gravitational forces. The hydrodynamic force is calculated during the half time-step Riemann solution and is the sum of all forces on the boundary faces (equation 23). For a boundary with mass m_b, this results in an acceleration

$$\begin{eqnarray} \boldsymbol{a}_{\rm {hy}} = \frac{1}{m_{\rm {b}}} \sum _{j} \hat{\boldsymbol{n}}_{j} A_{j} P_{{\rm con},j}. \end{eqnarray}$$

(25)

We combine this with |$\boldsymbol{a}_{\rm {gr}}$| to get the total acceleration of the boundary

$$\begin{eqnarray} \boldsymbol{a}_{\rm {b}} = \boldsymbol{a}_{\rm {gr}} + \boldsymbol{a}_{\rm {hy}}. \end{eqnarray}$$

(26)

At the drift step, the boundary velocity kick

$$\begin{eqnarray} \delta \boldsymbol{v}_{\rm {b}} = \boldsymbol{a}_{\rm {b}} \delta t \end{eqnarray}$$

(27)

is used to update the velocities of all boundary and edge cells. Unlike normal gas cells, the boundary and edge cells do not receive kicks in the direction of the fluid velocity, nor the velocity corrections described in CWQ17.

4 TEST PROBLEMS

4.1 Supersonic wind

We begin by simulating a sphere moving through air at supersonic speed. The air is initialized with atmospheric conditions ρ = 0.001225 g cm⁻³, T = 288 K, γ = 1.4, and mean molecular weight μ = 29 g/mol. Our simulation box has dimensions (10, 6, 6) cm and our spherical boundary has a radius of 1 cm. We use an HLLC Riemann solver for this test and do not include gravity. For all simulations presented in this paper, we use an adiabatic ideal gas equation of state.

We initialize the positions of our mesh-generating points using a pre-computed glass distribution of particles embedded in a 3D cube. This glass distribution is periodically replicated to produce a sufficient number of particles for our simulation. Low spatial resolution is problematic for this test case as large mesh cells cause the surface of the boundary to be rough and jagged, which can eventually lead to unphysical values for the state vector |$\boldsymbol{\mathcal {U}}$|⁠. In fact, we find in all test cases that the stability of the code depends largely on the smoothness of the boundary surface. When the boundary surface is well-resolved, we find that our code is robust against unphysical results, even in extreme situations. For this test case, our simulation contains 1.5 × 10⁶ gas and 1.7 × 10⁴ boundary cells.

Our boundary moves through the air with Mach number M = 2, shown in Fig. 1. The top left-hand panel shows a density projection at t = 15.5 μs, where we can see the formation of a bow shock. The bottom left-hand panel shows a density slice through the z = 0 plane overplotted with streamlines, where we see turning of the streamlines across the shock as well as a stagnation point behind the boundary. In the panels on the right, we show the same plots at t = 38.7 μs. Here, the shock has strengthened and we see a void in the path cleared by the sphere. In addition, the turbulence behind the boundary has now evolved into vortices. Qualitatively, these are the results that one would expect from this test case. This is encouraging, as the problem of an object moving through a dense gas appears in many astrophysical settings.

Figure 1.

The supersonic wind test case at t = 15.5 μs (left) and t = 38.7 μs (right). (Top) Density projections along the z-axis, showing the bow shock as it strengthens. (Bottom) Density slices through the z = 0 plane with streamlines. Here, the stagnation point and vortices are visible.

4.2 Sod shock tube

We now perform a Sod shock tube test using a simulation box of dimension (16, 1, 1) cm centred at (0, 0, 0). On the left side of the box (x < 0), the gas has density ρ = 1 g cm⁻³ and temperature T = 2.5 K, while the right side (x ≥ 0) has ρ = 0.25 g cm⁻³ and T = 1.8 K. The gas is initially stationary |$(\boldsymbol{v}= 0)$| with γ = 1.4. We use a periodic simulation box, so that the shock propagates through the tube in both directions. These conditions are similar to those used by CWQ17 to validate manga. We choose our moving boundary to be a piston of mass m = 0.1 g (1 per cent of the mass of the gas in the tube) that extends from x = 1 cm to x = 2 cm and spans the tube in the y- and z-directions. The piston is free to move but is initially at rest. We use a rectangular mesh – an appropriate choice because of our rectangular piston – with 160 cells in the x-direction and 10 cells in the y- and z-directions for a total of 16 000 cells.

In Fig. 2, we show the density as a function of position at t = 19.4, 42.6, and 85.2 μs. We initially see the formation of the shock front, contact discontinuity, and rarefaction wave. At t = 34.8 μs, the shock front impacts the piston (shaded in red) and reflects back to the left. This impact also pushes the piston to the right, inducing a transmitted shock.

Figure 2.

Gas density in our Sod shock tube test at t = 19.4, 42.6, and 85.2 μs. (Top) We initially see a shock front propagating to the right, followed by a contact discontinuity, as well as an expanding rarefaction wave. (Middle) The shock front then collides with the piston (shaded in red) and reflects back to the left. (Bottom) The piston is pushed to the right, creating a transmitted shock; the reflected shock passes the contact discontinuity.

We validate the motion of the piston using an analytical result for the impact of a shock wave on a movable wall. Following Meyer (1957), the motion of the wall is described by

$$\begin{eqnarray} \frac{du_{\rm w}}{dt} = \frac{p_{\rm T}-p_{\rm L}}{p_{\rm T 0}-p_{1}} \left(\frac{du_{\rm w}}{dt} \right)_{t=t_{\rm i}}. \end{eqnarray}$$

(28)

Here u_w is the wall velocity, p is the gas pressure, and t_i is the time of shock impact. The subscripts L and T correspond to the leading and trailing edges of the wall (Fig. 3). The regions 1 and 2 are the states on the right and left sides of the wall before disruption by the transmitted or reflected shocks. Furthermore, p_T0 is the pressure on the trailing edge just after shock reflection. We use the isentropic relations as well as the Riemann invariants

$$\begin{eqnarray} P = \frac{2}{\gamma -1}c_{\rm s}+v \quad \text{and}\quad Q = \frac{2}{\gamma -1}c_{\rm s}-v \end{eqnarray}$$

(29)

to obtain the relations

$$\begin{eqnarray} p_{\rm L} = p_{1} \left(\frac{Q_{1}+u_{\rm w}}{Q_{1}}\right)^{2\gamma / (\gamma -1)} \end{eqnarray}$$

(30)

and

$$\begin{eqnarray} p_{\rm T} = p_{\rm T 0} \left(\frac{P_{2}-u_{\rm w}}{P_{2}}\right)^{2\gamma / (\gamma -1)}. \end{eqnarray}$$

(31)

We can express the initial wall acceleration pressure as

$$\begin{eqnarray} \left(\frac{du_{\rm w}}{dt} \right)_{t=t_{\rm i}} = (p_{\rm T 0}-p_{1}) \frac{A}{m}, \end{eqnarray}$$

(32)

where A is the cross-sectional area of the wall. Putting this all into equation (28), we arrive at

$$\begin{eqnarray} \frac{du_{\rm w}}{dt} = \frac{A}{m} \left[ p_{\rm T 0} \left(\frac{P_{2}-u_{\rm w}}{P_{2}}\right)^{2\gamma / (\gamma -1)} - p_{1} \left(\frac{Q_{1}+u_{\rm w}}{Q_{1}}\right)^{2\gamma / (\gamma -1)} \right], \end{eqnarray}$$

(33)

which can be numerically integrated twice with the boundary conditions u_w = 0 and x = 1 cm at t = t_i to find the path of the wall. We obtain the relevant pressures and Riemann invariants for equation (33) from our simulation output. The results are shown in Fig. 4 alongside those obtained in manga, showing that the two are consistent.

Figure 3.

Wave diagram showing the paths of the wall and shock waves. The regions 1 and 2, which denote the states on the right and left sides of the wall before disruption by the transmitted or reflected shocks, are shown as well as the pressures relevant to equation (28).

Figure 4.

An analytical solution for the response of a movable wall to a shock wave (blue) following Meyer (1957) compared with the results from our Sod shock tube test case in manga (red). On the horizontal axis is the time after shock impact t_i.

4.3 Sedov–Taylor blast wave

In this test problem, we use the same setup as in the Sod shock tube test. However, we now initialize the gas with density ρ = 1 g cm⁻³ and temperature T = 1 K except for the region |x| < 0.1 cm, where we set T = 1000 K. This results in a point explosion that is described by the Sedov–Taylor solution. It is known that the Sedov–Taylor solution is self-similar in 3D, and we briefly show that this also holds in 1D. Following Taylor (1950), the appropriate similarity assumptions for an expanding blast wave of constant total energy in 1D are

$$\begin{eqnarray} p/p_{0} = R^{-1} f_{1}, \quad \rho / \rho _{0} = \psi \quad \text{and}\quad u = R^{-1/2} \phi _{1}. \end{eqnarray}$$

(34)

Here, R is the radius of the blast wave, u is the fluid velocity, p₀ and ρ₀ are the ambient pressure and density, and f₁, ψ, and ϕ₁ are functions only of similarity variable η = x/R. In the absence of gravity, we can now cast the Euler equations (1), (2), and (3) into the 1D forms

$$\begin{eqnarray} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = -\frac{p_{0}}{\rho } \frac{\partial (R^{-1} f_{1})}{\partial x}, \end{eqnarray}$$

(35)

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + u \frac{\partial \rho }{\partial x} + \rho \frac{\partial u}{\partial x} = 0 , \end{eqnarray}$$

(36)

and

$$\begin{eqnarray} \left(\frac{\partial }{\partial t} + u\frac{\partial }{\partial x}\right)(p \rho ^{-\gamma }) = 0. \end{eqnarray}$$

(37)

These equations are satisfied if |$\frac{dR}{dt} = A R^{-1/2}$|⁠, where A is a constant. If we additionally define a² = γp₀/ρ₀, f = f₁a²/A², and ϕ = ϕ₁/A, then equations (35), (36), and (37) become

$$\begin{eqnarray} \frac{1}{2}\phi - \eta \phi ^{\prime } + \phi \phi ^{\prime } = \frac{1}{\gamma \psi }(f - f^{\prime }), \end{eqnarray}$$

(38)

$$\begin{eqnarray} -\eta \psi ^{\prime } + \phi \psi ^{\prime } + \phi ^{\prime }\psi = 0 , \end{eqnarray}$$

(39)

and

$$\begin{eqnarray} -\gamma \eta f \psi ^{\prime } + f \psi ^{-\gamma } + \eta f^{\prime } \psi ^{-\gamma } + \gamma f \phi \psi ^{\prime } - f^{\prime } \phi \psi ^{-\gamma } = 0, \end{eqnarray}$$

(40)

respectively. Here, a prime denotes a derivative with respect to η. We see that equations (38), (39), and (40) are independent of R; they depend only on similarity variable η and the gas properties p₀, ρ₀ and γ. Therefore, this is a self-similar solution. Furthermore, because equations (38), (39), and (40) are all first-order, only one boundary condition is needed for each of f, ϕ, and ψ to find a solution. These are given by the Rankine–Hugoniot shock jump conditions, which for a sufficiently strong shock take on the asymptotic forms

$$\begin{eqnarray} \psi = \frac{\gamma +1}{\gamma -1}, \quad f = \frac{2\gamma }{\gamma +1}, \quad \text{and}\quad \phi = \frac{2}{\gamma +1} \end{eqnarray}$$

(41)

at the shock front (η = 1).

We again choose a piston as our moving boundary, extending from x = 1 to x = 1.5 cm with mass m = 0.5 g. We plot the density as a function of position in Fig. 5 at t = 6 μs, showing only particles with |y| < 0.1 and |z| < 0.1 cm in order to reduce the particle noise. The blast wave quickly forms and propagates outward until it impacts the piston, which is shaded in red, pushing it to the right and creating a transmitted shock. A gap is visible between the left side of the piston and the mesh-generating points comprising the shock front, due to the edge cell centres maintaining their original distance from the boundary. To analyse the conservation properties of our simulation, we consider the fractional error in the conservation of momentum of the system ϵ_p, including the contributions from both the gas and the piston. Here, we find ϵ_p = 4.68 × 10⁻¹⁰. However, we also see from the spacing of the mesh-generating points that the spatial resolution of the blast wave is poor and thus conduct a resolution test to determine its effect on ϵ_p. Rather than changing the cell size, we exploit the self-similarity of the Sedov–Taylor solution. We perform two additional runs in which the piston is initially placed further to the right so that the number of cells N_c resolving the blast wave is greater. The mass and width of the piston are taken to be proportional to R to preserve the self-similarity of the system. Specifically, we use pistons of mass m = 1 and m = 1.5 g, extending from x = 2 to x = 3 and from x = 3 to x = 4.5 cm, respectively. These are shown in Table 1 along with ϵ_p, which we calculate when the left edge of the piston has moved from its initial position by 10 per cent. We find that higher resolution improves momentum conservation as |$\epsilon _{p}\propto N_{\rm c}^{-1.7}$|⁠.

Figure 5.

Gas density in the Sedov–Taylor blast wave test at t = 6 μs, restricted to particles with |y| < 0.1 and |z| < 0.1 cm to reduce the particle noise. The blast wave is colliding with a piston (shaded in red) and pushing it to the right. The gap between the left side of the piston and the mesh-generating points comprising the shock front is due to the edge cell centres maintaining their original distance from the boundary.

Table 1.

Piston initial position, piston mass m, number of cells resolving the blast wave N_c, and momentum error ϵ_p in our Sedov–Taylor blast wave resolution test. We calculate ϵ_p when the left edge of the piston has moved from its initial position by 10 per cent.

x Position (cm)	m (g)	Nc	ϵp
1–1.5	0.5	10	4.68 × 10−10
2–3	1	20	1.27 × 10−10
3–4.5	1.5	30	7.44 × 10−11

4.4 Common envelope evolution

In this test problem, we demonstrate the applicability of moving boundaries to astrophysical systems by simulating CEE, a brief phase in the life of a binary star system. In CEE, two stars – a giant star and a much smaller companion – share a common gaseous envelope (for a review, see Ivanova et al. 2013). As the companion and the core of the giant orbit each other, loss of orbital energy and angular momentum drive the two stellar centres close to one another and eject the envelope from the system. CEE is a critical process in the lives of these binary systems and is responsible for progenitors of Type Ia supernova (potentially), X-ray binaries, double white dwarfs, double neutron stars, and possibly the merging double black hole and double neutron star systems discovered by Advanced LIGO (Ivanova et al. 2013; Belczynski et al. 2016).

In particular, multiple groups have run 3D numerical simulations of CEE. Several numerical methods have been tried, including SPH (Passy et al. 2012; Nandez, Ivanova & Lombardi 2015), Eulerian grid codes (Passy et al. 2012; Ricker & Taam 2012; Chamandy et al. 2019a), and moving-mesh codes (PC19 Ohlmann et al. 2016). In all of these simulations, the radius of the companion is considered negligibly small compared to the size of the giant, and thus the companion is modelled by a point particle that interacts only through gravity. This approximation is accurate when the companion in question is a compact object such as a black hole or neutron star; however, in many systems of interest the companion is a small star with size of order 1 R_⊙. In such cases, the drag force resulting from the companion encountering the envelope may significantly alter the dynamics of the spiral-in. This is where moving boundaries come into play, as the companion can be modelled as a spherical boundary. Chamandy et al. (2019b) and De et al. (2019) investigated the drag force on the companion in CEE simulations using adaptive mesh refinement codes with point particle companions, and in future work, we hope to make comparisons to their results.

The application of manga to CEE was the main focus of PC19, and here we use the same binary system. We summarize the simulation setup here but refer the reader to PC19 for a more detailed discussion. We first use mesa (Paxton et al. 2011, 2013, 2015, 2018, 2019) to evolve a 2 M_⊙ star with metallicity Z = 0.02 from the pre-main sequence to the red giant phase. We stop when the stellar radius reaches 52 R_⊙ with a helium core of mass M_c = 0.36 M_⊙. Because of the great difference in density between the helium core and the hydrogen envelope, we model the core as a dark matter particle with gravitational softening length R_c = 1.99 R_⊙.

We construct an appropriate particle mesh for the star using the same glass distribution as in Section 4.1. We assume that each particle is of approximately equal volume and re-scale them to the appropriate mass based on the computed M(r) from mesa. The total number of particles representing the star is 8.5 × 10⁴. Outside of the star, we include a low-density atmosphere of 10⁻¹³ g cm⁻³ with temperature 10⁵ K that extends out to the total box size of 3.5 × 10¹⁴ cm (5000 R_⊙), with periodic boundary conditions at its edges. The total number of particles within the simulation box is 5.4 × 10⁵.

We place the 1 M_⊙ companion in an initially circular Keplerian orbit at the red giant radius a = 52 R_⊙. We assume that, prior to the onset of CEE, tidal forces from the companion have spun up the rotation of the giant (Soker 1996). Thus, we set the rotation of the giant such that it is tidally locked to the binary orbit.

To create the moving boundary, we assume that the companion is a sphere of radius 4 R_⊙, which is unrealistically large. We use such a large companion as a proof of concept for the methodology of moving boundaries. As discussed in Section 3.2, it is important to have a sufficiently high density of mesh-generating points near the boundary to accurately produce the desired boundary shape. To this end, we refine our mesh near the initial position of our boundary. Note that this needs only be done during the initialization; because the boundary cells and edge cells move identically, a boundary that is highly resolved initially will remain so indefinitely. We perform our refinement within a cube of side length 10 R_⊙ centred on the companion. We use the same glass distribution but scaled such that the cube contains 1000 mesh-generating points. Our boundary comprises 270 of these points, and the rest of the particles in the cube are given the same properties as the external atmosphere. This somewhat cuts into the red giant envelope; however, the outer layers of the envelope are very rare and the total mass of the giant is not significantly altered.

We simulate the binary for 33 d and show several density projections in Fig. 6 at t = 1, 5, 10, 15, 20, and 30 d. For the first 14 d, the companion plunges into the envelope, reducing the separation between the core and the companion by a factor of 5 (Fig. 7). We see tidal tails of gas ejected from the system by both the core and the companion during this plunge. The orbit then continues to shrink for the remainder of the simulation period as orbital energy is transferred to the gas. Spiral shocks facilitating this transfer can be seen in the lower panels of Fig. 6. The simulation ends well before the outflow reaches the edge of the simulation box.

Figure 6.

Density projections along the z-axis in our CEE simulation. The + sign marks the red giant core.

Figure 7.

Separation a between the stellar core and companion. The companion quickly plunges into the envelope, reducing the separation between the core and the companion by a factor of 5. The orbit then continues to shrink at a slower rate.

The total energy of each gas particle is given by

$$\begin{eqnarray} E_{\rm tot,i}= m_{i} \left[ \frac{1}{2} (\boldsymbol{v}_{i} - \boldsymbol{v}_{\rm CM})^{2} + \phi _{i} + I_{\rm {e},\mathit{ i}} \right] \end{eqnarray}$$

(42)

(Nandez et al. 2015), where m_i, |$\boldsymbol{v}_{i}$|⁠, |$\boldsymbol{v}_{\rm CM}$|⁠, ϕ_i, and I_{e, i} are the mass, velocity, centre of mass (CM) velocity of the bound material, gravitational potential, and specific internal energy of each particle. Gas particles with a negative total energy are bound to the binary, while those with a positive total energy are unbound. The kinetic energy is computed relative to the velocity of the CM of the bound matter, that is, the bound gas as well as the red giant core and companion. However, because the total energy is needed in order to determine which gas particles are bound, we use an iterative scheme to find the velocity of the CM, which is described in PC19. With the total energies of all particles known, we can find the mass of the unbound gas m_{g, u} as a fraction of the total mass of the envelope m_env,

$$\begin{eqnarray} f_{\rm unb}= \frac{ m_{\rm g,u}}{ m_{\rm env} }. \end{eqnarray}$$

(43)

We refer to this as the ejection efficiency, shown in Fig. 8. The ejection efficiency has reached nearly 20 per cent by the end of the simulation period but is still steadily increasing. The first periastron passage at t ≈ 14 d temporarily causes an increase in f_unb.

Figure 8.

The fraction of mass of the envelope that has acquired enough energy to be unbound from the system. It shows fairly linear behaviour over most of the simulation period, with the exception of the first periastron passage at t ≈ 14 d.

This simulation exhibits features similar to those shown in the previous work listed above but with a lower resolution and a much shorter simulation period. Therefore, the numerical results presented here should be considered only as a validation of our methodology. Nevertheless, this test case shows that using moving and reactive boundaries in CEE simulations is feasible and produces realistic results. We leave a detailed study for future work.

5 DISCUSSION AND FUTURE DIRECTIONS

We describe the implementation of moving and reactive boundary conditions into the moving-mesh code manga. ALE schemes are ideal for moving boundaries because the mesh cells are designed to move with respect to one another; S10 demonstrates this by carrying out a simple 2D simulation of a curved boundary that moves through a fluid along a fixed path. We describe the hydrodynamics of boundary conditions, which can be realized by modification of the flux function |$\boldsymbol{\mathcal {F}}$| and state vector |$\boldsymbol{\mathcal {U}}$| at the boundary surface in both the Riemann solver and the gradient estimation. To realize reflective boundary conditions, the state vector and flux function on either side of a boundary face must be equal except for the fluid velocities normal to the face, which must be equal and opposite. For outflow conditions, the fluid velocities are simply equal. We show that all components of the HLL flux across a boundary face are zero except for the momentum normal to the face, meaning that no other conserved variable can be transported across the boundary. We require that the gas cells adjacent to the boundary do not move with respect to the boundary, so that the surface of the boundary does not deform. However, the motion of our boundaries is reactive to external forces. We describe the modifications to our Riemann and gravity solvers to compute the aerodynamic drag and gravitational force on the boundary.

We test our methodology in four cases. The first is a supersonic spherical boundary with Mach number M = 2 moving through air, forming a bow shock on the leading edge. We see the expected turbulence pattern behind the sphere: a stagnation point appears before evolving into vortices. The next two are the Sod shock tube test and a Sedov–Taylor blast wave, where in both cases the shock front collides with a movable piston. We validate the motion of the piston in the Sod shock tube case by comparison with the analytical results of Meyer (1957). For the case of a blast wave, we show that the Sedov–Taylor solution is self-similar in one dimension and exploit this fact to find a relationship between the number of cells used to resolve the blast wave and the fractional momentum error. Finally, we demonstrate an astrophysical application of our methodology by simulating CEE in a binary system. All previous work on numerical simulations of CEE have modelled the companion object as a point particle, and moving boundaries allow us to model it as an object of finite size. Though this test simulation has low spatial resolution, it shows the expected behaviour of a common envelope phase.

We find that the stability of the code depends largely on the smoothness of the boundary surface. This typically necessitates a high density of mesh-generating points in the vicinity of the boundary at the start of the simulation. In our CEE simulation, we refine the mesh in a cubical region surrounding the initial position of our boundary. Because the number of cells in this cube after the refinement is small compared to the total number of cells, this does not significantly slow our code.

Moving boundaries have yet to find success in astrophysics due to the difficulty of implementation in most codes as well as the relative lack of applications as compared to other fields. However, there are many cases in which they prove useful, including the study of stellar pulsations as well as any astrophysical object moving through a dense gas. The first and most straightforward use for this methodology will be to conduct in-depth simulations of CEE using a moving boundary. High-resolution simulations of CEE using the mesa equation of state could shed light on the effects of companion size and inform future work. In particular, the drag force on the companion is expected to depend upon the size of the companion and could alter the dynamics of the spiral-in phase (Chamandy et al. 2019b; De et al. 2019).

Furthermore, if the creation and destruction of mesh-generating points is implemented into manga, additional applications will become possible. Because our reflective boundary conditions can be transformed into outflow conditions with a few sign changes in the code, our methodology can also be used for outflow surfaces such as black holes. Although black holes are negligibly small in many astrophysical systems, their size can be significant in cases such as the tidal disruption of a star by a supermassive black hole (SMBH). If the SMBH is modelled as a point particle, it can accrete gas that exerts an (unphysical) outward pressure (Spaulding & Chang, in preparation). Modelling the SMBH using an outflow boundary avoids this problem as the accreted gas is removed from the simulation entirely.

ACKNOWLEDGEMENTS

We thank Philip Chang for useful and detailed discussions. We are supported by the NASA Astrophysical Theory Program (ATP) through NASA grant NNH17ZDA001N-ATP. We used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (NSF) grant no. ACI-1053575. We acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper (URL: http://www.tacc.utexas.edu). We use the yt software platform for the analysis of the data and generation of plots in this work (Turk et al. 2011).

Computational resources supporting this work were also provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at the Ames Research Center. This material is based upon work supported by NASA under Award No. RFP19_7.0 issued through the Wisconsin Space Grant Consortium and the National Space Grant College and Fellowship Program. Any opinions, findings and conclusions, or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Aeronautics and Space Administration.