2.1. The RHD for the plasma

A plasma in the presence of radiation photons can be modelled with a set of RHD equations (Castor Reference Castor2004) in the form

(2.1) \begin{equation} \begin{cases} \partial _{t} \rho _{g} +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{g}\boldsymbol{u}_{g} ) =0 , \\[3pt] \partial _{t} ( \rho _{g}\boldsymbol{u}_{g} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{g} \boldsymbol{u}_{g}\boldsymbol{u}_{g} )+\boldsymbol{\nabla }(p_{g}+p_{r} ) =0 , \\[3pt] \partial _{t} ( \rho _{g} e_{g} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{g} e_{g} \boldsymbol{u}_{g} ) +p_{g}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{g}=\boldsymbol{\nabla }\boldsymbol{\cdot }( \lambda _{g}\boldsymbol{\nabla }T_{g} ) +\omega _{r,g}\big( T_{r}^{4}-T_{g}^{4} \big) , \\[3pt] \partial _{t} E_{r} +\boldsymbol{\nabla }\boldsymbol{\cdot }( E_{r} \boldsymbol{u} ) +p_{r}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=\boldsymbol{\nabla }\boldsymbol{\cdot }\big( \chi _{r,g}\boldsymbol{\nabla }T_{r}^{4} \big) + \omega _{r,g}\big( T_{g}^{4}-T_{r}^{4}\big ), \end{cases} \end{equation}

where the subscripts $g$ and $r$ denote plasma (i.e. ionised gas) and radiation, respectively, $t$ is the time, and the plasma quantities include the mass density $\rho _{g}$ , the macroscopic velocity $\boldsymbol{u}_{g}$ , the pressure $p_{g}$ , the temperature $T_{g}$ , the specific internal energy $e_{g}$ , and the thermal conduction coefficient $\lambda _{g}$ . The two components of the plasma, i.e. electrons and ions, are assumed to have identical macroscopic fluid velocity and temperature, and act as a single fluid in this model. Radiation has been modelled under a diffusion approximation, which is valid in an optically thick regime (Castor Reference Castor2004) where the mean free paths of radiative photons are much shorter than the characteristic length scale of hydrodynamics. The radiation photons act as an ideal gas in an optically thick regime (Castor Reference Castor2004), and provide a radiation pressure $p_r$ . They are assumed to follow the Planckian (black-body) distribution with radiation temperature $T_r$ that is allowed to differ from the local plasma temperature $T_g$ . The radiation energy $E_r = a T_{r}^{4}$ (where $a=7.5657\times 10^{-16}\ \mathrm{J}\ \mathrm{m}^{-3}\ \mathrm{K}^{-4}$ is the radiation constant) given by the Planckian distribution is linked with $p_r$ as $p_{r} = E_{r}/3$ due to the Eddington approximation (Castor Reference Castor2004). Once $T_r$ differs from $T_g$ , the radiation photons exchange energy with the plasma at an energy exchange rate $\omega _{r,g}$ . Under the diffusion approximation, the energy transported by radiation photons can be effectively separated into two components: one part acts as an energy diffusion process, characterised by the coefficient $\chi _{r,g}$ , while the remaining portion contributes to energy convection alongside the plasma flow, as demonstrated in the last equation of (2.1). The simplified radiation model adopted here, which relies on the Planckian distribution approximation, is generally sufficient to capture the essential hydrodynamic behaviour. This approach has been successfully applied in LPP-EUV simulations (Yuan et al. Reference Yuan, Ma, Wang, Chen, Cui, Zi, Yang, Zhang and Leng2022), and proves adequate for the pre-pulse simulations performed in our study.

The last two equations in (2.1) are the energy equations for the plasma and the radiation, respectively. Alternatively, the combination of these two equations yields the equation for the total energy of the system:

(2.2) \begin{align} \partial _{t} \big [ \rho _{g}\big ( e_{g}+| \boldsymbol{u}_{g} |^{2}/2 \big )+E_{r}\big ] &+\boldsymbol{\nabla }\boldsymbol{\cdot }\big [ \rho _{g}\big ( e_{g}+| \boldsymbol{u}_{g} |^{2}/2 \big ) \boldsymbol{u}_g +E_{r} \boldsymbol{u}_{g}\big ] + \boldsymbol{\nabla }\boldsymbol{\cdot }\big [ ( p_{g}+p_{r} )\boldsymbol{u}_g \big ] \nonumber \\ &=\boldsymbol{\nabla }\boldsymbol{\cdot }( \lambda _{g}\boldsymbol{\nabla }T_{g} ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\big( \chi _{r,g}\boldsymbol{\nabla }T_{r}^{4} \big). \end{align}

Table 1.

The SG parameters for two liquids considered in the present work.

Table 1. Long description

A table comparing properties of water and liquid tin. The table has three columns and two rows. The columns are labeled Fluid, γ, P∞,f (Pa), and cv,f (J kg-1 K-1). The row labels are Water and Liquid tin. Row 1: Fluid, Water; γ, 4.4; P∞,f (Pa), 6 x 10^8; cv,f (J kg-1 K-1), 1816. Row 2: Fluid, Liquid tin; γ, 30; P∞,f (Pa), 1.421 x 10^9; cv,f (J kg-1 K-1), 210.

2.2. Euler equations for the liquid

The liquid is assumed to be inviscid and weakly compressible, and its dynamics is governed by the Euler equations

(2.3) \begin{equation} \begin{cases} \partial _t \rho _{l} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _{l} \boldsymbol{u}_{l}) = 0, \\[3pt] \partial _t (\rho _{l} \boldsymbol{u}_{l}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _{l} \boldsymbol{u}_{l}\boldsymbol{u}_{l}) + \boldsymbol{\nabla }p_{l} = 0, \\[3pt] \partial _t \big[\rho _{l} \big(e_{l}+ | \boldsymbol{u}_{l} |^{2}/2\big) \big] + \boldsymbol{\nabla }\boldsymbol{\cdot }\big[\rho _{l} \big(e_{l}+ | \boldsymbol{u}_{l} |^{2}/2 \big) \boldsymbol{u}_{l}\big] + \boldsymbol{\nabla }\boldsymbol{\cdot }(p_{l}\boldsymbol{u}_{l}) = 0, \end{cases} \end{equation}

where subscript $l$ denotes liquid, and $\rho _{l}$ , $\boldsymbol{u}_{l}$ , $p_{l}$ and $e_{l}$ are density, velocity, pressure and specific internal energy of the liquid, respectively. In order to be consistent with the RHD equations (2.1), the energy equation in (2.3) is reformulated as the transport equation for the internal energy:

