2.1 Scaling analysis

We consider a spherical droplet with radius $R$ and density $\unicode[STIX]{x1D70C}_{0}$ that is subjected to an ablation pressure on the illuminated side with an amplitude $p_{e}$ and a duration $\unicode[STIX]{x1D70F}_{e}$ , see figure 2. The total impulse received by the droplet is given by $J\sim p_{e}\unicode[STIX]{x1D70F}_{e}$ . In this work we explore the effect of decreasing the pulse duration while keeping the total impulse transferred to the droplet constant; i.e. decreasing $\unicode[STIX]{x1D70F}_{e}$ at constant $J$ .

Figure 2.

Sketch of the problem: an ablation pressure of amplitude $p_{e}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)$ and duration $\unicode[STIX]{x1D70F}_{e}$ acts on the surface of a droplet with radius $R$ . As a result, a pressure field is induced inside the droplet over a depth $\ell _{e}\sim c\unicode[STIX]{x1D70F}_{e}$ . The colour bar denotes the pressure amplitude (blue is ambient and yellow is peak pressure). The spherical coordinate system ( $r$ , $\unicode[STIX]{x1D703}$ , $\unicode[STIX]{x1D719}$ ) used is indicated (the azimuthal angle $\unicode[STIX]{x1D719}$ is not shown but rotates around the $z$ -axis).

During the pulse the pressure disturbance on the surface of the droplet penetrates over a length scale $\ell _{e}\sim c\unicode[STIX]{x1D70F}_{e}$ , where $c$ is the speed of sound in the droplet. If $\ell _{e}$ is short compared to $R$ all momentum is initially concentrated inside a thin layer, as is illustrated in figure 2. By contrast, if $\ell _{e}>R$ , all fluid inside the droplet has experienced a change in momentum directly after the pulse. The ratio between $\ell _{e}$ and $R$ is quantified by the acoustic Strouhal number and is a dimensionless pressure pulse duration

(2.1) $$\begin{eqnarray}St=\frac{\ell _{e}}{R}=\frac{c\unicode[STIX]{x1D70F}_{e}}{R}.\end{eqnarray}$$

To investigate the effect of short pulse durations, we are interested in the limit $St\lesssim 1$ .

Figure 3.

A phase diagram showing lines of constant impulse transfer to the droplet. The blue dotted-dashed line is the isoline $St\,Ma=0.01$ , the orange dashed line is the isoline $St\,Ma=0.1$ and the red dotted line is the isoline $St\,Ma=0.5$ . The plot shows three distinct regimes can be observed at constant impulse: a strongly compressible regime, a weakly compressible regime and an incompressible regime. The weakly compressible regime is the focus of the present work.

When $\unicode[STIX]{x1D70F}_{e}$ is decreased at constant $J$ , $p_{e}$ rises. From momentum conservation in this thin layer it follows that the typical velocity induced inside $\ell _{e}$ is given by $u_{e}\sim p_{e}/(\unicode[STIX]{x1D70C}_{0}c)$ , where $\unicode[STIX]{x1D70C}_{0}$ is the density of the liquid droplet. Hence we observe that a large $p_{e}$ induces large velocities in $\ell _{e}$ , which is quantified by the acoustic Mach number and is a dimensionless pressure pulse amplitude

(2.2) $$\begin{eqnarray}Ma=\frac{u_{e}}{c}=\frac{p_{e}}{p_{0}},\end{eqnarray}$$

where $p_{0}=\unicode[STIX]{x1D70C}_{0}c^{2}$ is the base pressure of the droplet. When $Ma$ is large the fluid response inside the droplet is nonlinear and shock waves dominate the flow. If $Ma$ small, the flow inside the droplet can be considered linear.

The product $Ma\,St$ sets the total dimensionless impulse received by the droplet

(2.3) $$\begin{eqnarray}St\,Ma=\frac{p_{e}\unicode[STIX]{x1D70F}_{e}}{\unicode[STIX]{x1D70C}_{0}Rc}=\frac{U}{c},\end{eqnarray}$$

where (1.1) is used to express the centre-of-mass velocity $U$ of the droplet as whole. This product is often referred to as the global Mach number of the droplet.

