2.1 Basic assumptions

The physical system considered consists of four subsystems, representing the different materials involved. Its first part represents its boundary, formed by the perfectly rigid, impervious, axisymmetric, relatively slender bottleneck, typically made of glass. The accordingly impervious axisymmetric stopper made of cork forms the second. The third is given by the expanding gas, a solution of ${\rm CO}_2$, slowly degassed by the liquid champagne and so accounting for its pressurisation prior to the opening of the bottle, air, alcohol and water vapour, which propels the stopper out of the bottle in the axial direction. It is confined by the surface of the contained champagne, during the process considered taken as a stationary, plane and impermeable boundary, and before the opening sealed by the stopper. For what follows, the bottleneck is then precisely defined as the void portion of the bottle between its opening and the liquid level. The initially quiescent ambient air represents the fourth subsystem. This, and the expanding gas, are taken as Newtonian ideal gases having constant specific heats, heat conductivities and dynamic viscosities (perfect gas), all of negligibly different magnitudes. Correspondingly, they form non-distinguishable phases, the mixing of which is entirely ignored, as are the temperature dependences of those quantities. Even though reasonable for the gas phase, this simplification admittedly disregards the phase transitions observed by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019).

The actual opening of the bottle effectively means the rapid overpowering of stiction (self-retention) between the stopper and the bottleneck, compressing the former radially, in favour of much lower sliding friction setting in. The so triggered force imbalance has the gas expand and starts the proper motion of the stopper. Here, this very first stage of the opening process is abstracted as instantaneous with respect to the resolved time scales. It thus provides the initial condition. Hence, the stopper is initially accelerated by the gas but decelerated by wall friction. The expanding gas remains contained in the bottle as long as the stopper has not entirely left the bottle. Simultaneously, it starts to move the outer air. Once it has left the orifice behind, it rapidly expands towards its original (decompressed) shape as it is pushed further by the now forming (under-expanded) gas jet and the thereby generated air flow. The gas and the air completely dissolve into one another, and the resulting flow and the motion of the stopper are presumed to remain axisymmetric.

We neglect the characteristic ‘dome’ of the stopper for a champagne bottle, which, in combination with the frustum at its base, gives it the characteristic mushroom-like appearance. However, this renders its discretisation unnecessarily complex whilst taking into account that the drag-increasing bumper would not alter substantially the fascinating flow patterns between the stopper and the bottle opening, tied in with the emergence of the Mach disc and fading out once the stopper is sufficiently remote from the opening. That surprisingly complex shock regime lies at the heart of our numerical analysis.

Before we specify the modelling of the four subsystems and their interactions, we summarise the essential assumptions and simplifications (i)–(vi):

1. (i) The ambient air and the pressurised gas mixture are indistinguishable single-phase ideal gases, characterised by constant (temperature- and pressure-independent) thermophysical properties (hence, the pressurised gas treated as air).

2. (ii) The bottleneck is beheld as a slender truncated cone, tapered towards the bottle opening, with sufficient accuracy. The mushroom-shaped cap of the stopper and the fastening wire cage (muselage) are both disregarded (as being irrelevant for the process under consideration).

3. (iii) The expanded (unloaded) stopper is also taken as a slender truncated cone, even though real sparkling wine corks are so-called agglomerated ones, composed of up to three rings of slightly different radii (Reference MargalitMargalit, 2012).

4. (iv) The cork used is typically modelled as a homogeneous and isotropic hyperelastic material of negligibly small lateral contraction: for this special case of an Ogden material model see Reference Sasso, Lattanzi, Farotti, Sarasini, Sergi, Tirillò and ManciniSasso et al. (2022), Reference Fernandes, Pascoal and Alves de SousaFernandes, Pascoal, and Alves de Sousa (2014) and references therein as well as Reference OgdenOgden (1997) and Reference DillDill (2006). Here, both the weak visco-elasticity usually found to be at play and the degradation of elasticity due to long-time ageing are ignored. Since the compressive reactive force the weakly tapered bottleneck exerts onto the stopper sliding along it appears to be much larger than all the other forces at play, the stopper can deform only radially while it remains rigid in the axial direction. Therefore, it assumes its original undeformed shape once its reversible decompression has stopped after a finite, properly defined relaxation time.