(2.4) \begin{equation} \partial _t (\rho _{l} e_{l}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _{l} e_{l} \boldsymbol{u}_{l}) + p_{l}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{l} = 0. \end{equation}

The liquid is assumed to be very opaque to radiation, and the self-emission of the liquid is negligibly weak due to its low temperature. Although placed in an radiative environment, the radiation energy of and the radiation transport inside the liquid are neglected, and thus are absent in (2.3). So is heat conduction in the liquid, which is also negligible compared to that in the plasma.

2.3. The EOS

The plasma and the liquid have different material properties and thus obey different forms of EOS. We employ the stiffened gas (SG) model (Shyue Reference Shyue1998), which is able to describe the EOS of different materials in a unified form by tuning parameters, for both fluids.

The EOS of the liquid reads

(2.5) \begin{equation} \begin{cases} p_{l}+\gamma _{l}{P}_{\infty ,l} = (\gamma _{l} -1)\rho _{l} e_{l},\\[3pt] \rho _{l} e_{l} = \rho _{l} c_{v,l} T_{l}+{P}_{\infty ,l}, \end{cases} \end{equation}

where the specific heat ratio $\gamma _{l}$ , the reference pressure ${P}_{\infty ,l}$ , and the specified heat capacity at constant volume $c_{v,l}$ (Le Métayer et al. Reference Le Métayer, Massoni and Saurel2004) are three constant parameters that need to be assigned. The values of $\gamma _{l}$ and ${P}_{\infty ,l}$ are usually chosen such that the sound speed inside the liquid,

(2.6) \begin{equation} C_{l}=\sqrt {\gamma _{l}\left (p_{l}+ {P}_{\infty ,l}\right )/\rho _{l}}, \end{equation}

is close to reality (Shyue Reference Shyue1998). Typical choices of $\gamma _{l}$ , ${P}_{\infty ,l}$ and $c_{v,l}$ for two liquids in the present work are listed in table 1, to recover a reference state under the standard atmospheric pressure: $T_l = 300\ \mathrm{K}$ , $\rho _{l} = 1000\ \mathrm{kg\ m^{-3}}$ , $C_{l} = 1450\ \mathrm{m\,s}^{-1}$ for water; $T_l = 593\ \mathrm{K}$ , $\rho _{l} = 6980\ \mathrm{kg\ m^{-3}}$ , $C_{l} = 2471\ \mathrm{m\,s}^{-1}$ for liquid tin.

A weakly coupled plasma behaves like an ideal gas (Krall & Trivelpiece Reference Krall and Trivelpiece1973), corresponding to a particular SG with the reference pressure ${P}_{\infty ,g}$ set to zero. Its EOS reads

(2.7) \begin{equation} \begin{cases} p_{g} = (\gamma _{g} -1)\rho _{g} e_{g},\\[3pt] \rho _{g} e_{g} = \rho _{g} c_{v,g} T_{g}, \end{cases} \end{equation}

where the specific heat ratio $\gamma _{g}$ is set to $5/3$ . Here, $c_{v,g}$ is the specified capacity at constant volume considering the average ionisation degree ( $Z_{g}$ ) as

(2.8) \begin{equation} c_{v,g} = \frac {R (1+Z_{g}) }{(\gamma _{g}-1)M_{g}}. \end{equation}

Here, $R=8.314\ \mathrm{J}\,(\mathrm{mol}\boldsymbol{\, }\mathrm{K} )^{-1}$ is the universal gas constant, and $M_{g}$ is the molar mass of the element ( $118.71\ \mathrm{g}\,\mathrm{mol}^{-1}$ for tin). The sound speed in the plasma is written in a form similar to (2.6):

(2.9) \begin{equation} C_{g}=\sqrt {\gamma _{g} p_{g}/\rho _{g}}. \end{equation}

Note that $C_g$ can be an order of magnitude larger than $C_l$ , e.g. a typical tin plasma emitting the interested EUV lights has $T_g$ ranging from $10^{5}$ to $5\times 10^{5}$ K (Nishihara et al. Reference Nishihara2008), yielding $C_g$ from $11\,000$ to $25\,000$ $\mathrm{m\,s}^{-1}$ .

3.1. Jump conditions at interfaces

In this work, we model the interface between plasma and liquid as a sharp contact discontinuity with no phase transition or chemical reaction, as illustrated in figure 1( $a$ ). The two phases are spatially separated and immiscible, with distinctly different material properties, thereby leading to jump conditions at interfaces. Since the interface is essentially a contact discontinuity, the normal velocity and pressure are continuous across the interface in the absence of surface tension, i.e. $\boldsymbol {u}_l \boldsymbol{\cdot }\boldsymbol {n} =\boldsymbol{u}_g \boldsymbol{\cdot }\boldsymbol {n}$ and $p_{l}=p_{g}$ , where $\boldsymbol {n}$ is the unit normal vector to the interface, while the density is discontinuous, $\rho _{l} \ne \rho _{g}$ . In contrast, surface tension induces a pressure jump across the interface. Apart from these, the thermal fluxes in the plasma are assumed to vanish at interfaces, due to the comparatively negligible thermal conduction coefficient of the liquid. Similarly, the radiation energy fluxes are also expected to vanish at the interface because of high opacity of the liquid. In the present study, all these conditions are well handled by the proposed diffuse interface model introduced in § 3.2 together with the mixing rules introduced in §§ 3.3 and 3.4.

3.2. Diffuse interfaces

In diffuse interface models, the sharp interface between two immiscible fluids is replaced by a diffuse interface with finite thickness where the two fluids are mixed, and the jump conditions at the interface are realised by mixing rules in the diffuse interface region. Inspired by the concept of the diffuse interface, specifically the transport five-equation model (Allaire et al. Reference Allaire, Clerc and Kokh2002), we propose a diffuse interface model for liquid–plasma interactions, where the volume fraction of the liquid ( $\alpha _l$ ) is used to represent the interface position, as illustrated in figure 1( $b$ ). We can see that $\alpha _l$ rapidly changes from 1 on the liquid side to 0 on the plasma side within a thin layer, which is referred hereafter to as the diffuse interface region. In the diffuse interface region, the fluid is treated as a homogeneous mixture of the liquid and the plasma.

