Computational Fluid Dynamic simulation of the supersonic CO2 expansion during champagne cork popping

Behind the iconic “pop!" accompanying the uncorking of a champagne bottle hides a gas ow of surprising complexity. Its modeling is made delicate by its supersonic nature, its interaction with the cork stopper, the eminently unsteady character of the ow escaping from the bottle, and the continuous change of the geometry of the computational ow domain due to the displacement of the cork. Our numerical simulations reveal the formation, evolution and dissipation of shock wave patterns during the rst millisecond following cork popping. A rst crown-shaped shock wave pattern develops radially, then followed by the formation of a detached shock wave, or bow shock, induced by the presence of the cork in the axial path of the supersonic gas. The numerical simulations correctly reproduce the position of the bow shock previously observed through high-speed imaging.


Introduction
Nowadays, after more than three centuries and undergoing continuous re ning, champagne has undoubtedly become the most renowned French sparkling wine, praised world-wide for the neness of its carbon dioxide (CO 2 ) bubbles. 1 Indeed, during champagne making, bottles reach a high pressure because gas-phase CO 2 forms together with ethanol during a key step in-bottle fermentation process promoted by adding yeasts and sugar in the bottles hermetically sealed with a cork stopper. 1 This is basic application of Henry's law which states that the concentration of dissolved CO 2 in the liquid phase is proportional to the partial pressure of gas-phase CO 2 in the bottleneck. Moreover, because the solubility of CO 2 in the liquid phase is strongly temperature-dependent (the lower the temperature of the wine, the higher the gas solubility), the partial pressure of gas-phase CO 2 in the bottleneck is therefore also strongly temperaturedependent (the lower the temperature of the wine, the lower the partial pressure of gas phase CO 2 ).
Otherwise, because the driving force behind the cork popping process is the force exerted by gases under pressure in the bottleneck on the base of the cork stopper, the cork popping velocity is therefore de nitely under the in uence of the champagne temperature, as already shown in a previous article. 2 During the cork popping process, concomitantly with the cork stopper being expelled from the bottleneck under the action of pressure, a plume made of a gas mixture mainly composed of CO 2 (with traces of water vapor) freely expands out of the bottleneck through ambient air and experiences adiabatic cooling. [3][4][5] High-speed video imaging was used recently to visualize champagne cork popping, and especially the condensation processes following the adiabatic expansion of the gas mixture found in the bottlenecks of transparent champagne bottles stored at different temperatures. 5 The bottles' temperature (and therefore their inner pressure) was found to be a key parameter concerning the condensation processes that can occur above, and inside the bottlenecks. Most interestingly, for the bottles stored at 20°C (under a pressure close to 7.5 bar), the characteristic grey-white cloud of fog classically observed above the bottlenecks of bottles stored at lower temperatures was replaced by a more evanescent plume, surprisingly blue, starting from inside the bottleneck. 5 Depending on the strongly bottle temperature dependent saturation ratio experienced by gas-phase CO 2 after adiabatic expansion, it was emphasized that blue haze is the signature of a partial and transient heterogeneous freezing of gas phase CO 2 on ice water clusters homogeneously nucleated in the bottlenecks. 5 Even more recently, a comparison between the condensation phenomena accompanying cork popping from bottles stored at 20 and 30°C was made. 6 The initial bottleneck-to-ambient-pressure ratio (7.5 for the bottles stored at 20°C, against about 10 for those stored at 30°C) much exceeded the critical ratio needed for the CO 2 /H 2 O gas mixture to reach Mach 1, thus forming under-expanded supersonic jets expelled from the throat of the bottlenecks.
Moreover, during the very rst millisecond following the cork popping process, evanescent normal shock waves were unveiled in the plumes, until the reservoir-to-ambient-pressure ratio goes below a critical ratio needed for their formation, which were de facto seen vanishing. 6 Actually, the expansion of a high-pressure gas suddenly opening into a low-pressure environment gives rise to a so-called free jet, as opposed to a jet con ned by the walls of a duct. The free jet acquires supersonic speed when the speed of sound is reached at the outlet that separates the high and low pressure areas. This condition is met when the ratio of the upstream and downstream pressures is greater than approximately 2 (the precise value depends on the ratio of the speci c heats of the gas). The free jet can then reach very high Mach numbers, which corresponds to high speeds but also to temperatures which can be extremely low, close to absolute zero, 7 under the law of conservation of energy. This speci city of supersonic ows is used in particular to carry out studies of chemical reactivity 8 or spectroscopy 9 at very low temperature in the gas phase. The interaction of the supersonic free jet with the downstream static gas leads to the formation of a complex network of shock waves and boundary layers, also called diamond shocks by Ernst Mach who rst described them using shadowgraphy. 10,11 Note that at high pressure ratios, a normal shock wave forms, perpendicular to the ow, through which the gas becomes subsonic again. For a circular ori ce, this normal shock wave is called a Mach disk. Such a shock wave occurs at the outlet of Laval's nozzles which propel rockets, missiles or scramjets, but can just as easily form during the eruption of a volcano, [12][13][14] or during the formation of geysers observed on the surface of the Europa satellite. 15 The Mach disk is formed all the further from the gas outlet the greater the upstream and downstream pressure ratio. 16 The position of the Mach disk remains stationary as long as the upstream and downstream pressures do not change over time, which is not the case if the volume of the high pressure gas is nite, as is the case for example for experiments carried out in a shock tube 17 or for gas escaping from a bottle of champagne. To simulate this type of expansion, it is then necessary both to take into account the nite nature of the high pressure reservoir 18 and consider an unsteady ow. An additional di culty in the case of the bottle of champagne is linked to the presence of an obstacle on which the ow is impacting. Supersonic ow hitting a surface leads to the formation of a detached shock wave 19 through which the gas will adapt to the low speed and high pressure conditions imposed by the presence of the obstacle. This situation is very similar to the cold spray process device which consists in ejecting particles at supersonic speeds which collide with a substrate to be treated. 20 In this case, two shock wave systems interfere: a free jet type structure and a detached shock wave type structure. The situation is indeed quite similar during the champagne cork popping process. The gas mixture expelled from the bottleneck must adapt to the ambient pressure through a network of shock waves but also face the presence of the cork that is in its path. Nevertehless, a signi cant difference from the study by Pattison et al. 20 is that the distance between the gas outlet (the bottleneck) and the obstacle (the ying cork) varies over time.
The present study addresses the aforementioned issues by considering the champagne cork popping process as a supersonic free jet expelled from a nite reservoir and impacting a moving obstacle. Based on Computational Fluid Dynamics (CFD) calculations, the shock wave structure experienced by the gas ow expelled at supersonic speeds from the champagne bottleneck after ejection of the cork stopper was followed with time. The calculation performance and the stability of the two-dimensional axisymmetric method used were evaluated by comparing the numerical results of the transient gas expansion with the normal shock wave position previously observed through high-speed video imaging. 5,6 Results Shock wave pattern. According to gas dynamics and as con rmed by our simulations shown in Fig. 1, the gas escaping from the bottle reaches a supersonic regime. Before describing into details the ow eld properties revealed by our numerical simulations, let's brie y recall some basics of compressible ow 21 .
Two conditions have to be met by the gas ow to break the speed of sound. The rst one relates to the existence of a throat, i.e. a spatial zone delimiting a converging duct area from a diverging one. Due to the mass conservation, a subsonic ow needs to move in a converging duct to be accelerated.
Conversely, the decrease of density of a supersonic ow is so rapid that it cannot be compensated anymore by the velocity increase to maintain its mass ow (ṁ = ρAV). Both velocity (V) and area (A) must increase to compensate for the density (ρ) decrease of a supersonic ow. It follows that the speed of sound can only be reached if the gas ow passes through a throat. The shoulder of the bottle of champagne behaves as a converging duct whereas the diverging one is only virtual, it corresponds to the large atmospheric volume in which the high pressure gas is "freely" discharging once the cork popped. The second condition to reach a supersonic regime is that the ow must reach the speed of sound at the throat, which then behaves like a sonic throat characterized by a Mach number equal to 1, this condition is called chocked ow. In case this condition is not satis ed, the gas will decelerate downstream the throat. A chocked ow is obtained provided that the gas reaches the critical pressure p * at the throat which is given by the isentropic chocking pressure ratio where p 0 is the pressure prevailing in the bottle and where γ denotes the speci c heat ratio of the gas. Thus, a supersonic ow is developing provided that p * /p 00 .5, or equivalently, that the nozzle pressure ratio (NPR) p 0 /p * 2 .
As the initial pressure of a bottle of champagne far exceeds two atmospheres, it is certain that the gas it contains will reach a supersonic regime when discharging in the atmosphere. This is well con rmed by the complex shock wave pattern (shock diamonds) very speci c to supersonic ows and rst described by Mach, 10 which is clearly visible in our CFD simulations ( Fig. 1 (500 µs) and Fig. 2 (518 µs)).