One can use $Ma$ and $St$ to delineate different regimes in the droplet response, as illustrated in figure 3. For lines of constant $Ma\,St$ (hence constant impulse), we can identify three regimes. Firstly, when $St$ is small and $Ma$ is large, we are in a strongly compressible regime where nonlinear advective acceleration and nonlinear viscous dampening need to be taken into account to describe the flow. Secondly, for intermediate $Ma$ and $St$ , compressible effects are important but nonlinear effects are small, which renders this regime analytically accessible. This regime, which we term the weakly compressible regime, will be the main focus of this paper. Finally, when $St\gg 1$ and $Ma\ll 1$ , we enter the incompressible regime that was subject of previous studies where long pulse durations (large $St$ ) were considered (Gelderblom et al. Reference Gelderblom, Lhuissier, Klein, Bouwhuis, Lohse, Villermaux and Snoeijer2016). This regime is also relevant to droplet impact studies on rigid surfaces, see e.g. Richard, Clanet & Quere (Reference Richard, Clanet and Quere2002), Clanet et al. (Reference Clanet, Beguin, Richard and Quere2004), Yarin (Reference Yarin2006), Josserand & Thoroddsen (Reference Josserand and Thoroddsen2016), Philippi, Lagree & Antkowiak (Reference Philippi, Lagree and Antkowiak2016), Wildeman et al. (Reference Wildeman, Visser, Sun and Lohse2016).

2.2 The weakly compressible model

The compressible flow equations are given by (Batchelor Reference Batchelor1967, p. 164)

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D70C}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\left(\unicode[STIX]{x1D705}+\frac{1}{3}\unicode[STIX]{x1D707}\right)\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the density field, $\boldsymbol{u}$ the velocity field, $p$ the thermodynamic pressure, $\unicode[STIX]{x1D707}$ is the shear viscosity coefficient and $\unicode[STIX]{x1D705}$ the expansion viscosity coefficient. In this work we only focus on the flow inside the droplet and we do not consider the flow of the outer gas phase, since for typical experimental situations (tin droplet in vacuum and water droplet in air): $\unicode[STIX]{x1D70C}_{g}/\unicode[STIX]{x1D70C}_{l}\ll 1$ and $\unicode[STIX]{x1D707}_{g}/\unicode[STIX]{x1D707}_{l}\ll 1$ (Klein et al. Reference Klein, Bouwhuis, Visser, Lhuissier, Sun, Snoeijer, Villermaux, Lohse and Gelderblom2015; Kurilovich et al. Reference Kurilovich, Klein, Torretti, Lassise, Hoekstra, Ubachs, Gelderblom and Versolato2016). In order to find a closed analytic expression for the pressure and velocity field inside the droplet, we will first need to simplify (2.4) and (2.5) and finally solve the system obtained with appropriate surface boundary conditions for the droplet.

In the weakly compressible regime, i.e. for intermediate $Ma$ , $St$ (see figure 3), we can expand the pressure $p$ , density $\unicode[STIX]{x1D70C}$ and velocity $\boldsymbol{u}$ inside the droplet as a constant plus a small time-dependent part, analogous to Batchelor (Reference Batchelor1967, p. 166)

(2.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle p(\boldsymbol{x},t)=p_{0}+p_{1}(\boldsymbol{x},t),\\ \displaystyle \unicode[STIX]{x1D70C}(\boldsymbol{x},t)=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}_{1}(\boldsymbol{x},t),\\ \displaystyle \boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}_{1}(\boldsymbol{x},t),\end{array}\right\}\end{eqnarray}$$