Similar to the model of Allaire et al. (Reference Allaire, Clerc and Kokh2002), the proposed model assumes that the velocities and pressures of the two phases are in equilibrium (namely mechanical equilibrium) in the diffuse interface region – i.e. the mixture velocities are $\boldsymbol{u} = \boldsymbol{u}_l = \boldsymbol{u}_g$ , and the mixture pressures are $p = p_l = p_g$ – but can be in thermal non-equilibrium ( $T_l \ne T_g$ ). Such mechanical equilibrium at the interface ensures the proper propagation of acoustic waves across material interfaces while allowing the fluids to have different densities, internal energies and EOS. Following the continuous surface force formulation (Garrick, Owkes & Regele Reference Garrick, Owkes and Regele2017), the surface tension is incorporated as a volumetric force in the momentum equation. As a result, we can tentatively write the diffuse interface model for liquid–plasma interactions, which consists of an interface evolution equation, two phasic mass equations, a mixture momentum equation, a mixture energy equation and a radiation energy equation (note that the liquid-phase region is also encompassed by the radiation field, with the radiation energy being set to zero):

(3.1) \begin{equation} \begin{cases} \partial _{t} \alpha _{l}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\alpha _{l}= 0 , \\[3pt] \partial _{t} ( \rho _{l} \alpha _{l} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{l} \alpha _{l}\boldsymbol{u}) =0 ,\\[3pt] \partial _{t} ( \rho _{g} \alpha _{g} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{g} \alpha _{g}\boldsymbol{u}) =0 ,\\[3pt] \partial _{t} ( \rho \boldsymbol{u} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho \boldsymbol{u}\boldsymbol{u} )+\boldsymbol{\nabla }(p+p_{r}) =\sigma \kappa\, \boldsymbol{\nabla }\phi , \\[3pt] \partial _{t} ( \rho e ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho e \boldsymbol{u} ) +p\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}= -\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}+ S_{r}, \\[3pt] \partial _{t} E_{r} +\boldsymbol{\nabla }\boldsymbol{\cdot }( E_{r} \boldsymbol{u} ) +p_{r}\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{F}_{\!r} -S_{r}, \end{cases} \end{equation}

where $\rho$ and $e$ are the density and the specific internal energy of the mixture, respectively, $\sigma$ is the surface tension coefficient, $\kappa = -\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\nabla }\phi /| \boldsymbol{\nabla }\phi |)$ is the interface curvature, and $\phi$ is an order parameter to identify the interface region where the surface tension is applied.

In our diffuse interface model, the absence of an anti-diffusion term implies that the interface thickness (represented by the $\alpha$ contours) may increase over time. To restrict the surface tension force to the vicinity of the droplet interface, we therefore employ $\phi$ instead of $\alpha$ . The parameter $\phi$ is defined as

(3.2) \begin{equation} \phi =\begin{cases} 1 & \text{ if } \alpha _{l}\gt 0.9, \\[3pt] \min \left [ \max \left ( (\rho -\rho _{\textit{min}})/(\rho _{\textit{max}}-\rho _{\textit{min}}), 0 \right ) , 1 \right ] & \text{ if } 0.1\leqslant \alpha _{l}\leqslant 0.9, \\[3pt] 0 & \text{ if } \alpha _{l}\lt 0.1, \end{cases} \end{equation}

where $\rho _{\textit{min}}$ and $\rho _{\textit{max}}$ are set to $1\ \mathrm{g\,cm}^{-3}$ and $4\ \mathrm{g\,cm}^{-3}$ , respectively. As a result, the surface tension force is confined to the interface region bounded by $[\rho _{\textit{min}}, \rho _{\textit{max}}]$ and $0.1\leqslant \alpha _{l}\leqslant 0.9$ .

The energy flux due to thermal conduction $\boldsymbol{q}$ , the energy deposition rate to the mixture from radiation $S_{r}$ , and the energy flux due to radiation diffusion $\boldsymbol{F}_{\!r}$ , are applied to the mixture; their detailed expressions are provided in § 3.4. It is worth noting that the proposed model of (3.1) reduces to the original transport five-equation model for compressible two-phase flows (Allaire et al. Reference Allaire, Clerc and Kokh2002) in the absence of thermal conduction and radiation transport (i.e. $p_{r}$ , $E_{r}$ , $S_{r}$ , $\boldsymbol{q}$ and $\boldsymbol{F}_{\!r}$ are all zero). Equation (3.1) are not closed yet until the EOS of the mixture is supplied.

3.3. The EOS in the diffuse interface region

In the absence of phase transition and chemical reaction, we establish the EOS for the mixture in a similar manner as in Allaire et al. (Reference Allaire, Clerc and Kokh2002). In the diffuse interface region, the conservation of mass and energy straightforwardly leads to $\rho$ and $e$ in the forms

(3.3) \begin{align} \rho = \rho _{l}\alpha _{l}+\rho _{g}\alpha _{g}, \end{align}

(3.4) \begin{align} \rho e = \rho _{l}\alpha _{l}e_{l}+\rho _{g}\alpha _{g}e_{g}, \end{align}

respectively. The isobaric closure (Allaire et al. Reference Allaire, Clerc and Kokh2002) is adopted, namely assuming that the mixture is in mechanical (pressure) equilibrium:

(3.5) \begin{equation} p = p_{l}(\rho _{l}, e_{l})=p_{g}(\rho _{g}, e_{g}). \end{equation}

Substituting the EOS of the liquid (2.5) and that of the plasma (2.7) into (3.4) under the pressure equilibrium condition of (3.5) readily yields the EOS for the mixture:

(3.6) \begin{equation} \left (\displaystyle \frac {\alpha _{l}}{\gamma _{l}-1} +\frac {\alpha _{g}}{\gamma _{g}-1}\right )p=\rho e-\left ( \displaystyle \frac {\alpha _{l}\gamma _{l}{P}_{\infty ,l}}{\gamma _{l}-1}+\frac {\alpha _{g}\gamma _{g}{P}_{\infty ,g}}{\gamma _{g}-1} \right )\!, \end{equation}

which can also be expressed in a neater form similar to the SG model as

(3.7) \begin{equation} p+\gamma {P}_{{\infty }}=(\gamma -1)\rho e, \end{equation}

where the specific heat ratio $\gamma$ and the reference pressure ${P}_{\infty }$ of the mixture are defined as