( )
The CFD simulations reveal that the cork popping is followed by successive owing phases, the reader must refer to the simulations rows shown in Fig. 1 corresponding to 500, 917 and 1167 µs, respectively, as well as to the evolution of the Mach number of the gas versus time presented in Fig. 2 from 166 to 1577 µs.
During the early phase closely following the cork ejection (for example for a owing time of 500 µs as shown in Fig. 1; see also Fig. 2 (166, 249, 415 and 518 µs)), the cork is so close to the throat of the bottleneck that the pressure of the gas located between the cork and the bottle is still several bars, much higher than the critical pressure. The ow comprised in the small volume between the bottleneck and the cork stopper is not chocked along the bottle axis and remains subsonic. The gas is accelerated until it reaches a maximum speed of about 200 m.s −1 at the bottleneck. This speed corresponds to a Mach number of 0.78 (Fig. 3) which indicates that the speed of sound has not yet been reached at the bottleneck. This results in deceleration of the gas after the bottleneck to zero speed at the cork. Interestingly, the temperature of the gas at the bottleneck drops to -43°C before rising to -10°C at the cork, thus explaining the formation of water ice observed on the internal face of the cork during this initial phase of ejection. 6 The high pressure gas in the region under the cork has no alternative but to escape laterally into atmospheric air giving rise to a nice crown shaped supersonic expansion characterized by a typical succession of shock cells along which the expanding gas gradually equilibrates its static pressure to ambient pressure. Although the lateral diamond shocks were not revealed by the high speed camera, 5,6 their presence seems to be con rmed by another unpublished numerical simulation. 22 From Fig. 2, it can be seen that the Mach number of the supersonic jet crown reaches about 2.9 in the very early stage of the lateral gas expansion (< 249 µs).
A second phase begins once the cork is far enough from the bottleneck so that the pressure downstream is lower than the critical pressure. In that case, the bottleneck can behave as a sonic throat and the gas ow is accelerated downstream to reach supersonic speed up to Mach 1.6-1.8 (see second raw of Fig. 1, relative to 917 µs; see also panels 747, 913 and 996 µs of Fig. 2). The enthalpy of the gas is converted into kinetic energy leading to a drop in temperature down to -95°C then triggering the condensation of CO 2 , as observed in our previous work. 6 The uid molecules then impinges on the cork face and rebound leading to the formation of a bow shock which can be seen as the coalescence of re ected signals by the object. The shock is a compression front across which the ow properties jump. According to the second law of thermodynamics, a ow behind a normal shock wave is necessarily subsonic 21 . The distance between the bow shock and the object, i.e. the thickness of the shock layer, depends on the Mach number and cannot be determined using a simple analytical expression, it needs to be determined by numerical calculation methods as the one used in the present work.
As con rmed by our simulations, the ow behind the shock wave is subsonic, its temperature is raised to ∼ -60°C and the gas is recompressed to about 2.2 bar so that it further expands laterally up to recover a supersonic regime responsible for the formation of lateral shock cells. The simulations clearly show that the radial ejection of the gas during the initial phase (phase 1) is followed by a purely axial ejection once the normal shock wave is formed at the bottle outlet (phase 2). Gas escapes laterally only after passing through the shock wave and impacting on the plug.
During the third phase the pressure of the gas in the bottle has become too low to maintain a NPR > 2, the chocked conditions are not anymore satis ed, neither at the bottleneck nor at the edge of cork. The normal shock wave and the crown shaped lateral shock cells can no longer form (see last raw in Fig. 1 corresponding to a time ow of 1167 µs; see also Fig. 2 for times greater than 996 µs). The ow is now subsonic and a rarefaction wave propagates in the bottleneck.
Position of the detached shock. Figure 4  , where D is the diameter of bottle outlet. The Mach disk gets closer as gas escapes and cylinder pressure drops. However, it remains much further from the opening of the bottle than the bow shock, which shows that it is the cork that controls the formation of the shock waves.
The normal shock wave could be observed experimentally (Fig. 5) due to the light scattered by the solid CO 2 particles which formation is triggered by the very low temperature reached by the ow during its adiabatic cooling. 5,6 The calculated distance of the shock wave from the cork is generally in good agreement with the observations, with nevertheless a better agreement for longer owing times (Fig. 6). These small disagreements could be explained by the presence of condensation of water and CO 2 molecules which is not taken into account in our CFD simulations although responsible for an important release of heat of condensation in the ow. Not surprisingly, the shock detachment distance reduces when the Mach number increases, i.e. with increasing bottle pressure. Although di cult to model and predict, it is a well-known effect observed in many studies of gas dynamics. 19,23 Additional simulations were performed for lower initial pressure conditions of 3, 5 and 6 bar (Fig. 7). They do con rm observations: the normal shock wave is as close to the cork stopper that upstream pressure is strong.