5. (v) We model the instantaneous expansion of the stopper by taking it as cylindrical when it has just left the bottle and on the basis of the estimated characteristic time scale at play rather than founding it more rigorously. This seems admissible given the relatively short duration of the expansion process.

6. (vi) We neglect any (unavoidable) leakage of gas out of the bottle as long as the stopper has not completely been released from it. That is, we may assume dry sliding (Coulomb) friction between the cork and the glass of the bottle (see Reference Fatima Vaz and FortesFatima Vaz & Fortes, 1998).

Furthermore, the following conventions prove sensible. Tildes indicate dimensional quantities while otherwise their non-dimensional form is used; the subscript $0$ indicates the reference state of the quiescent ambient air; the subscripts $B$ and $C$ indicate the properties of, respectively, the bottleneck, including the pressurised gas initially at rest and provided by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019), and the cork stopper.

2.2 Geometry of bottleneck, unloaded stopper and computational domain

We introduce cylindrical coordinates $\tilde {r}$ and $\tilde {z}$, radially from and along the axis of the bottle outwards from its opening, respectively. Let us first describe the two truncated cones representing the bottleneck and the fully relaxed stopper as inferred from figure 1. The bottleneck of an axial extent $\skew{2}\tilde {l}_B$ is slightly tapered under an inclination angle $\theta _B$ towards its opening, having a diameter $\skew{6}\tilde {d}_0$; this defines its void volume $\tilde {V}_B$. The accordingly shaped relaxed stopper of an axial length $\skew{2}\tilde {l}_C$ has an angle of taper $\theta _C$ and the diameters $\skew{6}\tilde {d}_2$ at its base (rear) end and $\skew{6}\tilde {d}_1(<\skew{6}\tilde {d}_2)$ at its top (front) end, this initially exposed to the ambient air. The (average) thickness $\Delta \tilde {r}_B$ of the bottle glass completes the geometry of the systems bottleneck and stopper.

Hereafter, all lengths are advantageously made non-dimensional with $\skew{6}\tilde {d}_0$. At first, $r=\tilde {r}/\skew{6}\tilde {d}_0$, $z=\tilde {z}/\skew{6}\tilde {d}_0$, $d_{1,2}=\skew{6}\tilde {d}_{1,2}/\skew{6}\tilde {d}_0$, $l_B=\skew{2}\tilde {l}_B/\skew{6}\tilde {d}_0$. Using $V_B=\tilde {V}_B/\skew{6}\tilde {d}_0^3$, we then introduce the (small) slopes

(2.1a,b)\begin{equation} a_B = \tan{\theta_B} = \frac{1}{l_B} \left(\sqrt{\frac{3 V_B}{{\rm \pi} l_B}-\frac{3}{16}}-\frac{3}{4}\right),\quad a_C = \tan{\theta_C} = \frac{d_2-d_1}{2 l_C}.\end{equation}

The first relationship herein is needed to complete table 1, presenting the geometrical input data together with a mean value $\tilde {\rho }_C^m$ of the density of a fully relaxed and dry cork (Reference SaadallahSaadallah, 2020). We obtained the geometrical data by measuring customary sparkling-wine bottles and the associated stoppers, apart from the pairing for the larger value of $\tilde {V}_B$, which we could not verify. However, only this and the value of $\skew{6}\tilde {d}_0$ were proposed by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019). We extracted the value of $\theta _B$ from the measured one of $\skew{2}\tilde {l}_B$ for $\tilde {V}_B=20$ ml and (2.1a). In order to accomplish the task of comparing our simulations with the experiments, we then stipulated a congruent bottleneck and thus the same value of $a_B$ for $\tilde {V}_B=25$ ml. The different liquid level $\skew{2}\tilde {l}_B$ then ensues as the single real root of the cubic satisfied by $l_B$ and arising from (2.1a). We also note the so found inclination angles $\theta _B\doteq 2.24^\circ$ and $\theta _C\doteq 3.43^\circ$.

Table 1.

Geometrical input quantities, cf. figure 1, and density of relaxed cork (all suitably rounded).

The axisymmetric computational domain is given by $0\leq r\leq 2$ and $0\leq z+l_B\leq 8$. Numerical tests with enlarged domains, at the expense of considerably increased computational costs, suggest the chosen dimensions as being sufficiently large (cf. the discussions regarding the spatial resolution in § 4).