(3.8) \begin{equation} \left \{\begin{array}{l} \displaystyle \frac {1}{\gamma - 1} = \frac {\alpha _{l}}{\gamma _l - 1} + \frac {\alpha _{g}}{\gamma _g - 1}, \\[12pt] \displaystyle \frac {\gamma {P}_{\infty }}{\gamma - 1} = \frac {\alpha _{l} \gamma _l {P}_{\infty ,l}}{\gamma _l - 1} + \frac {\alpha _{g} \gamma _g {P}_{\infty ,g}}{\gamma _g - 1}. \end{array}\right . \end{equation}

The sound speed of the mixture ( $C$ ), which will be needed by the numerical solver discussed later in § 4 as a characteristic speed, is constructed in a similar way as in Allaire et al. (Reference Allaire, Clerc and Kokh2002):

(3.9) \begin{equation} \frac {\rho C^{2}}{\gamma -1}=\frac {\alpha _{l}\rho _{l}C_{l}^2}{\gamma _{l}-1}+\frac {\alpha _{g}\rho _{g}C_{g}^2}{\gamma _{g}-1}. \end{equation}

This formulation guarantees that the sound speed monotonically transits from one phase to the other across the diffuse interface region in the present work.

3.4. Fluxes and sources in the diffuse interface region

To complete (3.1), the energy fluxes and sources for the mixture (i.e. $\boldsymbol{q}$ , $\boldsymbol{F}_{\!r}$ and $S_r$ ) need to be formulated, in addition to the EOS. Furthermore, their expressions must satisfy the jump conditions described in § 3.1, i.e. these quantities vanish at the interfaces.

To obtain physically reasonable expressions for $S_r$ and $\boldsymbol{q}$ in the mixture energy equation (3.1), we start with the separate energy equation for each phase:

(3.10) \begin{align} \partial _{t} \big( \rho _{l}\alpha _{l} e_{l} \big) +\boldsymbol{\nabla }\boldsymbol{\cdot }\big( \rho _{l}\alpha _{l} e_{l} \boldsymbol{u} \big) +\alpha _{l}p\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = -\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\boldsymbol{q}_{ll}+\boldsymbol{q}_{lg}\big)+S_{r,l}+Q_{lg} , \end{align}

(3.11) \begin{align} \partial _{t} \big( \rho _{g}\alpha _{g} e_{g} \big) +\boldsymbol{\nabla }\boldsymbol{\cdot }\big( \rho _{g}\alpha _{g} e_{g} \boldsymbol{u} \big) +\alpha _{g}p\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = -\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\boldsymbol{q}_{gl}+\boldsymbol{q}_{gg}\big)+S_{r,g}+Q_{gl} , \end{align}

where the fluxes and source for the liquid are retained at this moment without losing generality; in other words, the liquid is also treated as a kind of plasma. Note that energy exchanges occur both within and between phases. Here, $\boldsymbol{q}_{ij}\ (i, j \in \{l, g\})$ denotes the thermal energy flux across the interface of any infinitesimal control volume, carried by the surrounding phase $j$ and applied to the internal phase $i$ , $S_{r,i}\ (i \in \{l, g\})$ is the radiation energy deposition rate to phase $i$ , and $Q_{ij}\ (i, j \in \{l, g\},\ i\ne j)$ is the energy exchange rate from internal phase $j$ to internal phase $i$ with $Q_{lg} = -Q_{gl}$ . For generality, we initially retain all terms for the liquid phase, effectively treating it as a plasma with its own transport coefficients – the thermal conductivity $\lambda _{l}$ , the radiative energy exchange rate $\omega _{r,l}$ , and the radiative energy diffusion coefficient $\chi _{r,l}$ – despite these being typically much smaller than their plasma-phase counterparts.

The radiative energy source for each phase is modelled as

(3.12) \begin{equation} S_{r,i}=\alpha _{i}\omega _{r,i}\big( T_{r}^{4}-T_{i}^{4} \big), \quad i \in \{l, g\}. \end{equation}

Therefore the source of the mixture is the sum of $S_{r,g}$ and $S_{r,l}$ , which reads

(3.13) \begin{equation} S_{r}=\alpha _{g}\omega _{r,g}\big( T_{r}^{4}-T_{g}^{4} \big)+\alpha _{l}\omega _{r,l}\big( T_{r}^{4}-T_{l}^{4} \big). \end{equation}

According to its definition, the plasma-to-plasma conductive heat flux $\boldsymbol{q}_{gg}$ should be proportional to the local area over which the internal plasma contacts the surrounding plasma. For a homogeneous mixture of plasma and liquid, the local fraction of contact area between the plasmas is equal to the product of the volume fractions on both sides, namely $\alpha _{g}^+ \alpha _{g}^-$ , where the superscripts $+$ and $-$ denote the internal and external sides of the contact surface, respectively. Given homogeneity, $\alpha _{g}=\alpha _{g}^+= \alpha _{g}^-$ and $\alpha _{l}=\alpha _{l}^+= \alpha _{l}^-$ . Therefore, we can have

(3.14) \begin{equation} \boldsymbol{q}_{gg} = -\alpha _{g}^2\lambda _{g}\boldsymbol{\nabla }T_{g}. \end{equation}

Similarly, the contact area fraction between the internal plasma and the surrounding liquid can be expressed as $\alpha _{g}^+ \alpha _{l}^-$ . The interphase thermal conduction is expected to be determined by the two phases, thus we phenomenologically assign the effective interphase thermal conduction coefficient as $\sqrt {\lambda _{g}\lambda _{l}}$ , a geometric average of $\lambda _{g}$ and $\lambda _{l}$ . Assuming a local temperature difference $T_{g}^{+}-T_{l}^{-}$ over a small length scale $\delta \boldsymbol{x}$ , the liquid-to-plasma flux $\boldsymbol{q}_{gl}$ can be estimated by

(3.15) \begin{equation} \boldsymbol{q}_{gl} = -\alpha _{g}^+\alpha _{l}^-\sqrt {\lambda _{g}\lambda _{l}} \big ( T_{g}^{+}-T_{l}^{-} \big ) /\delta \boldsymbol{x}. \end{equation}

Equation (3.15) can be reformulated in two equivalent ways:

(3.16) \begin{align} \boldsymbol{q}_{gl} &= -\alpha _{g}^+\alpha _{l}^-\sqrt {\lambda _{g}\lambda _{l}} \big(\boldsymbol{\nabla }T_{g}+\big(T_{g}^{-}-T_{l}^{-}\big) /\delta \boldsymbol{x} \big), \qquad \text{or}\nonumber \\[3pt] &= -\alpha _{g}^+\alpha _{l}^-\sqrt {\lambda _{g}\lambda _{l}} \big( \boldsymbol{\nabla }T_{l}+\big(T_{g}^{+}-T_{l}^{+}\big) /\delta \boldsymbol{x} \big). \end{align}

In these forms, the local temperature difference between the plasma and liquid phases ( $T_{g}^{+}-T_{l}^{+}$ or $T_{g}^{-}-T_{l}^{-}$ ) implies that $\boldsymbol{q}_{gl}$ is also related to the exchange rate of internal energy from plasma to liquid, $Q_{gl}$ , evaluated on the internal and external sides of the contact surface, respectively.