Discussion
Our CFD simulations reveal the existence of three distinct owing phases following the cork popping of a bottle of champagne. In the early phase, between 0 and 600 µs approximately (time 0 corresponding to the precise moment when the cork stopper separates from the bottleneck), the gas which leaves the bottle is axially blocked by the cork. The situation is then very different from a gas which would escape freely from a high-pressure reservoir into the atmosphere. In the case of an unobstructed free stream ( ) downstream, the gas outlet would behave like a sonic throat, causing the gas to accelerate beyond the speed of sound. In the presence of the stopper and during this initial phase, the gas which escapes from the bottle does not have the possibility of expanding su ciently along the axis of the cylinder, so that the Mach number at the bottleneck remains less than 1. The bottleneck therefore does not behave like a sonic throat: the gas reaches its maximum (subsonic) speed at the bottleneck after its acceleration under the effect of the progressive narrowing of the outlet section of the bottle, then decelerates downstream of the bottleneck under the effect of the sudden increase in the outlet section. In this rst phase of the process, the gas located between the cork stopper and the bottle outlet is characterized by a pressure of several bars and it will obviously seek to balance its static pressure with atmospheric pressure by escaping laterally through the annular space which grows as the cork is ejected. Numerical simulations show that this space acts actually as an annular sonic throat. The gas escapes radially at supersonic speed, balancing its pressure through a succession of right and oblique shock waves forming the usual structure of diamond shock waves clearly visible in the simulations (see Figs. 1 and 2), note that the axial symmetry of the bottle leads to a crown-shaped supersonic expansion.
In a second phase, typically between 600 and 1000 µs, the gas pressure at the bottle outlet has become low enough for the critical condition (Mach 1) to be reached at the bottleneck (see e.g. Fig. 3). Then gas escapes from the bottle at supersonic speed, this time directed along the axis of the cylinder. The gas then impacts the cork stopper and becomes subsonic again passing through a detached shock (bow shock) induced by the presence of the cork stopper. The position of this bow shock wave agrees well with the observations. It must be underlined that this shock wave corresponds to a bow shock rather than the characteristic Mach disk of non-adapted free jets that would form further away in the absence of a cork stopper.
Finally, in the last phase, typically beyond 1000 µs, the pressure in the bottle is too low to maintain the critical pressure at the bottleneck. The gas leaves the bottle at subsonic speeds and then decelerates downstream. The simulations clearly show the propagation of a rarefaction wave inside the bottle during this last phase.