where $t$ is the time, $\boldsymbol{x}=(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ are the spherical coordinates defined in figure 2. The thermodynamic pressure can formally be expressed as $p=p(\unicode[STIX]{x1D70C},S,T)$ , where $S$ is the entropy and $T$ is the temperature. In this work however, we neglect any thermoelastic waves that result from thermal expansion due to the local heating of the surface of the droplet. This simplification is supported by typical experiments found in the literature, where both the thermal expansion coefficient and the thermal diffusion coefficient are small (Klein et al. Reference Klein, Bouwhuis, Visser, Lhuissier, Sun, Snoeijer, Villermaux, Lohse and Gelderblom2015; Kurilovich et al. Reference Kurilovich, Klein, Torretti, Lassise, Hoekstra, Ubachs, Gelderblom and Versolato2016). We do note that with increasing laser intensity, thermoelastic waves could become increasingly important (Sigrist & Kneubuhl Reference Sigrist and Kneubuhl1978; Wang & Xu Reference Wang and Xu2001). However, this is typically accompanied by an increase in the acoustic Mach number which is outside the scope of this work. Following these simplifications, we can write a basic equation of state for the fluid inside the droplet

(2.7) $$\begin{eqnarray}p(\boldsymbol{x},t)=p_{0}+c^{2}\unicode[STIX]{x1D70C}_{1}(\boldsymbol{x},t),\end{eqnarray}$$

where $c$ is the speed of sound and $p_{0}=c^{2}\unicode[STIX]{x1D70C}_{0}$ is the ideal base pressure of the stationary droplet. Since we are interested in the flow directly after the pulse we introduce the following non-dimensionalization, following the scaling analysis of the previous paragraph

(2.8a-e ) $$\begin{eqnarray}\boldsymbol{u}=\frac{p_{e}}{\unicode[STIX]{x1D70C}_{0}c}\tilde{\boldsymbol{u}},\quad \boldsymbol{x}=R\tilde{\boldsymbol{x}},\quad t=\unicode[STIX]{x1D70F}_{e}\tilde{t},\quad p=p_{e}\tilde{p},\quad \unicode[STIX]{x1D70C}=\frac{p_{e}}{c^{2}}\tilde{\unicode[STIX]{x1D70C}},\end{eqnarray}$$

where the tildes refer to the dimensionless parameters. From now on we drop the tildes and work with the dimensionless parameters. A consistent approximation of (2.4) and (2.5) is given by the system (Blackstock Reference Blackstock2000, p. 97)

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p_{1}}{\unicode[STIX]{x2202}t}+St(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{1})=0, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{St}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{1}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p_{1}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{1}+\frac{1}{Re_{v}}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{1}), & \displaystyle\end{eqnarray}$$

where $Re=\unicode[STIX]{x1D70C}_{0}Rc/\unicode[STIX]{x1D707}$ is the Reynolds number with $\unicode[STIX]{x1D707}$ the dynamic viscosity and $Re_{v}=\unicode[STIX]{x1D70C}_{0}Rc/((\unicode[STIX]{x1D707}+\unicode[STIX]{x1D705})/3)$ the Reynolds number for volume changes, where $\unicode[STIX]{x1D705}$ is the bulk viscosity. Although the Reynolds number in experiments is typically large ( $Re\sim 10^{3}$ ) we will see later on that we need to retain the viscous terms in (2.10) to overcome singularities when converging pressure waves superimpose in the centre of the droplet. In the high Reynolds regime viscous effects however do not influence the droplet deformation, which was already reported in a previous incompressible and inviscid model by Gelderblom et al. (Reference Gelderblom, Lhuissier, Klein, Bouwhuis, Lohse, Villermaux and Snoeijer2016). We note that there is no explicit dependence on the acoustic Mach number in (2.9) and (2.10), since this is a low-order Mach expansion of the compressible Navier–Stokes equations and we scaled all our fields accordingly (2.8). By taking the divergence of (2.10) and using (2.9) we obtain a viscous wave equation for the acoustic field inside the droplet

(2.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}p_{1}}{\unicode[STIX]{x2202}t^{2}}-St^{2}\unicode[STIX]{x1D6FB}^{2}p_{1}=\frac{1}{Re_{a}}\left[\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}p_{1}}{\unicode[STIX]{x2202}t}\right],\end{eqnarray}$$

where $1/Re_{a}=St(1/Re+1/Re_{v})$ is an effective Reynolds number for the viscous dissipation in the acoustic wave (Blackstock Reference Blackstock2000, p. 97). Below we describe how this equation for $p_{1}(\boldsymbol{x},t)$ is solved for the problem at hand. In § 2.2.2 we show how the velocity field $\boldsymbol{u}_{1}(\boldsymbol{x},t)$ can be computed once $p_{1}(\boldsymbol{x},t)$ is known.

2.2.1 The pressure field