To keep the notational symmetry, the liquid-to-liquid and plasma-to-liquid energy fluxes due to thermal conduction can be expressed as

(3.17) \begin{equation} \boldsymbol{q}_{ll} = -\alpha _{l}^{2}\lambda _{l}\boldsymbol{\nabla }T_{l}, \quad \boldsymbol{q}_{lg} = -\alpha _{l}^+\alpha _{g}^-\sqrt {\lambda _{g}\lambda _{l}}\big( T_{l}^{+}-T_{g}^{-} \big) /\delta \boldsymbol{x}. \end{equation}

The mixture thermal flux is then contributed by all of the fluxes within and between phases, i.e. $\boldsymbol{q}=\boldsymbol{q}_{gg}+\boldsymbol{q}_{gl}+\boldsymbol{q}_{lg}+\boldsymbol{q}_{ll}$ . It is worth noting that the temperature differences between the plasma and liquid phases in $\boldsymbol{q}_{lg}$ and $\boldsymbol{q}_{gl}$ cancel upon summation. Consequently, the final expression for $\boldsymbol{q}$ can be simplified accordingly:

(3.18) \begin{equation} \boldsymbol{q} = -\alpha _{g}^{2}\lambda _{g}\boldsymbol{\nabla }T_{g}-\alpha _{g}\alpha _{l}\sqrt {\lambda _{g}\lambda _{l}}\,(\boldsymbol{\nabla }T_{g}+\boldsymbol{\nabla }T_{l})-\alpha _{l}^{2}\lambda _{l}\boldsymbol{\nabla }T_{l}. \end{equation}

This modelling approach has the desirable property that it reduces to Fick’s law $\boldsymbol{q}=- \lambda \boldsymbol{\nabla }T$ in the limit where the two phases are identical, thereby ensuring its consistency with the single-phase theory.

The energy flux due to radiation diffusion is analogous to the conductive heat flux to some degree. In the optically thick limit, where the radiation can be described as a diffusion field, the radiation diffusion fluxes are proportional to the contact area between phases and the respective radiation diffusion coefficients ( $\chi _{r,g}, \chi _{r,l}$ ). By direct analogy to the conductive heat flux form derived in (3.18), the radiation energy flux for the mixture can be written as

(3.19) \begin{equation} \boldsymbol{F}_{\!r} = -(\alpha _{g}^{2}\chi _{r,g}+2\alpha _{g}\alpha _{l}\sqrt {\chi _{r,g}\chi _{r,l}}+\alpha _{l}^{2}\chi _{r,l})\,\boldsymbol{\nabla }T_{r}^{4}. \end{equation}

Because of the comparatively negligible thermal conduction coefficient and high opacity of the liquid, we set $\lambda _{l}\approx 0$ , $\omega _{r,l}\approx 0$ and $\chi _{r,l}\approx 0$ in this paper. Accordingly, the radiative energy source, the conductive heat flux and the radiative diffusion flux yields a simplified model valid within the diffuse interface region:

(3.20) \begin{equation} \begin{cases} \boldsymbol{q}=-\alpha _{g}^{2}\lambda _{g} \boldsymbol{\nabla }T_{g},\\[5pt] S_{r}=\alpha _{g}\omega _{r,g}\big( T_{r}^{4}-T_{g}^{4} \big), \\[5pt] \boldsymbol{F}_{\!r}=-\alpha _{g}^{2}\chi _{r,g} \boldsymbol{\nabla }T_{r}^{4}. \end{cases} \end{equation}

This simplified formulation preserves a continuous transition for the energy fluxes and sources from the plasma to the liquid, governed by the plasma volume fraction $\alpha _{g}$ . At the same time, it satisfies the required jump conditions at the interfaces.

3.5. Summary of the model for plasma–liquid interactions

By applying the expressions of the energy fluxes and source in (3.20), we now rewrite the radiation two-phase flow model in (3.1) as

(3.21) \begin{equation} \begin{cases} \partial _{t} \alpha _{l}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\alpha _{l}= 0 , \\[5pt] \partial _{t} ( \rho _{l} \alpha _{l} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{l} \alpha _{l}\boldsymbol{u}) =0 ,\\[5pt] \partial _{t} ( \rho _{g} \alpha _{g} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho _{g} \alpha _{g}\boldsymbol{u}) =0 ,\\[5pt] \partial _{t} ( \rho \boldsymbol{u} ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho \boldsymbol{u}\boldsymbol{u} )+\boldsymbol{\nabla }\boldsymbol{\cdot }\left ((p+p_{r}) \boldsymbol {I}\right ) =\sigma \kappa\, \boldsymbol{\nabla }\phi , \\[5pt] \partial _{t} ( \rho e ) +\boldsymbol{\nabla }\boldsymbol{\cdot }( \rho e \boldsymbol{u} ) +p\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}= \boldsymbol{\nabla }\boldsymbol{\cdot }\big(\alpha _{g}^{2}\lambda _{g} \boldsymbol{\nabla }T_{g}\big)+ \alpha _{g}\omega _{r,g}\big( T_{r}^{4}-T_{g}^{4} \big) , \\[5pt] \partial _{t} E_{r} +\boldsymbol{\nabla }\boldsymbol{\cdot }( E_{r} \boldsymbol{u} ) +p_{r}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\alpha _{g}^{2}\chi _{r,g} \boldsymbol{\nabla }T_{r}^{4}\big) -\alpha _{g}\omega _{r,g}\big( T_{r}^{4}-T_{g}^{4} \big). \end{cases} \end{equation}

Although this model is derived for the mixture, it exhibits the correct asymptotic behaviour by reducing to standard single-phase models in the limit of $\alpha _{g}$ . It is straightforward that the system of equations in (3.21) converges to the Euler equation (2.3) for the pure liquid in the limit of vanishing plasma and radiation, while it recovers the RHD equations for a pure plasma (2.1) as $\alpha _{g} \rightarrow 1$ .

4.1. Hyperbolic step: convection of fluid and radiation

The convection equations (4.1) are discretised using a second-order finite-volume scheme with the numerical fluxes computed via the HLLC approximate Riemann solver (Toro Reference Toro2013). The characteristic wave speed $C_{\!s}$ is obtained from the eigenvalue of the hyperbolic system, computed from its Jacobian matrix; see Appendix B for details. The term for velocity divergence in the volume-fraction advection equation (4.1a ) is evaluated using the adapt-HLLC scheme (Johnsen & Colonius Reference Johnsen and Colonius2006). The surface tensions in (4.1d ) and (4.1e ) are treated as source terms and discretised by the central difference scheme. The time step $\Delta t$ is chosen adaptively such that the Courant–Friedrichs–Lewy (CFL) condition is strictly satisfied:

(4.3) \begin{equation} \Delta t=\mathrm{CFL}\times\min _{i\in D} \left \{ \frac {\sqrt [3]{V_{i}} }{\left | \boldsymbol{u}_{i} \right | +C_{s,i}} \right \} \!, \end{equation}

where the $\mathrm{CFL}$ number ( $=0.44$ ) is a prescribed constant, the subscript $i$ denotes a Cartesian cell in the computational domain $D$ , $C_{\!s}=\sqrt {C^{2}+4p_{r}/(3\rho )}$ is the local wave speed, and $V$ is the volume of the cell.

After solving the fluid and radiation convection explicitly, we update the conservative variables to an intermediate state, denoted by the superscript $*$ : $[\alpha_l^*, (\rho_l\alpha_l)^*, (\rho_g\alpha_g)^*, (\rho\boldsymbol{u})^*, (\rho e_{t})^*]^{\rm T}$ . The corresponding primitive variables at the intermediate state are also updated. From $(\rho e_{t})^*$ , the internal energy $(\rho e)^*$ and radiation energy $(E_{r})^*$ of the mixture can then be calculated; see details in Appendix C.

4.2. Parabolic step: thermal conduction and radiation diffusion

An implicit scheme is adopted for the time integration of the thermal conduction and radiation diffusion in (4.2). All the coefficients in the equations are linearised using the flow variables at the intermediate state:

(4.4) \begin{equation} \begin{cases} \displaystyle \frac {(\rho e)^{n+1}-(\rho e)^{*}}{\Delta t} = \boldsymbol{\nabla }\boldsymbol{\cdot }\big(\big(\alpha _{g}^{2}\lambda _{g}\big)^{*}\,\boldsymbol{\nabla }(T_{g})^{n+1}\big)+(\alpha _{g}\omega _{r,g})^{*} \,\big( T_{r}^{4}-T_{g}^{4} \big)^{n+1},\\[12pt] \displaystyle \frac { (E_{r} )^{n+1}- (E_{r} )^{*}}{\Delta t} =\boldsymbol{\nabla }\boldsymbol{\cdot }\big( \big(\alpha _{g}^{2}\chi _{r,g} \big)^{*}\,\boldsymbol{\nabla }\big(T_{r}^{4} \big)^{n+1}\big) + (\alpha _{g}\omega _{r,g})^{*}\, \big(T_{g}^{4}-T_{r}^{4}\big)^{n+1}. \end{cases} \end{equation}

The gradient and divergence terms are discretised using a central finite-difference scheme. By solving (4.4) iteratively, $\rho e$ and $E_{r}$ are updated from the intermediate state to the $n+1$ time step. Subsequently, by employing the EOS, all remaining conservative and primitive variables associated with energy, such as $p$ , $p_{r}$ , $T_{g}$ and $T_{r}$ , are calculated for the $n+1$ time step, while those related to mass and momentum are retained from the intermediate state. This completes the numerical solution of the radiation two-phase flow system for plasma–liquid interactions.

5.1. Radiative shock tube

Figure 2.

The profiles of a radiative shock problem at $0\ \mathrm{ns}$ (black solid lines and open squares) and $4.2\ \mathrm{ns}$ (red dashed lines and solid circles), in terms of (a) density, (b) plasma temperature, and (c) radiation temperature. The symbols and lines denote numerical and semi-analytical solutions, respectively. Note that the symbols are sampled at every $8$ th grid point for visual clarity.

To validate the radiation transport in our model, we consider a one-dimensional radiative shock tube problem, where the radiation energy fluxes and radiation pressure significantly affect the hydrodynamics. As illustrated in figure 2, a low-temperature low-density plasma (with density $\rho ^{L}$ , temperature $T^{L}$ and velocity $u^{L}$ ) is initially impacting a high-temperature high-density plasma (with density $\rho ^{R}$ , temperature $T^{R}$ and velocity $u^{R}$ ). This configuration is designed to produce a stationary radiative shock at $x=0 \,\unicode{x03BC} \mathrm{m}$ , the solution of which can be compared with the semi-analytical result from Lowrie & Edwards (Reference Lowrie and Edwards2008).

Because the fluid system only has one phase, i.e. plasma, we set $\alpha _{g}=1$ everywhere in the domain. The properties of the plasma are an adiabatic index $\gamma _{g} = 5/3$ , a molar mass $M_{g} = 2\ \mathrm{g\,mol}^{-1}$ , and a constant ionisation degree $Z_{g} = 1$ , which determine $c_{v,g} = 1.247\times 10^8\ \mathrm{cm}^{2}\ \mathrm{s}^{-2}\ \mathrm{K}^{-1}$ . Thermal conduction is neglected (i.e. $\lambda _{g}=0$ ). The diffusion coefficient and the exchange rate of the radiative energy are $\chi _{r,g} = 1.268\times 10^{7}\ \mathrm{cm}^{2}\,\mathrm{s}^{-1}$ and $\omega _{r,g} = 1.268\times 10^{13}\ \mathrm{s}^{-1}$ , respectively. The initial conditions are defined by a step function at $x=0\, \unicode{x03BC} \mathrm{m}$ , as illustrated in figure 2. The radiation is initially set in thermal equilibrium with the plasma (i.e. $T_{r}=T$ ). The initial flow parameters of the plasma are

(5.1) \begin{equation} \begin{cases} \rho ^{R} = 2.286\ \mathrm{g}\,\mathrm{cm}^{-3},\quad T^{R} = 2.411\times 10^6\ \mathrm{K},\quad u^{R} = 1.11\times 10^7\ \mathrm{cm}\,\mathrm{s}^{-1}, \\[3pt] \rho ^{L} = 1\ \mathrm{g}\,\mathrm{cm}^{-3},\quad T^{L} = 1.16\times 10^6\ \mathrm{K},\quad u^{L} = 2.536\times 10^7\ \mathrm{cm}\,\mathrm{s}^{-1}. \end{cases} \end{equation}

The computation is performed in a domain $ [-300,\ 300 ]\,\unicode{x03BC} \mathrm{m}$ , discretised by a uniform grid of $2000$ points.

Figure 2 presents numerical results at $t = 4.2\ \mathrm{ns}$ in terms of density, plasma temperature and radiation temperature. These results are superimposed with the semi-analytical solution by Lowrie & Edwards (Reference Lowrie and Edwards2008). Obviously, the numerical results are in excellent agreement with the theoretical prediction, with respect to the position of the stationary shock and the diffusion profiles on either side of the discontinuity. The comparison demonstrates that the proposed model can be reduced to an RHD model for single phase, and successfully captures complex plasma flow structures, including radiative shocks.