Methods
Computational domain and initial conditions. The axisymmetric geometry approximation is chosen because it is well suited to the axial symmetry of the ow eld corresponding to the gas ejection during champagne cork popping from a bottle held vertically. One of the calculation domains (Fig. 8) concerns the internal gaseous volume of approximately 35 cm 3 , under pressure, found within the bottleneck of the champagne sealed bottle (i.e., the gas phase between the liquid surface and the cork stopper). This gas mixture consists in gas-phase CO 2 , with traces of water vapor, nitrogen, ethanol vapor and a myriad of other volatile compounds. 1 Nevertheless, the partial pressure of gas-phase CO 2 is largely dominant compared with that of all other gaseous compounds. 1 For the sake of simplicity, the numerical simulations were therefore carried out by considering that the gas phase under pressure in the bottleneck of the corked bottle consisted in pure CO 2 . Moreover, the initial pressure and temperature were respectively set at 7.5 bar and 20°C, in order to comply with one of the experimental conditions imposed in previous studies. 5,6 The calculations are based on the exact geometric pro le of a bottle of champagne. In the headspace, the section of the bottleneck gradually reduces until it reaches an internal diameter close to 1.8 cm at the outlet of the bottle. The other part of the calculation eld concerns the external environment formed by the ambient air in which the CO 2 will expand freely.
The movement of the champagne cork is controlled with a user-de ned function (UDF) imported from the ANSYS Fluent software. This function xes the ejection speed of the cork stopper at 18.6 m.s −1 established from experimental observations made under the same owing conditions 5 (Fig. 9). Note that for the sake of simplicity, in our simulations the elasticity of the cork has not been taken into account, which explains why in our simulations the base of the cork does not adopt the typical mushroom shape of a champagne popped cork (see inserted picture in Fig. 4) but remains cylindrical once removed from the bottle (Figs. 1 and 2).
Calculation method. Our uid dynamics calculation study is based on the commercial ANSYS Fluent software. More speci cally, in the present study, Navier-Stokes axisymmetric equations and species transport equations are used, and all ow variables are used in the mixing forms. The k-epsilon turbulence model ("Realizable" option) is used, and the Near-Wall Treatment ("Enhanced Wall Treatment" option) is selected. The "pressure-velocity coupling" algorithm is chosen for the calculation process. The "PRESTO" model is used for the discretization of the pressure and to obtain a good precision of the pressure gradients near the walls. The "Second Order Upwind" model is selected for the discretization of momentum and energy equations. Considering the cork stopper velocity and the Courant-Friedrichs-Lewy (CFL) condition, the computational time interval is set to 10 −6 s, and the maximum number of iterations for each time step is set to 200. The external limits of the domain of integration were maintained at ambient temperature and pressure.
Dynamic mesh. The calculation domain is modi ed as the cork stopper is ejected. It is therefore necessary to use a dynamic mesh to simulate the gas ow throughout the ejection of the cork stopper.
Smoothing and remeshing methods are used and properly parameterized to maintain a high mesh quality and avoid any degeneration of the grid during its updating and deformation throughout the calculation process. The diffusion-based smoothing method is used for the mesh deformation with a limit distance and a diffusion parameter set to 0.1. The local remeshing method is used when the mesh violates the skewness or size criteria 24 . The minimum and maximum length scales parameters are set to 10 −6 m and 0.002 m, respectively, and the maximum cell skewness is xed to 0.7. Simulation of pure CO 2 expanding freely into ambient air (20°C, 1.013 bar) following the cork popping of a champagne bottle initially kept at 20 °C with an inner pressure of 7.5 bar. The rst, second and third panel rows relate respectively to owing times of 500, 917 and 1167 ms measured after the cork popping. The rst, second and third panel columns are relative to the ow velocity, static pressure and temperature, respectively.  Mach number axial pro le of a gas escaping from a bottle initially at 20°C and 7.5 bar for owing times of 500, 917 and 1167 ms. At 500 µs (phase 1) the exit pressure is too high so that the gas ow is not chocked along the bottle axis, it decelerates downstream the bottleneck. At 917 µs (phase 2), the ow is supersonic and a detached shock is forming in front of the cork. At 1167 µs (phase 3), the pressure in the bottle is too low and nozzle pressure ratio (NPR) is not high enough for maintaining a chocked ow on the bottle axis.   Calculated and observed positions of the normal shock wave. Here, the bottle was initially kept at 20°C and 7.5 bar, both simulation and picture corresponds to a owing time of 917 µs after the cork popping.
The origin of the x-axis corresponds to the surface of the liquid in the bottle.

Figure 6
Temporal evolution of the calculated and observed positions of the detached shock wave forming in front of the cork.