2.3 Gas/air flow

According to the premise (i) in § 2.1, we introduce the dynamic viscosity $\tilde {\eta }$, the heat conductivity $\tilde {\lambda }$, the mass-specific gas constant $\tilde {R}$ and the heat capacity ratio $\kappa$ for air under standard conditions as the required, uniform gas properties. Let the ambient air be initially at rest under the uniform pressure $\tilde {p}_0$ and the temperature $\tilde {T}_0$. From the common relations holding for an ideal gas, $\tilde {c}_{p}=\tilde {R}\kappa /(\kappa -1)$ is its specific heat capacity at constant pressure, $\tilde {\rho }_0=\tilde {p}_0/(\tilde {R}\tilde {T}_0)$ its density at rest, $\tilde {c}_0=(\kappa \tilde {R}\tilde {T}_0)^{1/2}$ the associated isentropic speed of sound and $\tilde {e}_0=\tilde {\rho }_0 \tilde {c}_{p}\tilde {T}_0/\kappa =\tilde {p}_0/(\kappa -1)$ the associated density of the internal energy. We then typically have $\tilde {R}\doteq 287.058$ J (kg K)$^{-1}$ and $\kappa \doteq 1.4$. These data, together with the characteristic, fixed inner diameter $\skew{6}\tilde {d}_0$ of the cross-section of the bottle opening, adopted by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019), establish the reference values for all the remaining input quantities. As the data put forward by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019) and Reference Benidar, Georges, Kulkarni, Cordier and Liger-BelairBenidar et al. (2022) shall validate our simulation results, we only consider the two combinations of the pressure $\tilde {p}_B$ and the temperature $\tilde {T}_B$ measured by Reference Liger-Belair, Cordier and GeorgesLiger-Belair et al. (2019) and that identify the state of contained gas initially at rest. All the essential input data and their values are summarised in table 2 and the resulting four combinations addressed by our simulations in table 3.

Table 2.

Thermophysical reference data (suitably rounded, slashes separate values of $\tilde {p}_B$ and $\tilde {T}_B$).

Table 3.

The combinations (A)–(D) of $\tilde {V}_B$ and $\tilde {T}_B$ used in the simulations.

Consequently, we consider the local values of five flow quantities, depending on space and time: the fluid velocity $\boldsymbol {v}$, the fluid density $\rho$, fluid pressure $p$, fluid temperature $T$ and the density of its total energy (the sum of internal and kinetic energy)

(2.2)\begin{equation} e_t = p + \kappa(\kappa-1)\rho|\boldsymbol{v}|^2 /2 ,\end{equation}

already made non-dimensional with, respectively, $\tilde {c}_0$, $\tilde {\rho }_0$, $\tilde {p}_0$, $\tilde {T}_0$ and $\tilde {e}_0$. Accordingly, the time $t$ is in natural manner made non-dimensional with the basic time scale $\skew{6}\tilde {d}_0/\tilde {c}_0\doteq 0.0524$ ms. As seen from the particular dimensionless form of $e_t$ in (2.2), $p$ also ensues as the density of the internal energy. Those flow quantities satisfy the thermal equation of state for an ideal gas

(2.3)\begin{equation} p = \rho T,\end{equation}

and the resulting full set of the non-dimensional Navier–Stokes equations. These describe conservation of mass (2.4), of momentum (2.5) and of the total energy density $e_t$, (2.6), all conveniently written in coordinate-free divergence form:

(2.4)\begin{gather} \partial_t\rho + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{v}) = 0, \end{gather}

(2.5)\begin{gather}\partial_t(\rho\boldsymbol{v}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho\boldsymbol{v}\boldsymbol{v}) = {-}\kappa^{{-}1}\boldsymbol{\nabla} p + {{Fr}}^{{-}2}\rho \boldsymbol{e}_g + {{Re}}^{{-}1}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varSigma}, \end{gather}