5.2. Blast wave with thermal conduction

Figure 3.

(a) Density, (b) relative density error $|\rho ^{N}-\rho ^{A}|/\rho ^{A}$ , (c) temperature, and (d) relative temperature error $|T^{N}-T^{A}|/T^{A}$ for the Reinicke–Meyer-ter-Vehn blast wave problem at $0.52\ \mathrm{ns}$ , where superscripts $N$ and $A$ represent numerical and semi-analytical results, respectively.

To validate the thermal conduction in our model, we consider the Sedov–Taylor point explosion with nonlinear thermal conduction, also referred to as the Reinicke–Meyer-ter-Vehn problem, which only involves the plasma phase and has the semi-analytical solution for the blast wave (Reinicke & Meyer-ter-Vehn Reference Reinicke and Meyer-ter-Vehn1991). To avoid the singularity of a point explosion, the initial condition for the simulation is taken from the semi-analytical solution right after the bang time. At this moment, the stationary and uniform ambient plasma is about to be swept through by a tiny and rapidly expanding spherical blast wave.

The plasma is modelled as an ideal gas with adiabatic index $\gamma _{g} = 1.25$ , molar mass $M_{g} = 1\ \mathrm{g\,mol}^{-1}$ and ionisation degree zero, giving a specific heat capacity at constant volume $c_{v,g} = 3.326\times 10^8\ \mathrm{cm}^{2}\ \mathrm{s}^{-2}\ \mathrm{K}^{-1}$ . The radiation transport is not excluded from this problem; accordingly, the radiation energy $E_{r}$ , the radiation pressure $p_{r}$ , the radiative energy exchange rate $\omega _{r,g}$ , and the radiative energy diffusion coefficient $\chi _{r,g}$ for the plasma are all set to zero. The thermal conduction coefficient follows a nonlinear power-law function of density and temperature as recommended by Reinicke & Meyer-ter-Vehn (Reference Reinicke and Meyer-ter-Vehn1991), i.e. $\lambda _{g} = \rho _g^{-2}T_g^{6.5}\ (\text{in }\mathrm{g}\boldsymbol{\, }\mathrm{cm}\boldsymbol{\, }\mathrm{s}^{-3}\rm\, K^{-1})$ , where $\rho _{g}$ is in $\mathrm{g}\,\mathrm{cm}^{-3}$ , and $T_{g}$ is in $\mathrm{K}$ . Axisymmetric simulations are performed in 2-D cylindrical coordinates ( $(r,z)$ plane). The initial condition corresponds to the semi-analytical solution at $t = 0.2\ \mathrm{ns}$ , when the shock front has radius $0.225\ \mathrm{cm}$ , and the heat front has radius $0.45\ \mathrm{cm}$ . The computational domain is $ [ 0, 1 ] \times [ 0, 1 ]\ \mathrm{cm}^2$ , and is discretised by $2048\times 2048$ grid points.

Figures 3( $a$ ) and 3( $c$ ) illustrate the distributions of density and temperature at $t = 0.52\ \mathrm{ns}$ . As predicted by Reinicke & Meyer-ter-Vehn (Reference Reinicke and Meyer-ter-Vehn1991), the shock front reached radius $0.45\ \mathrm{cm}$ and the heat front $0.9\ \mathrm{cm}$ . The corresponding relative errors in density and temperature are shown in figures 3( $b$ ) and 3( $d$ ), respectively. The maximum errors, which do not exceed $5\,\%$ , are localised at the discontinuities of the shock and heat front, and are attributed to finite grid resolution. The good agreement with the semi-analytical solution demonstrates the capability of the proposed model to accurately deal with nonlinear thermal conduction.

5.3. Bubble collapse induced by shock wave

The shock-induced collapse of a gas bubble in water, where a planar shock wave propagating through the water interacts with an underwater gas bubble, has been extensively studied by simulations (Nourgaliev, Dinh & Theofanous Reference Nourgaliev, Dinh and Theofanous2006; Hawker & Ventikos Reference Hawker and Ventikos2012; Bo & Grove Reference Bo and Grove2014). This classic problem is used here to verify that the governing equations of the proposed model correctly reduce to those for compressible two-phase flows. In this case, thermal conduction, radiation transport and surface tension are neglected. Consequently, all associated variables and coefficients – including radiation energy $E_{r}$ , radiation pressure $p_{r}$ , thermal conduction coefficient $\lambda _{g}$ , radiative energy diffusion coefficient $\chi _{r,g}$ , radiative energy exchange rate $\omega _{r,g}$ , and surface tension coefficient $\sigma$ – are set to zero. Regarding the EOS, the water is described by the SG model (parameters are listed in table 1), while the gas phase is treated as an ideal gas with $\gamma _{g}=1.4$ .

The simulation is performed in 2-D Cartesian coordinates. Initially, a circular bubble is centred at $ ( 0, 0 )\ \mathrm{mm}$ with radius $3\ \mathrm{mm}$ . The density and pressure inside the bubble are $1\ \mathrm{kg}\,\mathrm{m}^{-3}$ and $10^{5}\ \mathrm{Pa}$ , respectively. An incident shock wave is positioned at $5.4\ \mathrm{mm}$ to the left of the bubble centre initially, and propagates to the $x+$ direction. The right-hand (with the superscript $R$ ) and left-hand (with the superscript $L$ ) sides of the shock are initialised as follows:

(5.2) \begin{equation} \begin{cases} \rho ^{R} = 1000\ \mathrm{kg}\,\mathrm{m}^{-3},\quad p^{R} = 10^{5}\ \mathrm{Pa},\quad u^{R} = 0, \\[3pt] \rho ^{L} = 1323.65\ \mathrm{kg}\,\mathrm{m}^{-3},\quad p^{L} = 1.9\times 10^{9}\ \mathrm{Pa},\quad u^{L} = 681.58\ \mathrm{m}\,\mathrm{s}^{-1}. \end{cases} \end{equation}

The computational domain spans $ [ -15, 15 ]\times [ -15, 15 ]\ \mathrm{mm}^2$ , and the bubble radius ( $3\ \mathrm{mm}$ ) is initially discretised by $800$ grid points.