We solve (2.11) subject to a pressure boundary condition on the droplet surface. At the interface of the droplet the stresses between the liquid and the gas phase must be continuous. We assume that the magnitude of the pressure variations in the gas phase are much smaller than those inside the droplet. For typical experiments (Klein et al. Reference Klein, Bouwhuis, Visser, Lhuissier, Sun, Snoeijer, Villermaux, Lohse and Gelderblom2015; Kurilovich et al. Reference Kurilovich, Klein, Torretti, Lassise, Hoekstra, Ubachs, Gelderblom and Versolato2016) this assumption is justified since the density and the viscosity of the gas phase are much smaller. As a result, the stress condition on the interface significantly simplifies since the pressure in the gas phase may be considered constant and equal to $p_{0}$ . As a consequence, the time-dependent part of the pressure at the surface of the droplet must satisfy

(2.12) $$\begin{eqnarray}p_{1}(1,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=0,\end{eqnarray}$$

to ensure that all acoustic energy remains inside the droplet. Furthermore we assumed that the interface remains immobile and therefore the droplet spherical during the pulse. This latter assumption is justified when the pulse duration $\unicode[STIX]{x1D70F}_{e}$ is much smaller than the typical interface deformation time scale $\unicode[STIX]{x1D70F}_{int}=R/u_{e}$ , or in dimensionless form $St\,Ma\ll 1$ . Typically in experiments $St\,Ma\sim 10^{-2}-10^{-1}$ . An ablation pressure acting on the surface of the droplet is introduced through a Green’s function formalism (Morse & Feshbach Reference Morse and Feshbach1953, chap. 7). The general solution for the spatio-temporal pressure field inside the droplet is given by

(2.13) $$\begin{eqnarray}\displaystyle p_{1}(\boldsymbol{x},t) & = & \displaystyle \iiint _{V}\left[\frac{\unicode[STIX]{x2202}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t_{0}}p_{1}(\boldsymbol{x}_{0},t_{0})\right.\nonumber\\ \displaystyle & & \displaystyle \left.\left.-\,\left(\frac{\unicode[STIX]{x2202}p_{1}(\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t_{0}}-\frac{p_{1}(\boldsymbol{x}_{0},t_{0})}{Re_{a}}\unicode[STIX]{x1D6FB}_{\boldsymbol{x}_{0}}^{2}\right)G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})\right]\right|_{t_{0}=0}^{t^{+}}\,\text{d}V_{0}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{0}^{t^{+}}\,\text{d}t_{0}\oint _{S}\left[G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})\left(St^{2}\unicode[STIX]{x1D735}_{\boldsymbol{x}_{0}}p_{1}(\boldsymbol{x}_{0},t_{0})+\frac{1}{Re_{a}}\unicode[STIX]{x1D735}_{\boldsymbol{x}_{0}}\frac{\unicode[STIX]{x2202}p_{1}(\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t_{0}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,\unicode[STIX]{x1D735}_{\boldsymbol{x}_{0}}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})\left(St^{2}p_{1}(\boldsymbol{x}_{0},t_{0})+\frac{1}{Re_{a}}\frac{\unicode[STIX]{x2202}p_{1}(\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t_{0}}\right)\right]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S_{0},\end{eqnarray}$$

where $G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})$ is the Green’s function satisfying

(2.14) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t^{2}}-St^{2}\unicode[STIX]{x1D6FB}^{2}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})-\frac{1}{Re_{a}}\left[\unicode[STIX]{x1D6FB}^{2}\frac{\unicode[STIX]{x2202}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})}{\unicode[STIX]{x2202}t}\right]=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{0})\unicode[STIX]{x1D6FF}(t-t_{0}),\end{eqnarray}$$

where we use a spherical coordinate system for $\boldsymbol{x}$ and $\boldsymbol{x}_{0}$ . To find the general solution to (2.14), we first define a Fourier transformation

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{G}}(\boldsymbol{x},\unicode[STIX]{x1D714};\boldsymbol{x}_{0},t_{0})=\int _{-\infty }^{\infty }G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})\exp (-\text{i}\unicode[STIX]{x1D714}t)\,\text{d}t, & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }{\hat{G}}(\boldsymbol{x},\unicode[STIX]{x1D714};\boldsymbol{x}_{0},t_{0})\exp (\text{i}\unicode[STIX]{x1D714}t)\,\text{d}\unicode[STIX]{x1D714}. & \displaystyle\end{eqnarray}$$