(2.6)\begin{gather}\partial_t e_t + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{v}e_t) = \kappa(\kappa-1) (-\kappa^{{-}1}\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{v}p) + {{Fr}}^{{-}2}\rho \boldsymbol{e}_g \boldsymbol{\cdot} \boldsymbol{v} + {{Re}}^{{-}1}\boldsymbol{\nabla} \boldsymbol{\cdot} [\boldsymbol{\varSigma} \boldsymbol{\cdot} \boldsymbol{v} + (\kappa-1)^{{-}1}{{Pr}}^{{-}1}\boldsymbol{\nabla} T]). \end{gather}

Herein, $\boldsymbol {\varSigma }$ denotes the tensor of the viscous Cauchy stresses for a Newtonian fluid and

(2.7a–c)\begin{equation} {{Fr}} = \tilde{c}_0/\sqrt{\tilde{g}\tilde{d}_0}\doteq 817, \quad {{Re}} = \tilde{c}_0\tilde{d}_0\tilde{\rho}_0/\tilde{\eta}\doteq 4.07\times 10^5, \quad {{Pr}} = \tilde{\eta}\tilde{c}_{p}/\tilde{\lambda}\doteq 0.701,\end{equation}

respectively, the Froude, Reynolds and Prandtl numbers as the key groups at play. In the first, $\tilde {g}$ denotes the (constant) scalar gravitational acceleration, $\tilde {g}\doteq 9.81\ {\rm m}\ {\rm s}^{-2}$, acting in the direction of some (constant) unit vector $\boldsymbol {e}_g$. The figures of these three parameters then follow from the aforementioned input values (see tables 1 and 2). The set (2.2)–(2.6) governs the aforementioned five flow quantities in full, apart from the required boundary, coupling and initial conditions discussed below.

Most importantly, ${{Fr}}$ and ${{Re}}$ appear to be so (predictably) large that we may safely consider the limits ${{Fr}}\to \infty$ and ${{Re}}\to \infty$ in (2.5) and (2.6). Then (2.4)–(2.6) reduce to the Euler equations governing the adiabatic and inviscid flow of a weightless ideal gas. To specify them for our axisymmetric situation, we introduce the $r$- and $z$-components $u$ and $w$, respectively, of $\boldsymbol {v}$. Then (2.4), the two scalar momentum equations for the $r$-and the $z$-direction resulting from (2.5) and the energy equation (2.6) give, respectively,

(2.8)\begin{gather} \partial_t\rho + \partial_r(\rho u) + \partial_z(\rho w) = {-}\rho u/r, \end{gather}

(2.9)\begin{gather}\partial_t(\rho u) + \partial_r(\rho u^2 + p/\kappa) + \partial_z(\rho u w) = {-}\rho u^2/r, \end{gather}

(2.10)\begin{gather}\partial_t(\rho w) + \partial_z(\rho w^2 + p/\kappa) + \partial_r(\rho u w) = {-}\rho u w/r, \end{gather}

(2.11)\begin{gather}\partial_t e_t + \partial_r(u[e_t + (\kappa-1)p]) + \partial_z(w([e_t + (\kappa-1)p]) = {-}u[e_t + (\kappa - 1)p]/r. \end{gather}

The source terms on the right sides of (2.8)–(2.11) are typically due to continuity in the radial direction. Hence, special care is required for resolving the flow near the axis of symmetry $r=0$. We are interested in finding $q(z,r,t)$, where $q$ stands for any of the dependent variables $\rho$, $\rho u$, $\rho w$, $e_t$ and $p$. Eliminating $p$ with the aid of (2.2) closes the system of (2.8)–(2.11). Specifically, $e_t+(\kappa - 1)p$ in (2.11) is replaced by $\kappa [e_t-(\kappa -1)^2\rho |\boldsymbol {v}|^2 /2]$. Governing the four remaining quantities, these equations are then ripe for their proper numerical treatment. Most importantly, their divergence form, together with (2.2), allows for their weak solutions, i.e. the capturing of shock waves. The temperature field $T$, the local speed of sound $c=T^{1/2}$, made non-dimensional with $\tilde {c}_0$, and the local Mach number $M=(u^2+w^2)^{1/2} /c$ can be calculated a posteriori, using (2.3). It is noted that $\kappa$ is the only material-specific parameter entering the problem.

The inviscid-flow description reduces the kinematic boundary conditions to those of axial symmetry and the impermeability of the gas/liquid interface and the bottle wall:

(2.12a,b)\begin{equation} u(z,0,t) = w({-}l_B,r,t) = 0\quad (0< r<1/2), \quad \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}_B = 0 \quad \text{on the bottle surface},\end{equation}