Figure 4 demonstrates different stages of the shock-induced bubble collapse in the simulation, visualised by pressure contours and numerical schlieren. Upon shock impact, a reflected rarefaction wave forms in the liquid, and a transmitted shock develops inside the bubble, as shown in figure 4( $a$ ). At the same time, a re-entrant liquid jet forms at the left-hand side of the bubble, and propagates in the same direction as the incident shock, as shown in figure 4( $b$ ). Upon impact of the jet with the rear interface of the bubble, a water-hammer shock is generated (see e.g. $t=3.9 \,\unicode{x03BC} \mathrm{s}$ in figure 4 $c$ ). The characteristic parameters of this event include the collision time $t_c$ , the water-hammer shock pressure $p_w$ , and the jet speed $v_{\!j}$ at the instant of impact. Subsequently, the jet penetrates the bubble and merges with the surrounding liquid, while secondary jets gradually develop (see figure 4 $d$ ). These flow features closely resemble the numerical results reported in Hawker & Ventikos (Reference Hawker and Ventikos2012). A quantitative comparison on $t_c$ , $v_{\!j}$ and $p_w$ with the literature is provided in table 2. Good agreement is achieved for $t_c$ and $v_{\!j}$ , while $p_w$ is relatively lower in the present simulation, primarily due to coarser spatial resolution. The results confirm that the proposed radiation two-phase flow model correctly reduces to a standard compressible two-phase flow formulation when the radiation-specific terms are disabled. Furthermore, the good agreement with expected flow dynamics validates the core flow solver and the implementation of the EOS before the coupling with radiation transport and thermal conduction is introduced.

Table 2.

Shock-induced bubble collapse compared to prior simulations: collision time ( $t_{c}$ ), jet speed at collision ( $v_{\!j}$ ), and water-hammer shock pressure ( $p_{w}$ ).

Table 2. Long description

A table comparing shock-induced bubble collapse parameters across different studies. The table has five rows and three columns. The columns are labeled t subscript c (microseconds), v subscript j (meters per second), and p subscript wh (gigapascals). The rows are labeled Present, Lin et al. (2017), Hawker & Ventikos (2012), Bo & Grove (2014), and Nourgaliev et al. (2006). Row 1: Present, t subscript c 3.70, v subscript j 2837, p subscript wh 5.03. Row 2: Lin et al. (2017), t subscript c 3.70, v subscript j 2832, p subscript wh 5.90. Row 3: Hawker & Ventikos (2012), t subscript c 3.66, v subscript j 2810, p subscript wh 5.89. Row 4: Bo & Grove (2014), t subscript c 3.70, v subscript j 2830, p subscript wh not provided. Row 5: Nourgaliev et al. (2006), t subscript c 3.69, v subscript j 2850, p subscript wh 10.1.

Figure 4.

Numerical results of bubble collapse induced by a planar shock with respect to pressure (upper half) and numerical schlieren (lower half) at (a) $2.2\,\unicode{x03BC} \mathrm{s}$ , (b) $3.7\,\unicode{x03BC} \mathrm{s}$ , (c) $3.9\,\unicode{x03BC} \mathrm{s}$ , (d) $4.1 \,\unicode{x03BC} \mathrm{s}$ . The red lines represent the bubble interfaces (by $\alpha _{g} = 0.5$ ).

5.4. Droplet oscillation driven by surface tension

Due to the imbalance in surface tension, an initially elliptical droplet oscillates periodically, converting between potential and kinetic energy (Garrick et al. Reference Garrick, Owkes and Regele2017). The oscillation period can be analytically obtained (Fyfe, Oran & Fritts Reference Fyfe, Oran and Fritts1988):

(5.3) \begin{equation} t_{os}=2\pi \sqrt {\frac {(\rho _{l}+\rho _{g})R^{3}}{\big( o^{3}-o \big)\sigma }}, \end{equation}

where $o$ is the oscillation mode, and $R$ is the radius of the droplet at equilibrium. We employ this case to validate the surface tension model.

Figure 5.

Droplet oscillation: ( $a$ ) evolution of the droplet over an oscillation period, where $A$ and $B$ denote the initial semi-axes; ( $b$ ) kinetic energy of the droplet.

In the droplet oscillation, thermal conduction and radiation transport are absent; therefore, all corresponding variables and coefficients ( $E_{r}, p_{r}, \lambda _{g}, \chi _{r,g}, \omega _{r,g}$ ) are set to zero in the simulation. For the EOS, the droplet is modelled as water by the SG model (see table 1), while the gas phase is treated as an ideal gas with $\gamma _{g}=1.33$ . The surface tension coefficient of the droplet is fixed at $\sigma = 500\ \mathrm{mN\,m}^{-1}$ . The simulation is performed in 2-D Cartesian coordinates. Initially, the elliptical droplet is centred at $ ( 0, 0 )\ \unicode{x03BC} \mathrm{m}$ with the semi-axes $A = 1.2\,\unicode{x03BC} \mathrm{m}$ and $B = 0.8\,\unicode{x03BC} \mathrm{m}$ , resulting in droplet radius $R = 0.9798\,\unicode{x03BC} \mathrm{m}$ at equilibrium. The initial pressure of both the droplet and the ambient gas is set at $1.3\times 10^{5}\ \mathrm{Pa}$ , while the densities of the two phases are $1073.55\ \mathrm{kg}\,\mathrm{m}^{-3}$ and $0.74\ \mathrm{kg}\,\mathrm{m}^{-3}$ , respectively. To exploit symmetry of the oscillation, only a quarter of the physical domain ( $ [ 0, 5 ]\times [ 0, 5 ]\,\unicode{x03BC} \mathrm{m}^2$ ) is simulated, with uniform mesh resolution $2\ \text{nm}$ .

Figure 5( $a$ ) shows snapshots of the droplet shapes at different times over one oscillation period, corresponding to the phase values $0$ , $\pi /2$ , $\pi$ , $3\pi /2$ and $2\pi$ respectively. An oscillation mode $o=2$ is observed. The periodicity of the droplet oscillation is quantitatively illustrated in figure 5( $b$ ), in terms of the kinetic energy ( $\int (\rho _{l}\alpha _{l}\left | \boldsymbol{u} \right | ^{2}/2)\,\mathrm{d}x\,\mathrm{d}y$ ) versus time. A gradual decay of the oscillation amplitude is observed due to numerical dissipation, consistent with the findings in other studies (Garrick et al. Reference Garrick, Owkes and Regele2017; Long, Cai & Pan Reference Long, Cai and Pan2023). The oscillation period in simulation is $118.4\ \mathrm{ns}$ , with relative error $2.7\,\%$ compared to the theoretical prediction $115.3\ \mathrm{ns}$ . This demonstrates that the continuous surface tension incorporated in the present model accurately captures capillary effects.