Using (2.15), (2.14) can now be transformed into a Helmholtz equation

(2.17) $$\begin{eqnarray}-\unicode[STIX]{x1D714}^{2}{\hat{G}}(\boldsymbol{x},\unicode[STIX]{x1D714};\boldsymbol{x}_{0},t_{0})-\left(St^{2}+\frac{\text{i}\unicode[STIX]{x1D714}}{Re_{a}}\right)\unicode[STIX]{x1D6FB}^{2}{\hat{G}}(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{0})\exp (-\text{i}\unicode[STIX]{x1D714}t_{0}),\end{eqnarray}$$

where $\text{i}$ is the imaginary unit. A general solution to this equation can be found by expanding Green’s function into eigenfunctions, resulting in

(2.18) $$\begin{eqnarray}\displaystyle G(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t;r_{0},\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0},t_{0}) & = & \displaystyle \mathop{\sum }_{nlm}\unicode[STIX]{x1D713}_{nl}^{m}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\unicode[STIX]{x1D713}_{nl}^{m}(r_{0},\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0})\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }\left(\frac{\exp (\text{i}\unicode[STIX]{x1D714}(t-t_{0}))}{St^{2}\unicode[STIX]{x1D6FD}_{nl}^{2}+\displaystyle \frac{\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FD}_{nl}^{2}}{Re_{a}}-\unicode[STIX]{x1D714}^{2}}\right)\,\text{d}\unicode[STIX]{x1D714},\end{eqnarray}$$

where $\unicode[STIX]{x1D713}_{nl}^{m}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ are the eigenfunctions of the spherical Helmholtz equation

(2.19) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{nl}^{m}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\frac{\sqrt{2}j_{l}(\unicode[STIX]{x1D6FD}_{nl}r)Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})}{j_{l+1}(\unicode[STIX]{x1D6FD}_{nl})},\end{eqnarray}$$

$j_{l}$ are the spherical Bessel functions, $Y_{l}^{m}$ are the spherical harmonics and $\unicode[STIX]{x1D6FD}_{nl}$ are the zeros of the spherical Bessel functions. To evaluate the inverse Fourier transform (2.16), we use complex contour integration (see figure 4). It can be shown that the contribution of the arc is zero in the limit where the contour radius $a\rightarrow \infty$ . Furthermore, there should be no response of an impulse released at $t_{0}$ at earlier times $t<t_{0}$ (causality condition). To this end, we pick the Jordan curve illustrated in figure 4(a) for $t>t_{0}$ and the curve of figure 4(b) for $t<t_{0}$ , in the limit $a\rightarrow \infty$ . A closed form expression of the Green’s function is now given by