where $\boldsymbol {n}_B$ denotes a unit normal of the bottle wall at its respective position. The typical outflow conditions prescribed on the remaining (outer) boundaries of the computational domain complete the boundary conditions. We remark that the neglect of heat conduction and dissipation and the associated boundary layers renders the temperatures at the solid surfaces as just resulting from convection. Therefore, the (questionable) insulation conditions employed by Reference Benidar, Georges, Kulkarni, Cordier and Liger-BelairBenidar et al. (2022) are irrelevant in our inviscid approach.

The initial conditions are posed at some $t=t_0$ as we choose $t=0$ advantageously as that instance of time when the stopper has just entirely passed the bottle opening. That is, $z<0$ and $t<0$ indicate the regime of the sliding/compressed stopper that still seals the bottle, and $z>0$ and $t>0$ that of the freely moving released stopper and expanding jet. Hence, for

(2.13)\begin{equation} t = t_0({<}0){\colon}\quad u = w = 0,\quad \begin{cases} \rho = p = e_t = 1 & \text{(ambient air)}, \\ \rho = \tilde{\rho}_B/\tilde{\rho}_0,\quad p = e_t = \tilde{p}_B/\tilde{p}_0 & (\text{pressurised gas}), \end{cases}\end{equation}

cf. (2.2). The particular value of $t_0$ is extracted from the output data of the calculation initialised by (2.13).

2.4 Stopper–bottleneck and gas–stopper interactions

We first introduce the radial position $r=R(z,t)$ of the surface of the stopper, the axial one $z=Z(t)$ subject to $Z(0)=0$ at its base and the radial ones of the bottleneck, $r=r_B(z)$, and of the relaxed stopper in our resting frame of reference, $r=r_C(z-Z)$. From (2.1a,b)

(2.14a,b)\begin{equation} r_B(z) = 1/2-a_B\,z, \quad r_C(z-Z(t)) = d_2/2-a_C(z-Z(t)).\end{equation}

The fluid–stopper interaction is facilitated by virtue of the kinematic coupling condition and the dynamic one given by the equation of axial motion of the stopper:

(2.15)\begin{gather} ({\boldsymbol{v}}-{\boldsymbol{v}}_C)\boldsymbol{\cdot}{\boldsymbol{n}}_C = 0\quad \mbox{on the stopper's surface}, \end{gather}

(2.16a,b)\begin{gather}\kappa m_C\ddot{Z} = (F_b-F_t-F_B-F_{ls})(t),\quad m_C = \tilde{\rho}_C^m\tilde{V}_B/(\tilde{\rho}_0 \tilde{d}_0^3). \end{gather}

Here, $\boldsymbol {v}_C$ denotes the velocity on and $\boldsymbol {n}_C$ a unit normal to the stopper surface at the position considered, $m_C$ the non-dimensional mass of the stopper (cf. table 1) and $F_b$, $F_t$, $F_{ls}$ and $F_B$ are, respectively, the scalar pressure forces at the base, top and the axial components of the forces the fluid and the bottleneck exert at the respective portions of its lateral surface. We term $F_B$ the bottle force as it comprises the $z$-components of the compressive normal force and the tangential Coulomb friction force at play. Given the tapering of the bottleneck and the stopper, see figure 1, all forces in (2.16a) are taken as positive first.

The forces $F_b$ and $F_t$ have the straightforward explicit representations

(2.17a,b)\begin{equation} F_b(t) = 2{\rm \pi}\int_0^{R(Z(t),t)} p(Z(t),r,t)r\,{\rm d} r, \quad F_t(t) = 2{\rm \pi}\int_0^{R(Z(t) + l_C,t)} p(Z(t) + l_C,r,t)r\,{\rm d} r.\end{equation}

Let $\tilde {\sigma }_n$ and $\sigma _n=\tilde {\sigma }_n/\tilde {p}_0$ denote the stress acting normally on the portion of the lateral surface of the stopper yet compressed by the bottleneck, and $\mu$ the Coulomb friction coefficient. Hence, $\mu \sigma _n$ is the tangentially acting frictional stress and $(\sin {\theta _B}+\mu \cos {\theta _B})\sigma _n$ the $z$-component of the net stress exerted by the bottleneck on the lateral surface of the stopper yet inside the bottle (see figure 1a). The presumptions (ii) and (iii) in § 2.1 allow for approximating both the bottleneck and the stopper by very long cylinders and the static constitutive law proposed in the premises (iv) by taking $\tilde {\sigma }$ as a monotonically increasing function of the local relative compression, $\varepsilon (z,t)$. With $H(\cdot )$ denoting the Heaviside step function

(2.18)\begin{gather} F_B(t) = 2{\rm \pi} H({-}t)\int^0_{Z(t)} \left(-\dfrac{{\rm d} r_B}{{\rm d} z} + \mu\right) \sigma_n\left(1-\dfrac{r_B(z)}{r_C (z-Z(t))}\right) r_B(z)\,{\rm d} z, \end{gather}

(2.19)\begin{gather}F_{ls}(t) = 2{\rm \pi} H(t)\int^{Z(t) + l_C}_{Z(t)} \dfrac{\partial R(z,t)}{\partial z} p (z,R(z,t),t)R(z,t)\,{\rm d} z. \end{gather}

Relaxing the assumption of a constant negative slope ${\rm d} r_B/{\rm d} z$ of the bottleneck in (2.18) signifies the straightforward advancement towards a more general approach. The movement of the stopper decreases its portion inside the bottle and thus $F_B$, to remain zero for $t\geq 0$. We also anticipate that $F_{ls}=0$ for $t<0$: then, inside the bottle, the stopper slides along the solid wall and, outside, contributions to $F_{ls}$ vanish as its cylindrical shape, when it has just passed the opening, assumed in (v) in § 2.1, gives $\partial R/\partial z=0$.

Closing the relationship $\sigma _n(\varepsilon )$ in supplement A completes the statement of the problem. Supplement B presents the evaluations of (2.18) and (2.19) with the localised variation of $R$ examined in supplement C.

3.1 Numerical scheme: overview

Let us first recall briefly the layers of Clawpack, where the more methodically interested reader is referred to the online software manual by Clawpack Development Team (2022), the outline given by Reference WagnerWagner (2021) and the references in § 1.

Godunov's finite-volume discretisation and Roe's linearisation of a Riemann problem of most general type, posed by some hyperbolic system of transport equations expressed in (physically admissible) conservative form, here specified by the Euler equations (2.4)–(2.6), lie at the core of Clawpack. The stable spatial/temporal resolution of weak solutions, i.e. discontinuities/shocks, requires the numerically stable evaluation of the local waves and corresponding speeds. This is achieved by constructing the global solution from the waves resolved in the directions normal to the interfaces of adjacent finite-volume cells (dimensional splitting) and a positively conservative formulation of the actual Riemann solver (Reference EinfeldtEinfeldt, 1988; Reference Einfeldt, Munz, Roe and SjögreenEinfeldt, Munz, Roe, & Sjögreen, 1991; Reference Harten, Lax and van LeeHarten, Lax, & van Lee, 1983). In practice, the latter requirement can vary greatly in complexity, where the implemented Harten–Lax-van-Leer–Einfeldt (HLLE) scheme by this group of authors provides probably the most prominent one. In each time step, the dependent variables are updated successively via the so obtained scheme and a single call of a second-order Runge–Kutta method. The latter is potentially required to accomplish the time integration of any source terms in the set of equations. Here, their appearance as the right sides of (2.4)–(2.6) resorts to conservation of mass in the radial direction (and thus merely the use of polar coordinates).

3.2 Boundary and coupling conditions in stable scheme

In the following, we give a synopsis of the numerical implementation, adopting two formal ingredients. Supplement D puts forward specific technical details. Further information will be provided upon request.

At first, we introduce the solution vector $\boldsymbol {Q}=(\rho,\rho u,\rho w,e_t)$ and let subscripts indicate the type or the index of the evaluated cell and superscripts the index of the time step considered. In order to start the simulation, every cell must possess initial values, contained in $\boldsymbol {Q}^0=\boldsymbol {Q}$ at $t=t_0$. Hence, $\boldsymbol {Q}^0$ equals either $\boldsymbol {Q}_{{in}}^0=(\tilde {\rho }_B/\tilde {\rho }_0, 0, 0, \tilde {p}_B/\tilde {p}_0)$ or $\boldsymbol {Q}_{{out}}^0=(1, 0, 0, 1)$, where the subscripts ‘in’ and ‘out’ refer to the pressurised gas in the bottle and the ambient one, respectively. The cells occupied by the bottle and stopper are filled with non-numeric values so that they are not updated during the subsequent time steps (black regions in figure 1). Secondly, we take into account the method of dimensional splitting in the subsequent outline. It drastically alleviates the discretisation of the kinematic coupling conditions (2.15) as they can be treated separately for the radial and the axial directions, i.e. as if for independent families of one-dimensional waves.

In Clawpack, boundary conditions are realised by including ghost cells near the edges of the domain and at the interfaces between solid and gaseous phases (for a graphical depiction see supplement D.1). The kinematic constraints (2.12a,b) and (2.15) are treated as reflective walls, and extrapolating boundary conditions simulate an undisturbed flow out of the spatial domain (outflow conditions). The extrapolation condition assigns the first and last internal value to all ghost cells at the left and right or base and top edges, respectively, $\boldsymbol {Q}_j=\boldsymbol {Q}_{N_G+1}$, $\boldsymbol {Q}_{N_l+j}=\boldsymbol {Q}_{N_l}$ ($\,j=1,\ldots,N_G$). Here, $N_G$ denotes the number of ghost cells and $N_l$ the last internal cell index. With regard to (2.15), ghost cells inside the stopper result from mirroring fluid cells at its surface. Velocities are assigned to the ghost cells such that the average value of their and the original flow components perpendicular to the surface match that of its local normal speed ($\dot {Z}$ or $\dot {R}$). All other quantities are simply mirrored along the interface. Its index is set to $m-1/2$ consistently, where $m$ indicates the cell centre: $\rho _{m-1+j}=\rho _{m-j}$, $p_{m-1+j}=p_{m-j}$, $\boldsymbol {v}_{m-1+j}=2\boldsymbol {v}_{C,m-1/2}-\boldsymbol {v}_{m-j}$ ($\,j=1,\ldots,N_G$). While this allocation describes the filling of ghost cells right of a reflective wall, these equations can simply be rearranged to apply to a wall to the left. As a particular novelty, $p$ is mirrored instead of $e_t$, which, as part of the solution, is the intuitive and traditionally used choice (see the aforementioned literature on Clawpack). However, if the method is performed in this classical manner, unacceptable negative values of the pressure inevitably arise inside the HLLE algorithm when the value of the kinetic contribution in (2.2) exceeds a certain threshold.

The maximum value of the variable time step must be smaller than that given by the Courant–Friedrichs–Lewy (CFL) condition (Reference Courant, Friedrichs and LewyCourant, Friedrichs, & Lewy, 1928): $u\Delta t/\Delta r+w\Delta t/\Delta z\leq C_{max}$. For explicit time stepping methods, $C_{max} =1$. A grid with quadratic cells is chosen ($\Delta r=\Delta z$) to minimise distortions. A variety of different cell distances and time steps is used in order to assess the consistency of the discretisation and the temporal convergence of the scheme (§ 4).

The (continuous) motion of the stopper is resolved for discrete but variable time steps $\Delta t_n=t_{n+1}-t_n$ using Taylor approximations. Setting $\tau =\Delta t_n/\Delta t_{n-1}$, we approximate the base position of the stopper $Z(t_n)$ and the corresponding velocity $\dot {Z}(t_n)$ with $Z^n$, $\dot {Z}^n$ as

(3.1) \begin{align} Z^{n+1} &= -\tau\,Z^{n-1}+(1+\tau)\,(Z^n+\ddot{Z}^n\,\Delta t_n\,\Delta t_{n-1}/2)+\mathcal{O}(\Delta t^3)\,, \end{align}

(3.2) \begin{align} \dot{Z}^{n+1} &= \tau^2\,\dot{Z}^{n-1}+(1-\tau^2)\,\dot{Z}^n+\Delta t_n (1+\tau)\,\ddot{Z}^n+\mathcal{O}(\Delta t^3)\,. \end{align}