(2.20) $$\begin{eqnarray}G(\boldsymbol{x},t;\boldsymbol{x}_{0},t_{0})=\mathop{\sum }_{nlm}\unicode[STIX]{x1D713}_{nl}^{m}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\unicode[STIX]{x1D713}_{nl}^{m}(r_{0},\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0})\exp (-\unicode[STIX]{x1D705}_{nl}(t-t_{0}))\frac{\sin (\unicode[STIX]{x1D702}_{nl}(t-t_{0}))}{\unicode[STIX]{x1D702}_{nl}}{\mathcal{H}}(t-t_{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{nl}=\unicode[STIX]{x1D6FD}_{nl}^{2}/2Re_{a}$ and $\unicode[STIX]{x1D702}_{nl}=(\sqrt{4St^{2}\unicode[STIX]{x1D6FD}_{nl}^{2}-\unicode[STIX]{x1D6FD}_{nl}^{4}/Re_{a}^{2}})/2$ . The resulting spatio-temporal pressure field using (2.13) reads (without any initial condition)

(2.21) $$\begin{eqnarray}\displaystyle & & \displaystyle p_{1}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=\mathop{\sum }_{nl}\unicode[STIX]{x1D6FD}_{nl}\frac{j_{l}(\unicode[STIX]{x1D6FD}_{nl}r)}{j_{l+1}(\unicode[STIX]{x1D6FD}_{nl})}\left(1-\frac{j_{l-1}(\unicode[STIX]{x1D6FD}_{nl})}{j_{l+1}(\unicode[STIX]{x1D6FD}_{nl})}\right)\frac{2l+1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\int _{0}^{t}\exp (-\unicode[STIX]{x1D705}_{nl}(t-t_{0}))\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{\sin (\unicode[STIX]{x1D702}_{nl}(t-t_{0}))}{\unicode[STIX]{x1D702}_{nl}}{\mathcal{H}}(t-t_{0})P_{l}(\cos \unicode[STIX]{x1D6FE})\left(St^{2}p_{i}(1,\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0},t_{0})+\frac{1}{Re_{a}}\frac{\unicode[STIX]{x2202}p_{i}(1,\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0},t_{0})}{\unicode[STIX]{x2202}t_{0}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\sin (\unicode[STIX]{x1D703}_{0})\,\text{d}\unicode[STIX]{x1D703}_{0}\,\text{d}\unicode[STIX]{x1D719}_{0}\,\text{d}t_{0},\end{eqnarray}$$

where ${\mathcal{H}}$ is the Heaviside theta function, $P_{l}$ are the Legendre polynomials and $\cos (\unicode[STIX]{x1D6FE})=\cos (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D703}_{0})+\sin (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703}_{0})\cos (\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{0})$ . In § 3, we will use (2.21) using a particular pressure boundary condition specified by $p_{i}(1,\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D719}_{0},t_{0})$ . We note that this condition is integrated into the solution at $r=1^{-}$ , i.e.  $r=1-\unicode[STIX]{x1D716}$ where $\unicode[STIX]{x1D716}\rightarrow 0$ .

Figure 4.

Jordan curves used to evaluate the inverse Fourier transform (2.18) in the complex plane where $a\rightarrow \infty$ . (a) Contour used for $t>t_{0}$ , which includes the poles. (b) Contour used for $t<t_{0}$ to obey the causality condition for the Green’s function.

2.2.2 The velocity field

The velocity field inside the droplet is given by (2.10). Since there is no initial rotation present in the fluid and there are no rotational forces acting at later times, the velocity field remains irrotational and is given by a scalar potential. A straightforward time integral over the pressure gradient (the first term on the right-hand side) based on the spherical Bessel functions (2.21) results in a divergent series. To overcome this problem, we solve the velocity field in a different function basis. To this end, we first define the time integral over the thermodynamic pressure as the pressure impulse

(2.22) $$\begin{eqnarray}J_{1}(\boldsymbol{x},t)=\int _{0}^{t}p_{1}(\boldsymbol{x},t^{\prime })\,\text{d}t^{\prime }.\end{eqnarray}$$

The governing equation for the pressure impulse can be obtained by integration of (2.11) in time

(2.23) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}J_{1}=\frac{1}{St^{2}}\frac{\unicode[STIX]{x2202}p_{1}}{\unicode[STIX]{x2202}t}-\frac{1}{St^{2}Re_{a}}\unicode[STIX]{x1D6FB}^{2}p_{1},\end{eqnarray}$$

where we used that both the pressure field $p_{1}$ and its derivative vanish at $t=0$ . It now becomes apparent that the natural basis functions for the pressure impulse are in fact harmonic functions which results in a convergent series.

The general solution for the pressure impulse inside the droplet is therefore given by

(2.24) $$\begin{eqnarray}\displaystyle J_{1}(\boldsymbol{x},t) & = & \displaystyle \iiint _{V}G(\boldsymbol{x};\boldsymbol{x}_{0})\left[\frac{1}{St^{2}}\frac{\unicode[STIX]{x2202}p_{1}(\boldsymbol{x}_{0},t)}{\unicode[STIX]{x2202}t}-\frac{1}{St^{2}Re_{a}}\unicode[STIX]{x1D6FB}^{2}p_{1}(\boldsymbol{x}_{0},t)\right]\,\text{d}V_{0}\nonumber\\ \displaystyle & & \displaystyle -\,\oint _{S}[G(\boldsymbol{x},\boldsymbol{x}_{0})\unicode[STIX]{x1D735}_{\boldsymbol{x}_{0}}J_{1}(\boldsymbol{x}_{0},t)-J_{1}(\boldsymbol{x}_{0},t)\unicode[STIX]{x1D735}_{\boldsymbol{x}_{0}}G(\boldsymbol{x},\boldsymbol{x}_{0})]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S_{0},\end{eqnarray}$$

where the Green’s function satisfies the Poisson equation in spherical coordinates

(2.25) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}G(\boldsymbol{x};\boldsymbol{x}_{0})=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{0}).\end{eqnarray}$$

Completely analogous to the boundary conditions on $p_{1}$ , the boundary condition on $J_{1}$ is $J_{1}(1,t)=0$ , which yields

(2.26) $$\begin{eqnarray}G(\boldsymbol{x},\boldsymbol{x}_{0})=\frac{1}{4\unicode[STIX]{x03C0}}\left(\frac{1}{\sqrt{r^{2}r_{0}^{2}+1-2rr_{0}\cos (\unicode[STIX]{x1D6FE})}}-\frac{1}{\sqrt{r^{2}+r_{0}^{2}-2rr_{0}\cos (\unicode[STIX]{x1D6FE})}}\right),\end{eqnarray}$$

where $\cos (\unicode[STIX]{x1D6FE})=\cos (\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D703}_{0})+\sin (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703}_{0})\cos (\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{0})$ . We evaluate (2.24) numerically using Mathematica 11.1 to speed up the calculations as compared to evaluating (2.24) analytically (Wolfram Research Inc. 2017). The velocity field now reads

(2.27) $$\begin{eqnarray}\boldsymbol{u}_{1}(\boldsymbol{x},t)=-St\unicode[STIX]{x1D735}J_{1}(\boldsymbol{x},t)-\left(\frac{1}{Re}+\frac{1}{Re_{\unicode[STIX]{x1D708}}}\right)\unicode[STIX]{x1D735}p_{1},\end{eqnarray}$$

where $\boldsymbol{u}_{1}(\boldsymbol{r},0)=\mathbf{0}$ and $p_{1}(\boldsymbol{r},0)=0$ .

2.3 The lattice Boltzmann method

To validate our analytic expression for the pressure field inside the droplet, we employ axisymmetric compressible multiphase lattice Boltzmann simulations (Succi Reference Succi2001). Instead of solving the compressible Navier–Stokes equations directly, the numerical method solves a probability distribution function for the position and momentum of particles inside a domain. The evolution of the particle distribution in space and time is given by the Boltzmann equation which is solved numerically on a lattice with a reduced set of velocity vectors. It can be systematically shown that this method is able to correctly solve the fully compressible mass (2.4) and momentum conservation (2.5) equations (Shan, Yuan & Chen Reference Shan, Yuan and Chen2006).

The pressure in our numerical simulation is given by the van der Waals equation of state. This non-ideal equation of state allows us to successfully simulate a liquid–gas system where the interface dynamics is automatically described by the non-ideal pressure. In the vicinity of equilibrium, the van der Waals equation of state behaves as an ideal gas. Therefore, we can directly compare the pressure field inside the droplet with our analytic expression (2.21) in the weakly compressible regime.

Our numerical set-up consists of an initially stationary liquid droplet surrounded by a gas, with a liquid–gas density ratio of $\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}\sim 170$ . We identify the position of the interface by tracking the derivative of the density field where spikes reveal the location of the interface. We hit our droplet by applying a force directly on the liquid–gas interface. By tuning the amplitude (Mach number) and duration (Strouhal number) of this force, we simulate different regimes in the phase space; see figure 3. Further details on this numerical method can be found in Reijers, Gelderblom & Toschi (Reference Reijers, Gelderblom and Toschi2016).

Unfortunately, the compressible multiphase lattice Boltzmann method lacks the advantages of an adaptive grid refinement technique at the moment of writing. Therefore simulations of higher density ratios of $\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}>1000$ , simulations in the strongly compressible regime (see figure 3) or even long-time simulations in the weakly compressible regime are not possible. Only early time dynamics can be simulated with the computational resources available to us, which we use to validate the analytic expression for the pressure in the next section.