# Langmuir instability: wave scattering on circulation flow

S. S. Vergeles1,2, I. A. Vointsev1,2 1L. D. Landau Institute for Theoretical
Physics RAS, Prosp. Akademika Semenova 1A, Chernogolovka, Moscow region,
142432, Russian Federation
2National Research University Higher School of Economics, Faculty of Physics,
Myasnitskaya 20, 101000 Moscow, Russia

(August 28, 2024)

###### Abstract

We consider a classical problem about dynamic instability that leads to the
Langmuir circulation. The problem statement assumes that there is initially a
wind-driven shear flow and a plane surface wave propagating in the direction
of the flow. The unstable mode is a superposition of shear flow and surface
wave both modulated in the horizontal spanwise direction and the circulation
that is vortices in form of near-surface rolls with axes aligned along the
shear streamlines and of transverse size corresponding to the modulation
period. The novelty of our approach is we account for the scattering of the
initial surface wave on the slow current component of the unstable mode that
produces the modulated wave, the interference of the scattered and the initial
waves that produces the modulation of the Stokes drift in the spanwise
direction and the subsequent additional amplification of the circulation by
the vortex force produced by the nonlinear interaction of initial shear flow
and the modulated part of the Stokes drift. Previously, it was shown by S.
Leibovich & A.D.D. Craik that the third part of the mechanism can maintain the
Langmuir circulation. We calculate the growth increment which is larger than
that obtained by A.D.D. Craik. Considering the wave scattering, we describe
fast wave motion as a potential flow with relatively weak vortical correction.
Application of the technique can be expanded on others flows where fast
oscillating surface waves coexist with a slow current.

##  I Introduction

Langmuir circulation (LC) is an important mechanism which is responsible for
enhancement of turbulent mixing in the ocean surface boundary layer, see e.g.
[1] and reviews [2, 3]. The circulation arises on the background of a plane
surface wave propagating along the streamlines of a vertical shear flow,
typically driven by air wind co-directed with the wave. It includes vortical
near-surface rolls with axes oriented in the wave-wind direction and co-
directional shear flow modulated in the horizontal spanwise direction.
Although the most interesting question is the description of established LC
[4, 5], the process of emergence and growth of the circulation is also of
certain interest [6, 7], as it provides a tool for analysis of experimental
[8] and numerical [9] data. This emergence could be described by an
instability mechanism. If one drops out wave scattering on the circulation,
the following positive feedback leads to the instability. Consider a small
spanwise modulation of the shear flow. The flow possesses vertical vorticity
in contrast to the initial shear flow. The wave affects the modulated current,
the impact is described by vortex force [10], that is the cross product of
Stokes drift produced by the wave motion and vorticity of the slow current.
Thus the vorticity component and the Stokes drift produce spanwise vortex
force that drives the roll flow. The latter interacts with the initial shear
flow via the Lamb term in Navier-Stokes equation that drives the spanwise
modulation of the shear flow, so the feedback loop is closed.

Treatment [6] does not consider the wave scattering on the modulated flow,
although it exists and leads to the spanwise modulated wave according to the
further analytical investigations [11]. Subsequently, this modulation was
observed in the large variety of experimental and numerical works. For
example, large-eddy simulations of interacting turbulent shear flow and
surface waves [12] showed rather noticeable variation of wave field in
crosswind direction that can not be omitted when describing oceanic boundary
layer. In [13] the strong LC-like modulation of wave field was proved during
laboratory experiments. The results of recent numerical calculations also
indicate that the scattering of waves by the flow and further interference of
waves cannot be neglected [14]. A general conclusion is that the influence of
established LC on the surface motion is an important factor in the evolution
of the ocean surface boundary layer. Within analytical approach, effect from
the interference of two waves propagating at small angles of equal magnitudes
and opposite signs to the shear flow was considered separately in [10, 15] and
was shown to cause linear growth of LC if it is absent or weak and to maintain
the circulation when the growth is saturated due to the turbulent viscosity.
In this paper we show that the scattering leads to the enhancement of the
positive feedback so the instability develops faster than it is predicted in
[6].

In [11], the wave modulation by LC was considered and the instability problem
was examined when spanwise modulation period is short compared to the
wavelength. In the limit, no wave scattering occurs since the wavenumber of
the modulation exceeds the wavenumber of the plane (fundamental) wave and the
scattered wave would have almost the same frequency because the time scale
associated with vortical flow should be assumed to be much greater compared to
the wave period. Here we consider the opposite case of wavelength short
compared to the modulation period, so the fundamental wave does scatters on
the circulation and the scattered wave should be taken into account. The
stream function used in [11] is convenient when there is a single wave vector
characterising the wave motion. As we consider essentially three-dimensional
wave flow which consists of waves spreading in different directions, we use
potential approximation for the flow. On the way, the relatively weak vortical
correction to the wave flow is described in terms of vorticity [10]. This
technical solution allows us to consistently take into account wave scattering
on a vortex flow. The only requirement for the applicability of the developed
general mathematical apparatus is that the wave frequency has to be much
greater than the velocity gradient and the rate of change of the vortex flow.
The time scale separation allows one to consider the wave motion and the
vortical flow as independent subsystems in zero approximation which weakly
interact due to hydrodynamic nonlinearity. In the problem of Langmuir
instability, we show that the spanwise modulation of the wave is essential at
any strength of the shear flow. In particular, we obtain the increment for an
unstable mode larger than that found in [6]. The established analytical
results are applicable if the shear rate is much less than the wave frequency.

The general mathematical scheme is developed in Section II. Derivation of
vortex force is reproduced in Subsection II.1 and wave scattering process is
considered in Subsection II.2. Using this technique, the Langmuir instability
mechanism is considered in Section III. The unperturbed flow and the
perturbation structure are defined in Subsections III.1, III.2. Equations
governing the linearized development of the perturbation are derived in
Subsection III.3. Then we consider the limit of small Langmuir number
$\mathrm{La}\ll 1$ and present complete analytical solution for the problem in
the case of linear shear in Subsection III.4. Some calculations are swept into
Appendix.

##  II General description of interacting waves and current

We consider an incompressible flow with a free surface, along which surface
gravity waves can propagate. We consider deep-water case with equation
$z=\zeta(x,y,t)$ describing surface shape. Velocity field is a sum of surface
gravity wave term $\bm{u}$ and slow current $\bm{V}$ excited against the
background, i.e. $\bm{v}=\bm{u}+\bm{V}$. The characteristic time scale $T$ of
slow current, $T\sim|\partial_{t}{\bm{V}}|/|{\bm{V}}|$ or $T\sim
1/|\mathop{\mathrm{grad}}{\bm{V}}|$ is much larger than the inverse gravity
wave frequency $\omega$, $\varepsilon=1/\omega T\ll 1$, so we call ${\bm{u}}$
high-frequency component of the whole flow. In particular, this means that the
surface elevation is determined by the wave motion only, so the Froude number
of the slow current is considered to be small,
$\operatorname{Fr}\sim|{\bm{V}}|\cdot|\mathop{\mathrm{grad}}{\bm{V}}|/g\ll 1$,
where $g$ is the gravitational acceleration. Note that the spatial scale of
the current $L\sim|{\bm{V}}|/|\mathop{|\mathrm{grad}}{\bm{V}}|$ is not assumed
to be necessarily large in relation to the wavelength, but the smallness of
parameter $\varepsilon$ is an essential condition for our analytical approach
to be applicable.

The fluid motion is described by Navier-Stokes equation which can be written
in the Lamb form [16]

$\partial_{t}{\bm{v}}=[\bm{v},\,\bm{\varpi}]-\operatorname{grad}\left(p+\frac{v^{2}}{2}\right)+\nu\Delta\bm{v}+\bm{g},$
(1)

where $p$ is fluid pressure, ${\bm{\varpi}}=\mathop{\mathrm{curl}}{\bm{v}}$ is
vorticity and mass density is equal to one. The equation should be
supplemented with boundary conditions that are kinematic boundary condition

$\partial_{t}\zeta=v_{z}-v_{\alpha}\partial_{\alpha}\zeta,$ (2)

where indices $\alpha,\beta,\ldots$ runs values $\\{x,y\\}$ and dynamic
boundary condition

$\displaystyle\big{(}p-2\nu
n_{i}n_{k}\partial_{i}v_{k}\big{)}\big{|}_{z=h}={}$ $\displaystyle 0,$ (3)
$\displaystyle\delta^{\scriptscriptstyle\perp}_{li}\,(\partial_{i}v_{k}+\partial_{k}v_{i})n_{k}\big{|}_{z=h}={}$
$\displaystyle 0,$ (4)

where $\delta^{\scriptscriptstyle\perp}_{li}=\delta_{li}-n_{l}n_{i}$ is the
projector onto the plane tangent to the free surface, ${\bm{n}}$ is normal
unit vector on the surface and indices $i,j,k,\ldots$ runs values
$\\{x,y,z\\}$. Hereinafter, summation is assumed over the repeating indices.

In the absence of the current and the viscosity, the wave non-breaking flow
${\bm{u}}$ is irrotational so can be described in term of flow potential
$\phi$. The interaction of the wave with the vortical current leads to the
fast component of the flow oscillating with the wave frequency ceases to be
pure irrotational. There is an analytical description of the phenomenon based
on the wave flow representation through stream function instead of the
potential $\phi$, which is good for two-dimensional problems [17]. When the
wave motion is substantially three-dimensional, the technique becomes
cumbersome since it is impossible to describe the wave flow with only one
scalar stream function so one should introduce more functions [11]. Our aim is
to describe the dynamics of the wave motion itself including its interaction
with the current. The result of the interaction at large distances are
scattered waves which are reasonable to be described in terms of the potential
$\phi$ as well. There is a local response on the interaction, which is
relatively small as $\varepsilon$ and can be described in terms of oscillating
vorticity $\bm{\varpi}^{u}=\operatorname{curl}{\bm{u}}$ and pressure part
$p^{\varpi}$ (see below). Thus, in our approach we keep the potential $\phi$
which determines the wave motion itself and the vortical correction
$\bm{u}^{\varpi}$ so that

$\bm{u}=\bm{u}^{\phi}+\bm{u}^{\varpi},\quad\bm{u}^{\phi}=\operatorname{grad}\phi,\qquad
u^{\varpi}\sim\varepsilon\cdot u^{\phi}.$ (5)

For the division to be completely defined, we assume that $\bm{u}^{\varpi}$ is
not related to the surface dynamics in the linear approximation

$u^{\varpi}_{z}|_{z=0}=0,$ (6)

since the surface shape dynamics is closely related to the potential wave
motion.

###  II.1 Influence of the wave motion on the slow current

The impact of the waves on the current is established [10, 15] and is called
vortex force. On the way of its derivation one should find the high-frequency
part ${\bm{\varpi}}^{u}$ of the vorticity, that is also used when finding the
scattered waves as well. So here we briefly repeat the derivation of the
vortex force $\bm{f}^{V}$.

The projection of Navier-Stokes equation (1) onto the slow motion is

$\partial_{t}\bm{V}+(\bm{V}\,\nabla)\bm{V}=-\operatorname{grad}\left\langle
p+\frac{u^{2}}{2}\right\rangle+\nu_{\scriptscriptstyle
T}\Delta{\bm{V}}+\langle[\bm{u}^{\phi},\,\bm{\varpi}^{u}]\rangle.$ (7)

Here we have changed molecular viscosity by turbulent viscosity
$\nu_{\scriptscriptstyle T}$ assuming that there is some small-scale turbulent
flow on the background of flow ${\bm{V}}$ and the wave motion ${\bm{u}}$. The
last term in (7) stems from the high-frequency part of the flow, which is pure
potential in the main approximation, so it was taken ${\bm{u}}^{\phi}$ instead
of ${\bm{u}}$. To find the high-frequency part ${\bm{\varpi}}^{u}$ of the
vorticity, we linearize the vorticity equation

$\partial_{t}{\bm{\varpi}}=\mathop{\mathrm{curl}}[{\bm{v}},{\bm{\varpi}}]+\nu\Delta{\bm{\varpi}}$
(8)

with respect to the wave motion, that leads to

$\big{(}\partial_{t}+(\bm{V}\,\nabla)\big{)}\bm{\varpi}^{u}-(\bm{\varpi}^{u}\nabla)\bm{V}-\nu\Delta\bm{\varpi}^{u}=\operatorname{curl}\left[\bm{u}^{\phi},\,\bm{\Omega}\right],$
(9)

where the vorticity of the slow flow
${\bm{\Omega}}=\operatorname{curl}{\bm{V}}$, so the full vorticity
${\bm{\varpi}}={\bm{\Omega}}+{\bm{\varpi}}^{u}$. The second term in left-hand
side of equation (9) can be estimated as
$(\bm{\varpi}^{u}\nabla)\bm{V}\sim\Omega\varpi^{u}$ so it should be neglected
compared to $\partial_{t}\bm{\varpi}^{u}$. The right-hand side serves as
driving force for ${\bm{\varpi}}^{u}$. We assume that the spatial scale of the
right-hand side of the equation is much greater than the thickness of the
viscous sublayer. Thus, the influence of the viscosity is negligible and the
equation becomes local in space. Its solution is

$\bm{\varpi}^{u}=\operatorname{curl}\left[\bm{s},\,{\bm{\Omega}}\right],\quad\bm{s}(t,{\bm{r}})=\int\limits^{t}dt^{\prime}\,\bm{u}^{\phi}(t^{\prime},\,\bm{r}(t^{\prime})),$
(10)

where $\dot{\bm{r}}(t)={\bm{V}}(t,{\bm{r}})$ is the Lagrangian trajectory
produced by the slow flow, ${\bm{r}}(t)\equiv{\bm{r}}$ and ${\bm{s}}$ is a
particle displacement during the wave oscillations. Using (10), one can
rewrite the last term in (7) into more convenient form

$\langle[\bm{u}^{\phi},\,\bm{\varpi}^{u}]\rangle=\bm{f}^{V}+\frac{1}{2}\operatorname{grad}(\bm{\Omega}\cdot\bm{\mathcal{A}}),\qquad\bm{f}^{V}=[\bm{U}^{s},\,\bm{\Omega}],$
(11)

where $\bm{\mathcal{A}}=\langle[{\bm{u}}^{\phi},{\bm{s}}]\rangle/2$ and the
Stokes drift

$\bm{U}^{s}=\operatorname{curl}\bm{\mathcal{A}}=\langle\operatorname{curl}\,[{\bm{u}}^{\phi},{\bm{s}}]\rangle$
(12)

is determined by the potential part of the flow associated with the wave
motion. The gradient term in (11) should be included in effective pressure in
(7). Note that all time averaging $\langle\ldots\rangle$ in (7-12) should be
implemented along the Lagrangian trajectories ${\bm{r}}(t)$ defined after
(10).

Finally, the equation (7) takes the form

$\partial_{t}\bm{V}+(\bm{V}\,\nabla)\bm{V}=-\operatorname{grad}\overline{P}+\nu\Delta{\bm{V}}+\bm{f}^{V},$
(13)

where effective time-averaged pressure $\overline{P}$ contains pressure heads
from (7) and (11). Concerning the boundary conditions, we neglect the virtual
wave stress produced by the wave viscous damping [18, 19, 20], that is
reasonable if the vortex flow is strong enough so the velocity gradient much
exceeds the viscous damping rate of the wave, $|\nabla{\bm{V}}|\gg\nu k^{2}$.
Then the boundary conditions for the slow flow are stress-free rigid boundary,

$\partial_{z}V_{\alpha}\big{|}_{z=0}=0,\quad V_{z}\big{|}_{z=0}=0.$ (14)

###  II.2 Wave scattering process

The wave flow dynamics is determined by the Navier-Stokes equation linearized
in the wave amplitude

$\partial_{t}\bm{u}+\left(\bm{u}\nabla\right)\bm{V}+\left(\bm{V}\nabla\right)\bm{u}=-\nabla
p^{u}+\bm{g},$ (15)

where $p^{u}$ is high-frequency part of the pressure. We omitted the viscous
term in (15) so we deal with Euler equation due to the assumption that the
viscosity which alters the flow in the narrow viscous sublayer beneath the
surface is not essential in the wave dynamics. Our goal is to describe the
dynamics of the potential part of the wave flow (5) taking into account the
scattering process. Due to the incompressibility condition, the potential
still satisfies the Laplace equation

$\Delta\phi=0.$ (16)

Now we impose the boundary conditions. Kinematic boundary condition (2)
linearized in the wave amplitude is

$\left(\partial_{t}+V_{\alpha}\big{|}_{z=0}\partial_{\alpha}\right)\zeta=\partial_{z}\phi+\zeta\partial_{z}V_{z}\big{|}_{z=0},$
(17)

where we used condition (6) for $u^{\varpi}_{z}$. Because the viscosity was
neglected in (15), we need only one dynamical boundary condition for pressure
(3). To obtain this condition one should express the fast oscillating term in
pressure in terms of the potential $\phi$ and surface elevation $\zeta$. In
the case of pure wave flow, this is produced with the aim of Bernoulli
equation, which should be linearized in the wave amplitude in our
approximations, so we have $\partial_{t}p^{u}=-gz-\partial_{t}\phi$. We
generalize the equation in the presence of the slow current, that includes
replacement of the time derivative $\partial_{t}$ by substantial derivative
$(\partial_{t}+V_{i}\partial_{i})$. We denote by $p^{\varpi}$ the residual
part of the pressure, so by definition

$p^{u}\equiv-gz-\left(\partial_{t}+V_{i}\partial_{i}\right)\phi+p^{\varpi}.$
(18)

Thus, the dynamic boundary condition is

$\left(\partial_{t}+V_{\alpha}\big{|}_{z=0}\partial_{\alpha}\right)\phi\big{|}_{z=0}=-g\zeta+\left.p^{\varpi}\right|_{z=0},$
(19)

where the pressure part $p^{\varpi}$ satisfies

$\displaystyle\Delta p^{\varpi}=\partial_{i}\phi\,\Delta
V_{i}+2\partial_{j}u^{\varpi}_{i}\,\partial_{i}V_{j}$ (20)
$\displaystyle\partial_{z}p^{\varpi}\big{|}_{z=0}=\partial_{\alpha}\phi\partial_{z}V_{\alpha},\quad
p^{\varpi}\big{|}_{z\rightarrow-\infty}\rightarrow 0.$ (21)

In (20), the vortical part of the high-frequency flow ${\bm{u}}^{\varpi}$
should be restored using equation
${\bm{\varpi}}^{u}=\mathop{\mathrm{curl}}{\bm{u}}^{\varpi}$ from (10) with
boundary conditions (6) at the free surface and ${\bm{u}}^{\varpi}\to 0$ at
infinity. Note that if velocity ${\bm{V}}$ is pure potential, then
$\Delta{\bm{V}}=0$ and ${\bm{u}}^{\varpi}=0$, so $p^{\varpi}=0$ as well. The
contribution $p^{\varpi}$ into the pressure is related to the vortical part
${\bm{u}}^{\varpi}$ of the wave motion, but the relation is not linear in
${\bm{V}}$.

According to (10), the ratio of the terms in right-hand side of (20) is
estimated as

$\frac{2\partial_{j}u^{\varpi}_{i}\,\partial_{i}V_{j}}{\partial_{i}\phi\,\Delta
V_{i}}\sim\frac{V}{\omega/k}.$ (22)

Thus, the second term in (20) should be taken into account only if the slow
flow is large and the magnitude of its alternation in the horizontal plane is
comparable with the phase velocity $\omega/k$ of the waves. The limit can be
analysed with the ray approximation [21] as well.

Figure 1: Sketch for Langmuir instability development. It was chosen
$\theta=0.1\,\text{rad}$ and $\delta h/h^{\scriptscriptstyle(0)}=0.3$. In the
$Oyz$-slice, the red and blue areas correspond to positive and negative
$\delta\Omega_{x}$ respectively. Plane $y=0$ corresponds to the vertical
middle plane of the depicted cell.

##  III Langmuir instability mechanism

In this Section we consider the Langmuir circulation instability problem. We
assume that the unperturbed flow is a plane monochromatic wave propagating on
the background of co-directed vertically sheared flow. Thus, the unperturbed
flow is uniform in the spanwise direction, and the perturbation is modulated
in the direction with some period, which is much larger that the wavelength. A
sketch of the flow is depicted in Figure 1.

###  III.1 Unperturbed flow

We assume that the wave travels in the $x$ direction in Cartesian coordinate
system. Thus, the wave potential

$\phi^{{\scriptscriptstyle{(0)}}}=\real(\psi^{{\scriptscriptstyle{(0)}}}\exp[i\varphi^{{\scriptscriptstyle{(0)}}}+kz]\Big{)}$
(23)

with $\psi^{{\scriptscriptstyle{(0)}}}$ and $k$ being potential amplitude and
wave number respectively and the phase
$\varphi^{{\scriptscriptstyle{(0)}}}=kx-\omega t$. The corresponding surface
elevation rate is

$\zeta^{{\scriptscriptstyle{(0)}}}=\real(h^{{\scriptscriptstyle{(0)}}}\exp[i\varphi^{{\scriptscriptstyle{(0)}}}]\Big{)}.$
(24)

Without loss of generality we assume $h^{{\scriptscriptstyle{(0)}}}>0$. The
slow current is shear flow aligned in the same $x$ direction

$\displaystyle\bm{V}^{{\scriptscriptstyle{(0)}}}=\big{\\{}V^{{\scriptscriptstyle{(0)}}}_{x},\,0,\,0\big{\\}},\qquad\bm{\Omega}^{{\scriptscriptstyle{(0)}}}=\big{\\{}0,\,\Omega^{{\scriptscriptstyle{(0)}}},\,0\big{\\}},$
(25)

where
$\Omega^{{\scriptscriptstyle{(0)}}}(z)=\partial_{z}V^{{\scriptscriptstyle{(0)}}}_{x}(z)$.
The relative slowness of the flow means that
$\varepsilon\sim\Omega^{{\scriptscriptstyle{(0)}}}/\omega\ll 1$. One can think
about the resulting motion as generalized Gouyon wave, see e.g. [22]. Here we
apply our general scheme to find the influence of the shear onto the wave
demonstrating the reproduction of the well-known result [17].

We start from the Stokes drift, which $x$ component is

$U^{s{\scriptscriptstyle{(0)}}}_{x}=k\omega\big{(}h^{{\scriptscriptstyle{(0)}}}\big{)}^{2}e^{2kz}.$
(26)

The vortex force (11) appears to be pure potential

$\bm{f}^{V}=k\omega\Omega^{{\scriptscriptstyle{(0)}}}\big{(}h^{{\scriptscriptstyle{(0)}}}\big{)}^{2}e^{2kz}\big{\\{}0,\,0,\,1\big{\\}}$
(27)

and, according to the equation (7), does not produce any contribution to the
slow flow but only modifies the pressure.

Dispersion relation and, thus, phase velocity could be obtained from boundary
conditions (17) and (19) on the free surface:

$\displaystyle\big{(}-i\omega+ikV^{{\scriptscriptstyle{(0)}}}_{x}(0)\big{)}h^{{\scriptscriptstyle{(0)}}}=k\psi^{{\scriptscriptstyle{(0)}}},$
(28)
$\displaystyle\big{(}-i\omega+ikV^{{\scriptscriptstyle{(0)}}}_{x}(0)\big{)}\psi^{{\scriptscriptstyle{(0)}}}=-gh^{{\scriptscriptstyle{(0)}}}+p^{\varpi{\scriptscriptstyle{(0)}}},$
(29)

with pressure (note that $\bm{u}^{\varpi\scriptscriptstyle(0)}=0$ for the
geometry)

$p^{\varpi{\scriptscriptstyle{(0)}}}=2ik\psi^{{\scriptscriptstyle{(0)}}}\Big{(}V^{{\scriptscriptstyle{(0)}}}_{x}(0)-2k\int\limits_{-\infty}^{0}dz\,e^{2kz}V^{{\scriptscriptstyle{(0)}}}_{x}(z)\Big{)}.$
(30)

The solution of (28)-(29) with linear precision to
$\Omega^{{\scriptscriptstyle{(0)}}}/\omega\ll 1$ term is

$\omega=\sqrt{gk}+\omega^{\prime},\qquad\omega^{\prime}=2k^{2}\int\limits_{-\infty}^{0}dz\,e^{2kz}V^{{\scriptscriptstyle{(0)}}}_{x}(z),$
(31)

and the phase velocity $\omega/k$ appears the same as in [17]. One can choose
the reference system where $V_{x}^{\scriptscriptstyle(0)}|_{z=0}=0$. We assume
that the shear flow is produced by surface stress co-directed with $Ox$-axis,
as the Langmuir instability exists only in this case [23]. Then
$\Omega^{\scriptscriptstyle(0)}>0$ if one considers the simplest case of
constant sign shear and the correction to the wave frequency
$\omega^{\prime}<0$. If we suppose the linear shear,
${V}^{{\scriptscriptstyle{(0)}}}_{x}=\Omega^{{\scriptscriptstyle{(0)}}}z$,
then the correction to the wave frequency
$\omega^{\prime}=-\Omega^{{\scriptscriptstyle{(0)}}}/2$.

###  III.2 Perturbation

It is assumed in CL2 mechanism [6] of Langmuir instability that an unstable
mode includes a streamwise shear flow and a circulation in vertical-spanwise
plane both modulated in the spanwise direction. Here we add into the mode
structure the wave modulation, that can be thought as the result of scattering
of initial wave on the modulated slow flow. We call the modulated wave oblique
one after [24]. Hereinafter, we mark the contributions into the physical
quantities describing the unstable mode with symbol ‘$\delta$’, so the full
fields are now

$\displaystyle\bm{V}=\bm{V}^{{\scriptscriptstyle{(0)}}}+\delta\bm{V},\quad$
$\displaystyle\phi=\phi^{{\scriptscriptstyle{(0)}}}+\delta\phi,$ (32)
$\displaystyle\bm{\Omega}=\bm{\Omega}^{{\scriptscriptstyle{(0)}}}+\delta\bm{\Omega},\quad$
$\displaystyle\zeta=\zeta^{{\scriptscriptstyle{(0)}}}+\delta\zeta.$

The wave potential $\delta\phi$ and the surface elevation $\delta\zeta$ are:

$\displaystyle\delta\phi$ $\displaystyle=$
$\displaystyle\real(\delta\psi\,e^{\lambda
t}\exp[i\varphi^{\scriptscriptstyle(1)}+kz\big{]}\Big{)}$ (33)
$\displaystyle\delta\zeta$ $\displaystyle=$ $\displaystyle\real(\delta
h\,e^{\lambda t}\exp[i\varphi^{\scriptscriptstyle(1)}\big{]}\Big{)},$ (34)

where phase $\varphi^{\scriptscriptstyle(1)}=kx\cos\theta+ky\sin\theta-\omega
t$ contains spanwise modulation with wave number $k\sin\theta$. The
instability growth rate $\lambda$ is assumed to be much smaller compared to
the wave frequency, $\lambda\ll\omega$, so the wave is assumed to be quasi-
monochromatic. The correction to the oscillating part of the pressure could be
represented in the same form (see Appendix A)

$\delta{p}^{\varpi}=\real(\delta\tilde{p}^{\varpi}(z)\ e^{\lambda
t}\exp[i\varphi^{\scriptscriptstyle(1)}\big{]}\Big{)}.$ (35)

The dynamics and the spatial structure of the slow flow matches the dynamics
of the interference between the initial and scattered waves averaged over the
fast oscillations. Thus, the time dependence is exponential, $\exp(\lambda
t)$, and the spatial structure of the slow flow is determined by the phase
difference
$\varphi^{\scriptscriptstyle(1)}-\varphi^{\scriptscriptstyle(0)}=ky\sin\theta+kx(\cos\theta-1)$:

$\displaystyle\delta\bm{V}=\real(\delta\tilde{\bm{V}}(z)\,e^{\lambda
t}\exp[i\big{(}\varphi^{\scriptscriptstyle(1)}-\varphi^{\scriptscriptstyle(0)}\big{)}\big{]}\Big{)}.$
(36)

The same notations are considered for $\delta{\bm{\Omega}}$.

###  III.3 Getting equations

The next step is to derive the equations which define the dynamics of the
perturbation in the approximation linear in its amplitude. The full system of
equations to be linearized in the perturbation amplitude on the background of
the flow described in Section III.1 consists of equations (16,17,19,21) for
the fast oscillating motion and equations (13,14) for the slow flow. We note
that one should keep only the first summand in the right-hand side of Eq. (20)
since the amplitude of the perturbation is small. Note, that in the case we do
not need to restore ${\bm{u}}^{\varpi}$ from ${\bm{\varpi}}^{u}$. The slow
flow part of the perturbation is completely defined by $x$-components of
velocity and vorticity $\delta\tilde{V}_{x}(z)$, $\delta\tilde{\Omega}_{x}(z)$
and its fast oscillating part is set by the oblique wave parameters
$\delta\psi$, $\delta h$, so we establish equations for these quantities. We
examine the problem in the limit $\theta\ll 1$, so we neglect all corrections
which are relatively small as $\theta^{2}$. In particular, this means that
$\cos\theta=1$ in our approximation and $x$-derivatives of slow flow variables
are zero as the derivatives are proportional to $1-\cos\theta$. This
approximation corresponds to $x$-independence of the slow flow perturbation
adopted in [6] and allows one to introduce the stream function $\delta{\Psi}$
for the circulation in $Oyz$-plane:

$\displaystyle\delta{V}_{y}=\partial_{z}\delta{\Psi},\quad\delta{V}_{\\!z\,}=-\partial_{y}\delta{\Psi}=-ik\theta\delta{\Psi},$
$\displaystyle\delta{\Omega}_{x}=-(\partial_{z}^{2}-k^{2}\theta^{2})\delta{\Psi}.$
(37)

Getting the final answer, we choose the particular case of linear shear flow,
${V}^{{\scriptscriptstyle{(0)}}}_{x}=\Omega^{{\scriptscriptstyle{(0)}}}z$,
consider the limit of small Langmuir number $\mathrm{La}\ll\theta^{2}$, see
(59), and neglect the corrections which are relatively small as $\theta$ as
well. It is convenient to divide the solution into two steps. The first step
is determining the slow current part via the oblique wave parameters. The step
can be implemented analytically in the case of linear unperturbed shear
${V}^{{\scriptscriptstyle{(0)}}}_{x}=\Omega^{{\scriptscriptstyle{(0)}}}z$. On
the second step, we define the instability growth rate $\lambda$ from the
requirement that free surface equations (17) and (19) have nontrivial
solution. We also adopt below that all quantities describing the unstable mode
are complex, that is, in particular, one should drop $\mathrm{Re}$ in
definitions (33,34,35,36) for the perturbation and (23,24) for the initial
wave.

Let us first linearize in the perturbation amplitude equation (13) and
equation (8) written for the slow flow:

$\displaystyle\big{(}\lambda-\nu_{\scriptscriptstyle
T}\Delta+V^{{\scriptscriptstyle{(0)}}}_{x}\partial_{x}\big{)}\delta\bm{V}+\delta{V}_{\\!z\,}\bm{\Omega}^{{\scriptscriptstyle{(0)}}}+\operatorname{grad}\delta\overline{P}=\delta\bm{f}^{\scriptscriptstyle
V},$ (38) $\displaystyle\big{(}\lambda-\nu_{\scriptscriptstyle
T}\Delta+V^{{\scriptscriptstyle{(0)}}}_{x}\partial_{x}\big{)}\delta\bm{\Omega}-\Omega^{\scriptscriptstyle(0)}\big{[}\delta\Omega_{z}\bm{e}^{x}+\partial_{y}\delta\bm{V}\big{]}=$
$\displaystyle=\operatorname{curl}\delta\bm{f}^{\scriptscriptstyle V},$ (39)

where we took into consideration the variation of vortex force

$\delta\bm{f}^{\scriptscriptstyle
V}=\big{[}\bm{U}^{s{\scriptscriptstyle{(0)}}},\,\delta\bm{\Omega}\big{]}+\big{[}\delta\bm{U}^{s},\,\bm{\Omega}^{{\scriptscriptstyle{(0)}}}\big{]}$
(40)

as well. The correction to the Stokes drift velocity defined in (12) is

$\delta\bm{U}^{s}=\frac{1}{2}\,\partial_{l}\big{(}\delta
s_{l}^{*}\bm{u}^{\phi{\scriptscriptstyle{(0)}}}\big{)}+\frac{1}{2}\,\partial_{l}\big{(}s^{{\scriptscriptstyle{(0)}}*}_{l}\delta\bm{u}^{\phi}\big{)},$
(41)

where symbol ⋆ denotes complex conjugation. In the limit of small angle
$\theta$, the Stokes drift perturbation is

$\displaystyle\delta U^{s}_{x}={}$ $\displaystyle
2ik^{2}h^{{\scriptscriptstyle{(0)}}}\delta\psi\,e^{\lambda
t+2kz}\exp[i(\varphi^{\scriptscriptstyle(1)}-\varphi^{{\scriptscriptstyle{(0)}}}))\big{]},$
(42)

and $\delta U^{s}_{y}=\theta\,\delta U^{s}_{x}$, but $y$-component is not
interesting because it does not produce any contribution into vortex force
(40). The perturbation (40) of the vortex force appears to be directed in
$yz$-plane. Now we take $x$-components in equations (38,III.3) and obtain in
our approximations

$\displaystyle\big{[}\lambda-\nu_{\scriptscriptstyle
T}\big{(}\partial_{z}^{2}-k^{2}\theta^{2}\big{)}\big{]}\delta\tilde{V}_{x}+\Omega^{{\scriptscriptstyle{(0)}}}\delta\tilde{V}_{z}=0.$
(43) $\displaystyle\big{[}\lambda-\nu_{\scriptscriptstyle
T}\big{(}\partial_{z}^{2}-k^{2}\theta^{2}\big{)}\big{]}\delta\tilde{\Omega}_{x}+2ik^{2}\theta
U^{s{\scriptscriptstyle{(0)}}}\delta\tilde{V}_{x}=$
$\displaystyle=ik\theta\Omega^{\scriptscriptstyle(0)}\delta\tilde{U}^{s}_{x}.$
(44)

Note that $x$-component of square bracket in (III.3) was neglected as it is
equal to $\partial_{x}\delta V_{y}$. Equations should be supplemented with
boundary conditions on the free surface (14)

$\partial_{z}\delta V_{x}\big{|}_{z=0}=\partial_{z}\delta
V_{y}\big{|}_{z=0}=0,\quad\delta{V}_{\\!z\,}\big{|}_{z=0}=0.$ (45)

Next we rewrite equations (17,19),

$\displaystyle(\lambda-i\omega)\delta{h}={}$ $\displaystyle
k\delta\psi+h^{\scriptscriptstyle(0)}\partial_{z}\delta{V}_{\\!z\,}\big{|}_{z=0},$
(46) $\displaystyle(\lambda-i\omega)\delta\psi={}$
$\displaystyle-g\delta{h}+\delta{p}^{\varpi}\big{|}_{z=0}.$ (47)

Taking into account the boundary conditions for the slow current (45), the
pressure $\delta\tilde{p}^{\varpi}$ taken at the surface can be presented in
the form (see (67))

$\displaystyle\delta\tilde{p}^{\varpi}\big{|}_{z=0}$ $\displaystyle=$
$\displaystyle-\psi^{\scriptscriptstyle(0)}\int\limits_{-\infty}^{0}dz\,e^{2kz}(\partial_{z}^{2}-k^{2}\theta^{2})\big{(}i\delta\tilde{V}_{x}+\delta\tilde{V}_{z}\big{)}+$
(48) $\displaystyle+i\Omega^{\scriptscriptstyle(0)}\big{|}_{z=0}\delta{\psi}.$

###  III.4 Solution

Now we assume, that the perturbations growth rate $\lambda$ is large compared
to viscous decay rate of the slow flow component of the perturbation. Since
the Stokes drift in (43) penetrates on depth $\sim 1/k$, the vertical scale of
$\delta\Omega_{x}$ is the same. This means, that the viscous damping rate is

$\nu_{\scriptscriptstyle T}k^{2}\ll\lambda.$ (49)

Below we show that the inequality is equivalent to the limit of small Langmuir
number. In the limit, the viscous scale
$\delta_{\nu}\sim\sqrt{\nu_{\scriptscriptstyle T}/\lambda}$ for the slow flow
is less than the wave penetration depth $1/k$, $k\delta_{\nu}\ll 1$ and
$\lambda$ prevails the viscous terms in square brackets in (43,III.3). We
neglect the viscous term in equation (43) and obtain

$\delta\tilde{V}_{x}=-\frac{\Omega^{{\scriptscriptstyle{(0)}}}}{\lambda}\delta\tilde{V}_{z}=\frac{ik\theta\Omega^{{\scriptscriptstyle{(0)}}}}{\lambda}\delta\tilde{\Psi}$
(50)

due to definition (III.3). According to (50), the mutual signs of
$\delta\tilde{V}_{x}$ and $\delta\tilde{V}_{z}$ are fixed due to implied
$\Omega^{{\scriptscriptstyle{(0)}}}>0$ and, as it will be shown below, the
increment is real, $\lambda>0$. Downward flow in the circulation corresponds
to maximum in the amplitude of the shear flow. On the other hand, the downward
flow corresponds to divergent stream lines on the surface, i.e. areas where a
surface contamination is collected along $x$-oriented lines. Under the same
approximation, equation (III.3) gives

$\displaystyle(k^{-2}\partial_{z}^{2}-\theta^{2})\delta\tilde{\Psi}+\mu^{2}\,e^{2kz}\delta\tilde{\Psi}=\sqrt{\varepsilon}\,\mu\,e^{2kz}\delta{\psi},$
(51)

where the dimensionless quantities $\mu$ and $\varepsilon$ are

$\mu=\dfrac{\sqrt{2\Omega^{\scriptscriptstyle{(0)}}\omega}}{\lambda}\,kh^{\scriptscriptstyle{(0)}}\theta,\quad\varepsilon=\dfrac{2\Omega^{\scriptscriptstyle{(0)}}}{\omega^{\,}}\ll
1.$ (52)

As we are considering inviscid problem, the boundary conditions for the stream
function are

$\delta\tilde{\Psi}(0)=\delta\tilde{\Psi}(-\infty)=0$ (53)

instead of (45). Note that we gain the CL2-model equation [6] if we neglect
the variation of Stokes drift, i.e. let the right-hand side of (51) being
equal to zero.

Further we consider linear shear flow, so parameters $\mu$ and $\varepsilon$
are $z$-independent after $\Omega^{\scriptscriptstyle{(0)}}$. In the absence
of Stokes drift variation in (51), the eigenvalue problem (51,53) leads to the
solution $\delta\tilde{\Psi}={\mathrm{J}}_{\theta}(\mu\,e^{kz})$ and boundary
condition ${\mathrm{J}}_{\theta}(\mu)=0$ [6], where $\mathrm{J}_{n}(\mu)$ is
Bessel function of $n$-th order. Thus the perturbation grows rate in CL2-model
corresponds to the lowest root $\mu_{\text{\tiny{CL2}}}\approx 2.4$ in the
limit $\theta\ll 1$. Now we restore the right-hand side of (53). The equation
(51) can be solved analytically for the constant initial shear flow:

$\displaystyle\delta\tilde{\Psi}=\dfrac{\sqrt{\varepsilon}}{\mu}\Bigg{[}1-\Big{(}1+\dfrac{1-{\mathrm{J}}_{0}(\mu)}{{\mathrm{J}}_{0}(\mu)}\,e^{|\theta|kz}\Big{)}{\mathrm{J}}_{0}\big{(}\mu\,e^{kz}\big{)}\Bigg{]}\delta{\psi}.$
(54)

Expression (54) is not exact, there are corrections for the exact solution
which are relatively small as $\theta$. With the accuracy, one should put the
exponent $e^{|\theta|kz}\to 1$ when $kz\sim 1$ so the square bracket in (54)
is equal to
$1-{\mathrm{J}}_{0}\big{(}\mu\,e^{kz}\big{)}/{\mathrm{J}}_{0}(\mu)$. At large
depths $z\gtrsim 1/k\theta$, the stream function decays exponentially as the
square bracket in (54) is $(1-1/{\mathrm{J}}_{0}(\mu))\,e^{|\theta|kz}$. Thus,
the circulation penetrates $1/k\theta$ deep. However, the vorticity related to
the stream function according to (III.3) penetrates only $1/k$ deep,

$\delta\tilde{\Omega}_{x}=i\sqrt{\varepsilon}\,\dfrac{\mu{\mathrm{J}}_{0}(\mu\,e^{kz})}{{\mathrm{J}}_{0}(\mu)}\,e^{2kz}\omega
k\delta{h},$ (55)

where we used approximate equality $\delta\psi=-i(\omega/k)\delta h$. Note
that in CL2 model, the vorticity has the same vertical dependence in our
approximations, $\delta\tilde{\Omega}_{x}=(\mu
k)^{2}e^{2kz}{\mathrm{J}}_{0}(\mu e^{kz})$.

Equations (50,54) express the slow flow via the oblique wave parameters
$\delta h$, $\delta\psi$. To be able to write explicitly equations (46,46) in
terms of the parameters, we need to do this first for the residual part of
pressure $\delta{p}^{\varpi}$ using (48):

$\delta\tilde{p}^{\varpi}\big{|}_{z=0}=\big{(}(F(\mu)+1)\lambda+iF(\mu)\Omega^{{\scriptscriptstyle{(0)}}}\big{)}\delta{\psi},$
(56)

where a function

$F(\mu)=1+\Big{(}\mu-\dfrac{1}{\mu}\Big{)}\dfrac{{\mathrm{J}}_{1}(\mu)}{{\mathrm{J}}_{0}(\mu)}$
(57)

is even in $\mu$. It follows from (55) also, that
$\partial_{z}\delta\tilde{V}_{z}|_{z=0}=-i\sqrt{\varepsilon}({\mathrm{J}}_{1}(\mu)/{\mathrm{J}}_{0}(\mu))k^{2}\theta\,\delta\psi$
in (46). The linear system of equations (46,47) for $\delta h$, $\delta\psi$
has a nontrivial solution if $F(\mu)=0$ in the first precision with respect to
small parameters $\varepsilon$ and $\theta$. Among all existing solutions we
should choose that corresponding to positive and minimal real part of
$\lambda$. The numerical solution is

$\mu_{\star}\simeq
1.67,\quad\lambda=\dfrac{\sqrt{2\Omega^{\scriptscriptstyle{(0)}}\omega}}{\mu_{\star}}\,kh^{\scriptscriptstyle{(0)}}|\theta|.$
(58)

As $\mu_{\star}<\mu_{\scriptscriptstyle CL2}$, obtained growth increment
$\lambda$ (58) is greater than that found in [6]. Hence, the unstable mode
does contain the oblique wave component. Knowing that $\mu\sim 1$, let us
rewrite the inequality (49). We introduce friction velocity $u_{\ast}$
according to definition
$\Omega^{\scriptscriptstyle{(0)}}=u_{\ast}^{2}/\nu_{\scriptscriptstyle T}$.
Then (49) means that Langmuir number

$\mathrm{La}\equiv\frac{\nu_{\scriptscriptstyle
T}^{3}k^{2}}{\omega(h^{\scriptscriptstyle(0)})^{2}u_{\ast}^{2}}\ll\theta^{2}.$
(59)

Let us also test the properties of the solution with respect to changing the
sign of angle $\theta$, that leads to $\mu\to-\mu$ as well. Without loss of
generality we can adopt that $\delta h$ is pure real, $\delta h<0$, and the
phase difference between the initial and oblique waves along streamwise
direction is small, $k\theta^{2}x\ll 1$, so phase difference
$\varphi^{\scriptscriptstyle(1)}-\varphi^{\scriptscriptstyle(0)}=k\theta y$ in
(36). Then both $\delta\tilde{V}_{z}$ (III.3) and $\delta\tilde{V}_{x}$ (50)
are pure real and even in $\theta$, so $\delta V_{z}\propto\cos(k\theta y)$
(including the sign) is even function of $\theta$ after $\delta\tilde{V}_{x}$,
whereas $\delta\tilde{V}_{y}$ is pure imaginary and odd in $\theta$, so
$\delta V_{y}\propto-\theta\sin(k\theta y)$ (the sign is taken at the surface)
is even function of $\theta$ as well. This means that the unstable mode is
twice degenerate; two possible configurations have identical slow flow spatial
structure and oblique waves with opposite spanwise wave numbers. The derived
spatial distribution is schematically plotted in Figure 1.

##  IV Discussion

Presented analysis of Langmuir instability assumes that the flow has exactly
periodic structure in space. In reality, the periodicity should be lost at
some length $L_{c}$. For the analysis to be valid, the length should be
greater than the distance which wave passes during the characteristic time
$1/\lambda$ along streamwise direction, that is $L_{c,x}>\omega/(k\lambda)$.
Taking into account (58,52), one can rewrite the inequality in the form
$L_{c,x}>1/(\sqrt{\varepsilon}\,\theta\,kh^{\scriptscriptstyle(0)}k)$, that
much exceeds the period in spanwise direction. In spanwise direction, the
periodicity should be kept at smaller length $L_{c,y}\sim\theta L_{c,x}$. If
the periodicity appears to be disrupted at shorter distances, then the oblique
wave is not correlated with the Langmuir circulation and one should drop it
out during the Langmuir instability analysis. As a result, one arrives back to
CL2-model [6], where only $x$-independence of the flow at distances greater
than $1/(k\theta)$ is assumed.

Considering an extension of the scope of our approach, we believe, that the
developed mathematical scheme should be appropriate for the analysis of wave-
current interaction in flow confined in a basin, see e.g. [25, 26]. The
confinement preserves correlation between vortical flow and surface waves,
which are mostly standing ones in the case. At initial stage, when only
surface waves are excited, the leading mechanism of the interaction is virtual
wave stress [20, 26], which can be enhanced due to presence of surface
contamination in a form of liquid elastic film [27]. However, during the
development of large-scale vortical flow, oblique waves appears that is
evidence of the wave scattering by the vortical flow [25]. The presented
approach, combined with a more detailed analysis of experimental data, can
answer whether a wave-current interaction loop plays any role in formation of
large-scale flow, along with other mechanisms, including the two-dimensional
inverse energy cascade [28].

##  V Acknowledgments

The work was supported by Russian Science Foundation (project no.
23-72-30006).

##  Appendix A Contribution into pressure related to high-frequency part of
the vorticity

The equation on oscillating part of pressure follows from Navier-Stokes
equation for wave motion

$\partial_{t}\bm{u}^{\phi}+\left(\bm{u}^{\phi}\nabla\right)\bm{V}+\left(\bm{V}\nabla\right)\bm{u}^{\phi}+\partial_{t}\bm{u}^{\varpi}+\left(\bm{u}^{\varpi}\nabla\right)\bm{V}+\left(\bm{V}\nabla\right)\bm{u}^{\varpi}=-\nabla
p^{u}+\bm{g}$ (60)

Let us group potential and vortical terms

$\Bigl{(}\partial_{t}\bm{u}^{\varpi}+\left(\bm{u}^{\varpi}\nabla\right)\bm{V}+\left(\bm{V}\nabla\right)\bm{u}^{\varpi}\Bigr{)}-\left[\bm{u}^{\phi}\bm{\Omega}\right]=-\nabla\Bigl{(}p^{u}+gz+\left(\partial_{t}+V_{i}\partial_{i}\right)\phi\Bigr{)},$
(61)

using

$\left(\bm{u}\nabla\right)\bm{V}+\left(\bm{V}\nabla\right)\bm{u}=\nabla\left(\bm{u}\bm{V}\right)-\left[\bm{V}\bm{\omega}\right]-\left[\bm{u}\bm{\Omega}\right],$

valid for incompressible vector fields. The gradient term in right-hand side
of (61) we interpret as vortical part of pressure which we denote as
$p^{\varpi}$,

$\displaystyle p^{\varpi}\equiv
p^{u}+gz+\left(\partial_{t}+V_{i}\partial_{i}\right)\phi$ (62)
$\displaystyle\Delta p^{\varpi}=f^{\varpi},\qquad
f^{\varpi}=\partial_{i}\phi\,\Delta
V_{i}+2\partial_{j}u^{\varpi}_{i}\,\partial_{i}V_{j}.$ (63)

The boundary condition for $p^{\varpi}$ on the free surface is following from
$z$-component of equation (61). Taking into account the smallness of spatial
derivatives of $\bm{V}$, we present

$\left.\partial_{z}p^{\varpi}\right|_{z=0}=\partial_{\alpha}\phi\,\partial_{z}V_{\alpha},\qquad\left.p^{\varpi}\right|_{z\rightarrow-\infty}\rightarrow
0.$ (64)

The differential problem (63), (64) could be solved analytically in Fourier
space (${\bf q}$ denotes Fourier harmonic in $Oxy$-здате)

$\left(\partial^{2}_{z}-q^{2}\right)\delta{p}^{\varpi}_{\bf q}=f^{\varpi}_{\bf
q},\qquad\delta{p}^{\varpi}_{\bf
q}=-\dfrac{e^{qz}}{2q}\int\limits_{-\infty}^{0}dz^{\prime}e^{qz^{\prime}}f^{\varpi}_{\bf
q}(z^{\prime})-\int\limits_{-\infty}^{0}dz^{\prime}\dfrac{e^{-q|z-z^{\prime}|}}{2q}f^{\varpi}_{\bf
q}(z^{\prime})+\frac{e^{qz}}{q}\big{(}\partial_{\alpha}\phi\partial_{z}V_{\alpha}\big{)}_{\bf
q}.$ (65)

Now we turn to calculations in the sake of Langmuir instability analysis.
Using the estimation (22) we get $f^{\varpi}=\partial_{i}\phi\,\Delta V_{i}$
in (63). Thus,

$\delta{p}^{\varpi}_{\bf
q}\big{|}_{z=0}=-\dfrac{1}{q}\int\limits_{-\infty}^{0}dz\,e^{qz}\big{(}\partial_{i}\phi\Delta
V_{i}\big{)}_{\bf
q}+\dfrac{1}{q}\partial_{\alpha}\big{(}\phi\partial_{z}V_{\alpha}\big{)}_{\bf
q}.$ (66)

Using this expression, let us prove the equality (35). We substitute (33) and
(36) in (66) and by direct calculations we show that only harmonics
oscillating as $\exp[ikx\cos\theta+iky\sin\theta]$ appears nonzero (note that
we should set ${\bf q}=k\\{\cos\theta,\sin\theta\\}$). According to definition
(35) we obtain

$\displaystyle\delta\tilde{p}^{\varpi}\big{|}_{z=0}=-\psi^{{\scriptscriptstyle{(0)}}}\int\limits_{-\infty}^{0}dz\,e^{2kz}\big{(}\partial_{z}^{2}-k^{2}\sin^{2}\theta\big{)}\big{(}i\delta\tilde{V}_{x}+\delta\tilde{V}_{z}\big{)}+i\delta\psi\Omega^{\scriptscriptstyle{(0)}}\big{|}_{z=0}\cos\theta+\psi^{{\scriptscriptstyle{(0)}}}(\partial_{x}+ik)\partial_{z}\delta\tilde{V}_{x}\Big{|}_{z=0}$
$\displaystyle+ik\psi^{{\scriptscriptstyle{(0)}}}\sin\theta\partial_{z}\delta\tilde{V}_{y}\Big{|}_{z=0}.$
(67)

## References

  * Hamlington et al. [2014] P. E. Hamlington, L. P Van Roekel, B. Fox-Kemper, K. Julien, and G. P. Chini.  Langmuir–submesoscale interactions: Descriptive analysis of multiscale frontal spindown simulations.  _Journal of Physical Oceanography_ , 44(9):2249–2272, 2014.
  * Thorpe [2004] S.A. Thorpe.  Langmuir circulation.  _Annu. Rev. Fluid Mech._ , 36:55–79, 2004.
  * Teixeira [2019] Miguel A.C. Teixeira.  Langmuir circulation and instability.  In J. Kirk Cochran, Henry J. Bokuniewicz, and Patricia L. Yager, editors, _Encyclopedia of Ocean Sciences (Third Edition)_ , pages 92–106. Academic Press, Oxford, 2019.
  * Weller and Price [1988] R. A. Weller and J. F. Price.  Langmuir circulation within the oceanic mixed layer.  _Deep Sea Research Part A. Oceanographic Research Papers_ , 35(5):711–747, 1988.
  * Chang et al. [2019] H. Chang, H. S. Huntley, A.D. Kirwan, D. F. Carlson, J. A. Mensa, S. M., G. Novelli, T. M. Özgökmen, B. Fox-Kemper, B. Pearson, et al.  Small-scale dispersion in the presence of langmuir circulation.  _Journal of physical oceanography_ , 49(12):3069–3085, 2019.
  * Craik [1977] A. D. D. Craik.  The generation of langmuir circulations by an instability mechanism.  _Journal of Fluid Mechanics_ , 81(2):209–223, 1977.
  * Leibovich [2009] S. Leibovich.  Langmuir circulation and instability.  _Elements of Physical Oceanography: A derivative of the Encyclopedia of Ocean Sciences_ , 124(2):288, 2009.
  * Faller and Caponi [1978] A. J. Faller and E. A. Caponi.  Laboratory studies of wind-driven langmuir circulations.  _Journal of Geophysical Research: Oceans_ , 83(C7):3617–3633, 1978.
  * Sullivan and McWilliams [2019] P. P. Sullivan and J. C. McWilliams.  Langmuir turbulence and filament frontogenesis in the oceanic surface boundary layer.  _Journal of Fluid Mechanics_ , 879:512–553, 2019.
  * Craik and Leibovich [1976] A. D. D. Craik and S. Leibovich.  A rational model for langmuir circulations.  _Journal of Fluid Mechanics_ , 73(3):401–426, 1976.
  * Craik [1982] A. D. D. Craik.  Wave-induced longitudinal-vortex instability in shear flows.  _Journal of Fluid Mechanics_ , 125:37–52, 1982.
  * Kawamura [2000] T. Kawamura.  Numerical investigation of turbulence near a sheared air–water interface. part 2: Interaction of turbulent shear flow with surface waves.  _Journal of marine science and technology_ , 5:161–175, 2000.
  * Veron and Melville [2001] F. Veron and W. K. Melville.  Experiments on the stability and transition of wind-driven water surfaces.  _Journal of Fluid Mechanics_ , 446:25–65, 2001.
  * Fujiwara and Yoshikawa [2020] Y. Fujiwara and Y. Yoshikawa.  Mutual interaction between surface waves and langmuir circulations observed in wave-resolving numerical simulations.  _Journal of Physical Oceanography_ , 50(8):2323–2339, 2020.
  * Leibovich [1977] S. Leibovich.  On the evolution of the system of wind drift currents and langmuir circulations in the ocean. part 1. theory and averaged current.  _Journal of Fluid Mechanics_ , 79(4):715–743, 1977.
  * Lamb [1932] H. Lamb.  _Hydrodynamics_.  University Press, 1932.
  * Stewart and Joy [1974] R. H. Stewart and J. W. Joy.  Hf radio measurements of surface currents.  In _Deep sea research and oceanographic abstracts_ , volume 21, pages 1039–1049. Elsevier, 1974.
  * Longuet-Higgins [1953] M. S. Longuet-Higgins.  Mass transport in water waves.  _Mathematical and Physical Sciences_ , 245(903):535–581, 1953.
  * Nicolás and Vega [2003] J. A. Nicolás and J. M. Vega.  Three-dimensional streaming flows driven by oscillatory boundary layers.  _Fluid Dynamics Research_ , 32(4):119–139, 2003\.
  * Filatov et al. [2016] S. V. Filatov, V. M. Parfenyev, S. S. Vergeles, M. Yu. Brazhnikov, A. A. Levchenko, and V. V. Lebedev.  Nonlinear generation of vorticity by surface waves.  _Physical Review Letters_ , 116(5):054501, 2016\.
  * Garrett [1976] C. Garrett.  Generation of langmuir circulations by surface waves — a feedback mechanism.  _Journal of Marine Research_ , 34:117, 1976.
  * Abrashkin et al. [2022] A. A. Abrashkin, E. N. Pelinovsky, et al.  Gerstner waves and their generalizations in hydrodynamics and geophysics.  _Phys. Usp._ , 65:453–467, 2022.
  * Leibovich [1983] S. Leibovich.  The form and dynamics of langmuir circulations.  _Annual Review of Fluid Mechanics_ , 15(1):391–427, 1983.
  * Craik [1970] A. D. D. Craik.  A wave-interaction model for the generation of windrows.  _Journal of Fluid Mechanics_ , 41(4):801–821, 1970.
  * Filatov et al. [2018] S. V. Filatov, A. V. Orlov, M. Yu. Brazhnikov, and A. A. Levchenko.  Experimental simulation of the generation of a vortex flow on a water surface by a wave cascade.  _JETP Letters_ , 108:519–526, 2018.
  * Filatov et al. [2022] S.V. Filatov, A.V. Poplevin, A.A. Levchenko, and V.M. Parfenyev.  Generation of stripe-like vortex flow by noncollinear waves on the water surface.  _Physica D: Nonlinear Phenomena_ , 434:133218, 2022.
  * Parfenyev et al. [2019] V.M. Parfenyev, S.V. Filatov, M. Yu. Brazhnikov, S.S. Vergeles, and A.A. Levchenko.  Formation and decay of eddy currents generated by crossed surface waves.  _Physical Review Fluids_ , 4(11):114701, 2019\.
  * Colombi et al. [2022] R. Colombi, N. Rohde, M. Schlüter, and A. von Kameke.  Coexistence of inverse and direct energy cascades in faraday waves.  _Fluids_ , 7(5):148, 2022.


research print

Jessica K. Hodgins College of Computing
801 Atlantic Drive
Georgia Institute of Technology
Atlanta, GA 30332-0280 (404) 894-9763 (404) 894-0673<EMAIL_ADDRESS>
8

# Animating Explosions

Gary D. Yngve  James F. O’Brien  Jessica K. Hodgins  GVU Center and College of
Computing
Georgia Institute of Technology

###### Abstract

In this paper, we introduce techniques for animating explosions and their
effects. The primary effect of an explosion is a disturbance that causes a
shock wave to propagate through the surrounding medium. This disturbance
determines the behavior of nearly all other secondary effects seen in
explosions. We simulate the propagation of an explosion through the
surrounding air using a computational fluid dynamics model based on the
equations for compressible, viscous flow. To model the numerically stable
formation of shocks along blast wave fronts, we employ an integration method
that can handle steep pressure gradients without introducing inappropriate
damping. The system includes two-way coupling between solid objects and
surrounding fluid. Using this technique, we can generate a variety of effects
including shaped explosive charges, a projectile propelled from a chamber by
an explosion, and objects damaged by a blast. With appropriate rendering
techniques, our explosion model can be used to create such visual effects as
fireballs, dust clouds, and the refraction of light caused by a blast wave.

###### keywords:

Animation, Atmospheric Effects, Computational Fluid Dynamics, Natural
Phenomena, Physically Based Animation

I.3.5Computer GraphicsComputational Geometry and Object ModelingPhysically
based modeling; I.3.7Computer GraphicsThree-Dimensional Graphics and
RealismAnimation; I.6.8Simulation and ModelingTypes of SimulationAnimation

00footnotetext:  College of Computing, Georgia Institute of Technology,
Atlanta, GA 30332\<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
SIGGRAPH 2000, New Orleans, LA USA
This is an author’s preprint.
https://doi.org/10.1145/344779.344801

##  1 Introduction

Explosions are among the most dramatic phenomena in nature. A sudden burst of
energy from a mechanical, chemical, or nuclear source causes a pressure wave
to propagate outward through the air. The blast wave “shocks up,” creating a
nearly discontinuous jump in pressure, density, and temperature along the wave
front. The wave is substantially denser than the surrounding fluid, allowing
it to travel supersonically and to cause a noticeable refraction of light. The
air at the shock front compresses, turning mechanical energy into heat. The
waves reflect, diffract, and merge, allowing them to exhibit a wide range of
behavior.

An explosion causes a variety of visual effects in addition to the light
refraction by the blast wave. An initial chemical or nuclear reaction often
causes a blinding flash of light. Dust clouds are created as the blast wave
races across the ground, and massive objects are moved, deformed, or
fractured. Hot gases and smoke form a rising fireball that can trigger further
combustion or other explosions and scorch surrounding objects.

We present a physically based model of an explosion and show how it can be
used to simulate many of these effects. We model the explosion post-detonation
as compressible, viscous flow and solve the flow equations with an integration
method that handles the extreme shocks and supersonic velocities inherent in
explosions. We cannot capture many of the visual effects of an explosion in a
complex setting if we rely only on an analytical model of the blast wave; a
fluid dynamics model of the air is necessary to capture these effects. The
system includes a two-way coupling between dynamic objects and fluid that
allows the explosions to move objects. Figure 1 illustrates this phenomenon
with a projectile propelled from a chamber. We also use the pressure wave
generated by the explosion to fracture and deform objects. The user can
simulate arbitrarily complex scenarios by positioning polygonal meshes to
represent explosions and objects. The user controls the scale and visual
qualities of the explosion with a few physically motivated parameters.

Figure 1:  An image of a projectile propelled from a chamber by an explosion.
On the right is a cross-section of the three-dimensional fluid volume using a
colormap where hotter colors indicate higher densities.

Figure 2:  Shown here is a cross-section of pressures for a three-dimensional
explosion near an immovable wall. The timestep between frames is 10 ms. Many
of the behaviors of the blast wave can be seen, specifically the Mach stem
formed from the blast wave merging with the wave reflected off the ground
plane, and the diffracted wave formed when the blast wave crests over the
wall.

Our fluid model of an explosion simulates many phenomena of blast waves that
existing graphics techniques do not capture. Figure 2 shows a cross-section of
pressures for a three-dimensional explosion near a wall. The initial
disturbance in the first image interacts with the surrounding fluid and causes
a pressure wave to propagate through the medium. In the second image, the
blast wave has “shocked up,” as is evident by the large differences in
pressure across the shock front. The blast wave reflects off the wall and the
ground in the third image. In the fourth, the wave that reflected off the
ground merges with the initial blast wave to form a Mach stem, which has
pressure values twice that of the initial wave. In the final two images, the
blast wave crests over the wall and forms a weaker diffracted wave.

In the entertainment industry, explosions are currently created at full scale
in the real world, in miniature, or using heuristic graphics techniques[20, 9,
22]. Each of these methods has significant disadvantages, and we believe that
in many scenarios, simulation may provide an easier solution. When explosions
are generated and filmed at full scale in the real world, they often must be
faked to appear dangerous and destructive by using multiple charges and
chemicals with a low flashpoint. Because of the cost and danger of exploding
full-size objects, many explosions are created using miniatures. With
miniatures, the greatest challenge is often scaling the objects and the
physics to create a realistic effect. Current graphics techniques for creating
explosions are based on heuristics, analytical functions, or recorded data,
and although they produce nice effects for spherical blast waves, they are not
adequate for the complex effects required for many of the scenarios used in
the entertainment industry.

Physically based simulations of explosions offer several potential advantages
over these three techniques. In contrast to real physical explosions,
simulations can be used in an iterative fashion, allowing the director many
chances to modify or shape the effect. The rendering of the explosion is to a
large extent decoupled from the simulation, allowing the visual
characteristics of the dust clouds or fireball to be determined as a post-
process. Unlike heuristic or analytical graphical methods, physically based
simulations allow the computation of arbitrarily complex scenes with multiple
interacting explosions and objects.

The next section discusses relevant previous work in explosions, fluids,
flame, and fracture. The following section introduces the explosion model in
the context of computational fluid dynamics. The next two sections discuss
coupling between the fluid and solids and other secondary effects such as
refraction and fireballs. We close with a discussion of our results.

##  2 Previous Work

Explosions used in the entertainment industry tend to be visually rich.
Because of the inherent computational complexity of these explosions,
researchers largely neglected this field after the publication of particle
simulation techniques[16, 17]. Procedural methods can generate fiery, billowy
clouds that could be used as explosions[3].

Recently two papers specifically addressed explosions. Mazarak and colleagues
simulate the damage done by an explosion to voxelized objects[10]. They model
the explosion as an ideal spherical blast wave with a pressure profile curve
approximated by an analytic function based on the modified Friedlander
equation and scaled according to empirical laws[2]. The spherical blast wave
expands independent of existing obstacles, and forces are applied to objects
in the direction of the blast radius. Objects are modeled as connected voxels
and based on various heuristics, these radial forces may cause the voxels to
disconnect.

Neff and Fiume use data from empirical blast curves to model an explosion[12].
The blast curves relate the pressure and velocity of the blast wave to time
and are scalable. Unlike Mazarak and colleagues, they use a curve representing
the reflection coefficient to apply forces to objects based on the angle of
incidence of the blast wave. They assume quasi-static loading conditions where
the blast wave encloses the entire object and effects due to reflected waves
are ignored. They also model explosion-induced fracture in planar surfaces
using a procedural pattern generator.

An alternative to these analytic and empirical models is a computational fluid
dynamic simulation of the blast wave and the surrounding air. Foster and
Metaxas presented a solution for incompressible, viscous flow and used it to
animate liquids[5] and hot, turbulent gas[4]. They modeled fluid as a three-
dimensional voxel volume with appropriate boundary conditions. The fluid obeys
the Navier-Stokes equations; gas also follows an equation that represents
thermal buoyancy. Using an explicit scheme, they update velocities and
temperatures every timestep via Euler integration and readjust the values to
guarantee conservation of mass. The fluid is rendered by tracing massless
particles along the interpolated flow field. Their work with liquids included
dynamic objects that were moved by the fluid, although they assumed that the
objects were small enough not to influence the fluid. Recently Stam addressed
the computational cost of guaranteeing stability by introducing extra damping
to afford larger timesteps and using an implicit method to solve a sparse
system of equations[18]. Stam’s method is inappropriate for shocks and
explosions because his integration scheme achieves stability by encouraging
the fluid to dissipate.

In the dramatic effects produced by the entertainment industry, a fireball is
often the most salient visible characteristic of an explosion. Stam and Fiume
modeled flame and the corresponding fluid flow and rendered the results using
a sophisticated global illumination method[19]. The gases behaved according to
advection-diffusion equations; Stam and Fiume solve these equations
efficiently by reformulating the problem from a grid to “warped blobs.”
Illumination from gas is affected by emission and anisotropic scattering and
absorption. They only consider continuous emissions from blackbody radiation
and ignore line emissions from electron excitation. They develop a heuristic
for smoke emission due to the lack of a scientific analytic model.

Compressible flow has been studied for years in the computational fluid
dynamics community[1, 2, 8]. We have built on this work by taking the
governing equations and the donor-acceptor method of integration from this
literature. However, the reasons for simulating explosions, combustion,
detonation, and supersonic flow in engineering differ significantly from those
in computer graphics. Engineering problems often require focusing on one
element such as the boundary layer and simulating the other elements only to
the extent that they affect the phenomenon under study. For example,
engineering simulations are often two-dimensional and assume symmetry in the
third dimension. Because they are focused on a specific event, their
simulations may run for only a few microseconds. In computer graphics, on the
other hand, we need to produce a visually appealing view of the behavior
throughout the explosion. As a result, we need a more complete model with less
quantitative accuracy.

##  3 Explosion Modeling

An explosion is a pressure wave caused by some initial disturbance, such as a
detonation. In the results presented here, we assume that the detonation has
occurred and that its properties are defined in the initial conditions of the
simulation. This assumption is reasonable for most chemical explosions because
the detonation is complete within microseconds. We animate explosions by
modeling the pressure wave and the surrounding air as a fluid discretized over
a three-dimensional rectilinear grid. The following two subsections describe
the governing equations for fluid dynamics and the computational techniques
used to solve them. The remaining two subsections describe the parameters
available to the user for controlling the appearance of the explosion via the
boundary conditions and initial conditions.

###  3.1 Fluid Dynamics

In nearly all engineering problems, including the analysis of explosions,
fluids are modeled as a continuum. They are represented as a set of equations
in terms of density $\rho$ ($\mbox{$\mathrm{kg}$}/\mbox{$\mathrm{m}$}^{3}$),
pressure $P$ ($\mathrm{N}/\mathrm{m}^{2}$), velocity $v$
($\mathrm{m}/\mathrm{s}$), temperature $T$ ($\mathrm{K}$), the internal energy
per unit mass $N$ ($\mathrm{J}/\mathrm{kg}$), and the total energy per unit
mass $E=N+\frac{1}{2}\mbox{\boldmath$v$}^{2}$ ($\mathrm{J}/\mathrm{kg}$). The
equations that govern these quantities are defined in an Eulerian fashion,
that is, they apply to a differential volume of space that is filled with
fluid rather than to the fluid itself. In addition to the Navier-Stokes
equations, which model the conservation of momentum, the equations for
compressible, viscous flow include governing equations for the conservation of
mass and energy and for the fluid’s thermodynamic state[7].

We introduce several simplifying assumptions that make the equations easier to
compute but nevertheless allow us to capture the effects of compressible,
viscous flow. We discount changes in the vibrational energies of molecules and
assume air to be at chemical equilibrium; we ignore the effects from
dissociation or ionization. These assumptions, which are commonly used in the
engineering literature[1], allow us to reduce to constants many of the
coefficients that vary with temperature. The resulting deviation in the values
of the coefficients is negligible at temperatures below
$1000\,\mbox{$\mathrm{K}$}$; only minor deviations occur below
$2500\,\mbox{$\mathrm{K}$}$. Our implementation produces aesthetic results
with temperatures above $100000,\mbox{$\mathrm{K}$}$, although deviations in
constants could be on the order of a magnitude or two.

The first governing equation of fluid dynamics arises from the conservation of
mass. Because fluid mass is conserved, the change of fluid density in a
differential volume must be equal to the net flux across the volume’s
boundary, giving

$\frac{\partial\rho}{\partial t}=-\nabla\cdot(\rho{\mbox{\boldmath$v$}}).$ (1)

The second governing equation, commonly known as the Navier-Stokes equation,
concerns the conservation of momentum. For a Stokes fluid, where the normal
stress is independent of the rate of dilation, the equation for the $x$
component of the fluid velocity is given by

$\rho\frac{\partial\mbox{\boldmath$v$}_{x}}{\partial
t}=\rho\mbox{\boldmath$f$}_{x}-\nabla
P+\frac{\mu}{3}\nabla\cdot\left(\frac{\partial\mbox{\boldmath$v$}}{\partial
x}\right)+\mu\nabla^{2}\mbox{\boldmath$v$}_{x}-\rho(\mbox{\boldmath$v$}\cdot\nabla)\mbox{\boldmath$v$}_{x},$
(2)

where $f$ represents the body forces such as gravity and $\mu$ is the
coefficient of viscosity. The equations for the $y$ and $z$ components are
similar. The first two terms on the right-hand side of the equation model
accelerations due to body forces and forces arising from the pressure
gradient; the next two terms model accelerations due to viscous forces. The
last term is not a force-related term; rather it is a convective term that
models the transport of momentum as the fluid flows. This distinction between
time derivative (force) terms and convective terms will be important for the
integration scheme.

The final governing equation enforces the conservation of energy in the
system. The First Law of Thermodynamics dictates that the change in energy is
equal to the amount of heat added and the work done to the system. Accounting
for the amount of work done from pressure and viscosity and the heat
transferred from thermal conductivity yields

$\rho\frac{\partial N}{\partial
t}=k\nabla^{2}T-P\nabla\cdot\mbox{\boldmath$v$}+\Phi-\rho(\mbox{\boldmath$v$}\cdot\nabla)N,$
(3)

where $k$ is the thermal conductivity constant and $\Phi$ is the viscous
dissipation given by

$\Phi=-\frac{2\mu}{3}(\nabla\cdot\mbox{\boldmath$v$})^{2}+\frac{\mu}{2}\sum_{i,j\in\\{x,y,z\\}}\left(\frac{\partial\mbox{\boldmath$v$}_{i}}{\partial
j}+\frac{\partial\mbox{\boldmath$v$}_{j}}{\partial i}\right)^{2}.$ (4)

As with equation (2), the last term of equation (3) is a convective term and
models the transport of energy as the fluid flows.

In addition to the three governing equations, we need equations of state that
determine the relationship between energy, temperature, density, and pressure.
They are

$\displaystyle N=c_{V}T,\;\;P=\rho RT,$ (5)

where the coefficient $c_{V}$ is the specific heat at constant volume and $R$
is the gas constant of air.

Figure 3:  Pressure profile ($\mathrm{N}/\mathrm{m}^{2}$) over time
($\mathrm{ms}$) near an explosion.

Figure 4:  This figure illustrates the donor-acceptor method in which the
amount of mass transferred is proportional to the mass of the donor. The
voxels on the left show the transfer of mass and energy according to the flow,
indicated by the blue arrow. The two voxels on the right represent the
scenario with reversed flow of the same magnitude. Density is represented as
height, and unit energy is represented as color. Corresponding amounts of
energy are sent with the mass.

###  3.2 Discretization and Numerical Integration

The equations in the previous section describe the behavior of a fluid in a
continuous fashion. However, implementing them in a form suitable for
numerical computation requires that the space filled by the fluid be
discretized in some manner and that a stable method for integrating the
governing equations forward in time be devised.

Finite differences are used to discretize the space into a regular lattice of
cubical cells. These finite voxels take the place of the differential volumes
used to define the continuous equations, and the governing equations now hold
for each voxel. Fluid properties such as pressure and velocity are associated
with each voxel and these properties are assumed to be constant across the
voxel. The spatial derivatives used in the governing equations are
approximated on the lattice using central differences. For example, the $x$
component of the pressure gradient, $\nabla P$, at voxel $[i,j,k]$ is given by

$\frac{\partial P}{\partial x}\approx\frac{P_{[i+1,j,k]}-P_{[i-1,j,k]}}{2h},$
(6)

where subscripts in square brackets index voxel locations and $h$ is the voxel
width.

After the governing equations have been expressed in terms of discrete
variables using finite differences, they may be used as the update rules for
an explicit integration scheme. However, rapid pressure changes created by
steep pressure gradients moving at supersonic speeds would cause such a scheme
to diverge rapidly. (See Figure 3.) To deal with this problem, we improved the
basic integration technique using two modifications described in the fluid
dynamics literature[2, 8]. The first modification involves updating equations
(2) and (3) in two steps, first using only the temporal portion of the
derivatives and second using the convective derivatives. The second
modification is called the donor-acceptor method and is described in detail
below. It addresses problems that arise when mass, momentum, and energy are
convected across steep pressure gradients.

The modified update scheme operates by applying the following algorithm to
each voxel at every timestep:

  1. 1.

Approximate the fluid acceleration at the current time,
$\mbox{\boldmath$\widetilde{a}$}_{t}=\left({\partial{\mbox{\boldmath$v$}}}/{\partial{t}}\right)_{t}$,
using the non-convective (first four) terms of equation (2).

  2. 2.

Compute the tentative velocity at the end of the timestep,
$\mbox{\boldmath$\widetilde{v}$}_{t+\Delta t}=\mbox{\boldmath$v$}_{t}+\Delta
t\mbox{\boldmath$\widetilde{a}$}_{t}$, and the approximate average velocity
during the timestep
$\mbox{\boldmath$\bar{v}$}_{t}=(\mbox{\boldmath$\widetilde{v}$}_{t+\Delta
t}+\mbox{\boldmath$v$}_{t})/2$.

  3. 3.

Approximate change in internal energy, $N$, using the non-convective terms of
equation (3) and substituting $\mbox{\boldmath$\bar{v}$}_{t}$ for the fluid
velocity.

  4. 4.

Using $\mbox{\boldmath$\bar{v}$}_{t}$ for the fluid velocity, compute the new
density, $\rho_{t+\Delta t}$ with equation (1).

  5. 5.

Calculate the complete $\mbox{\boldmath${v}$}_{t+\Delta t}$ and $N_{t+\Delta
t}$ with equations (2) and (3) using the convective terms and the new value of
$\rho$.

  6. 6.

Use state equations (5) to update secondary quantities such as temperature.

Although this update scheme is more stable than a simple Euler integration,
sharp gradients in fluid density may still allow small flows from nearly empty
voxels to generate negative fluid densities and cause inappropriately large
changes to both velocity and internal energy. To prevent these problems, we
use a donor-acceptor method when computing $-\nabla\rho\mbox{\boldmath$v$}$ of
the convective terms in steps 4 and 5 above.

The donor-acceptor method transfers mass proportional to the mass of the voxel
in the upstream direction, or donor voxel. Suppose we have voxel $i$ and one
of its six neighbors $j$ in direction $d$ from $i$. Let
$v_{ij}=\frac{1}{2}(\mbox{\boldmath$v$}_{i}+\mbox{\boldmath$v$}_{j})\cdot\mbox{\boldmath$d$}$.
If $v_{ij}>0$, then flow is going from $i$ to $j$, $i$ is the donor, and
$\rho_{i}$ is used in equation (1) to compute the change in voxel $i$’s
density. Likewise, if $v_{ij}<0$, then $j$ is the donor, and the density of
$j$ is used when computing the change in voxel $i$’s density. These
calculations are repeated for the six neighbors of $i$ to obtain the new
density for $i$, $\rho_{t+\Delta t}$. The velocity and energy convection are
then scaled by ${\rho_{t}}/{\rho_{t+\Delta t}}$ to conserve momentum and
energy. Figure 4 illustrates the donor-acceptor method. The left and right
diagrams show flows in opposite directions of the same magnitude. The sent
mass is proportional to the mass of the donor and carries with it
corresponding amounts of energy.

Figure 5:  This figure shows cross-sections of pressure for three-dimensional
explosions of equal-volume charges in the shape of a sphere, cylinder, torus,
and wedge. The timestep between frames is 5 ms.

###  3.3 Boundary Conditions

The system has several types of boundary conditions that allow the fluid to
exhibit a wide range of behaviors. Free boundaries allow blast waves to travel
beyond the voxel volume as if the voxel volume were arbitrarily large. This
type of boundary allows us to model slow, long-term aspects of explosions,
such as fireballs and dust clouds. Hard boundaries force fluid velocity normal
to them to be zero while leaving all other fluid attributes unchanged. We
treat these boundaries as smooth surfaces, so tangential flow is unaffected.
We implement a third boundary condition to achieve faster execution. If a
voxel and its neighbors have pressure differences less than a threshold, the
voxel is treated as a free boundary and is never evaluated. This optimization
prunes out the majority of the volume while the blast wave is expanding.

###  3.4 Initial Conditions

The user specifies the pressure and temperature of the air, and the initial
values of other variables are determined from the state equations (5). The
detonation is initialized by specifying a region of the volume with higher
temperature or pressure. For example, a chemical explosion might have a
temperature of $2900\,\mbox{$\mathrm{K}$}$ and a pressure of
$1000\,\mbox{$\mathrm{atm}$}$ with the surrounding air at
$290\,\mbox{$\mathrm{K}$}$ and $1\,\mbox{$\mathrm{atm}$}$. The creation of the
explosion may be time-delayed or may be triggered when the fluid around the
charge reaches a threshold temperature.

The detonation may have an arbitrary geometry represented by a manifold
polygonal mesh. The mesh is voxelized to initialize the appropriate voxels in
the fluid simulation. By controlling the geometry, the user can produce a
variety of effects that could not be achieved with a spherical model. In blast
theory, planar, cylindrical, and spherical blast waves can be modeled by
analytic functions[2]; however nonstandard shapes can create surprising and
interesting effects. Figure 5 shows cross-sections of pressure for three-
dimensional explosions from equal-volume charges in the shape of a sphere,
cylinder, torus, and wedge. The inner blast wave of the torus merges to create
a strong vertical blast wave. The wedge concentrates its force directly to the
right, while leaving the surrounding area relatively untouched.

##  4 Interaction with Solids

People use explosions to impart forces on objects for both constructive and
destructive purposes. The movements of these objects and the resulting
displacement of air create many of the compelling visual effects of an
explosion. In this section, we present methods to implement a two-way coupling
between the fluid and solids. The coupling from fluid to solid allows us to
model phenomena such as a projectile being propelled by an explosion. The
coupling from solid to fluid can be used to model a piston compressing or the
shock wave formed as a projectile moves through the air supersonically. We
also extend previously published techniques for fracture to allow the pressure
wave to shatter objects.

To allow the two-way coupling, objects have two representations: a polygonal
mesh that is used to apply forces to the object from the fluid, and a volume
representation in voxels that is used to displace fluid based on the motion of
the object. We incorporate the coupling into the fluid dynamics code in the
following way:

  1. 1.

Apply forces on the objects from the fluid and compute the rigid body motion
of the objects.

  2. 2.

If the object has moved more than a fraction of a voxel, recompute the
voxelization of the object.

  3. 3.

Displace fluid based on object movement.

  4. 4.

Update the fluid using the techniques described in Section 3.2.

We explain the first three of these items in greater detail in the following
subsections.

###  4.1 Coupling from Fluid to Solid

An object embedded in a fluid experiences two separate sets of forces on its
surface, those arising from hydrostatic pressure and those arising from
dynamic forces due to fluid momentum. The forces due to hydrostatic pressure
act normal to the surface and are generated by the incoherent motions of the
fluid molecules against the surface. The dynamic forces are generated by the
coherent motion of the continuous fluid and can be divided into a force normal
to the surface of the object and a tangential shearing force. We neglect the
tangential shearing force because in the context of explosions, it is
negligible in comparison to the force due to hydrostatic pressure. We assume
that the object is in equilibrium under ambient air pressure and the
hydrostatic forces are computed using the overpressure $\bar{P}$, which is the
difference between the hydrostatic pressure $P$ and ambient pressure.

The magnitude of the normal force per unit area on the surface is given by the
dynamic overpressure:

$\bar{P}_{\mathit{dyn}}=\bar{P}+\frac{1}{2}\rho\left(\mbox{\boldmath$v$}_{\mathit{rel}}\cdot\mbox{\boldmath$\hat{n}$}\right)^{2},$
(7)

where $\mbox{\boldmath$v$}_{\mathit{rel}}$ is the velocity of the fluid
relative to the surface, and $\hat{n}$ is the outward surface normal.

We assume that the triangles composing each object are small enough that the
force is constant over each triangle. The force on a triangle with area $A$ is
then

$\mbox{\boldmath$f$}=-\mbox{\boldmath$\hat{n}$}A\bar{P}_{\mathit{dyn}},$ (8)

where the fluid properties are measured at the centroid of the triangle. The
forces are computed for all triangles of an object, and the translational and
angular velocities of the object are updated accordingly.

Figure 6:  A glass window is shattered by a blast wave. The blast wave
pressure is approximately $3\,\mbox{$\mathrm{atm}$}$ when it reaches the
window. The images show the scene at 0 ms, 13 ms, 40 ms, 67 ms, 107 ms, and
160 ms.

Figure 7:  A wall is knocked over by a blast wave from an explosion 3.5 m
away. The images are spaced 250 ms apart.

In addition to acting on rigid objects, the forces can also be applied to
flexible objects that deform and fracture[14]. The explosion simulation
results in pressures, velocities, and densities for each voxel in the
discretization of the fluid. The fracture simulation uses this information to
compute the forces that should be applied to a finite element model of the
objects in the scene. The force computation is similar to that for rigid
objects. This method was used to simulate the breaking window shown in Figure
6 and the breaking wall shown in Figure 7. This coupling is one-way in that
the fluid applies forces to the finite-element model, but the fluid is not
moved by the fragments that pass through it.

###  4.2 Coupling from Solid to Fluid

To allow the solid to displace fluid, the triangular mesh representing the
object is converted to voxels[13], which are then used to define the hard
boundaries in the fluid volume dynamically. The objects move smoothly through
the fluid, but because of the discrete nature of the voxelization, large
changes in the amount of fluid displaced may occur on each timestep. To
address this problem, the fluid displaced or the void created by the movement
of the objects is handled over a period of time rather than instantaneously.

The voxelization returns a value between zero and one representing the
proportion of the voxel that is not interior to any object. This value is
independent of geometric considerations about the exact shape of the occupied
volume. If any dimension of an object is smaller than the size of a voxel, the
appropriate voxels will have partial volumes, but because there are no fully
occupied voxels, the fluid will appear to move through the object. Nonzero
partial volumes below a certain threshold are set to zero to increase
stability. The implementation of partial volumes requires slight modifications
to the donor-acceptor method to conserve mass, momentum, and energy because
two adjacent voxels could have different volumes.

When any part of an object moves more than a fraction of a voxel, the object
is revoxelized, and the hard boundaries of the fluid are updated. When this
process occurs, the partial volume in a voxel might change, resulting in fluid
flow. We allow this flow to occur smoothly by sacrificing conservation of mass
and energy in the short-term. The voxelization determines the partial volumes
in an instantaneous fashion, but the fluid displacement routine maintains
internal partial volumes that change more slowly and are used to compute the
pressure, density, and temperature of the affected voxels. The internal
partial volumes change proportionally to the velocities of the moving objects,
and mass and energy are restored over time.

To compute a smooth change in the internal partial volume from $V_{1}$ to
$V_{2}$, we model an object moving into a voxel as a piston compressing or
decompressing fluid. We simplify the computation of the change in partial
volume by assuming that the piston is acting along one of the axes of the
voxelization. The appropriate axis is selected based on the largest axial
component of the velocity of the object, $\mbox{\boldmath$v$}_{p}$. The
displacement of the piston after $t$ seconds and the corresponding change in
partial volume of a voxel with width $h$ are

$\mbox{\boldmath$l$}=\mbox{\boldmath$v$}_{p}t,\;\;\Delta
V=V_{2}-V_{1}=h^{2}\mbox{\boldmath$v$}_{p}t.$ (9)

The displacement occurs linearly over

$t={\Delta V\over h^{2}\mbox{\boldmath$v_{p}$}}$ (10)

at a velocity of $\mbox{\boldmath$v$}_{p}$.

Given this model of the change in internal partial volume, we know that
$\rho_{1}V_{1}=\rho_{2}V_{2}$ because mass is conserved. However, the fluid is
compressible, so mechanical energy is not conserved (otherwise
$P_{1}V_{1}=P_{2}V_{2}$). To obtain the new pressures and densities of the
fluid, we use a thermodynamic equation relating the work done to the system
from changing the volume (or density),

$\frac{P_{2}}{P_{1}}=\left(\frac{\rho_{2}}{\rho_{1}}\right)^{\gamma}=\left(\frac{T_{2}}{T_{1}}\right)^{\gamma/(\gamma-1)},$
(11)

where $\gamma=1+R/c_{V}$ ($\gamma$ is about 1.4 for air and is closer to 1 the
more incompressible the fluid)[1]. Internal unit energy and total unit energy
are then updated by the state equations.

When the partial volume of a voxel changes from one nonzero value to another,
the resulting pressure changes cause fluid to move to or from a neighbor based
on the governing equations. However, when the partial volume of a voxel
changes from zero to nonzero or vice versa, the situation must be handled as a
special case by treating the affected voxel and one of its neighbors as a
single larger voxel with a nonzero partial volume. The neighbor is selected
based on the largest axial component of the object’s velocity,
$\mbox{\boldmath$v$}_{p}$. We calculate the internal partial volume for each
of the involved voxels $A$ and $B$ as $\widetilde{V}_{A_{1}}$ and
$\widetilde{V}_{B_{1}}$.

When $V_{A_{1}}$, the original partial volume of $A$, is zero and $V_{A_{2}}$,
the new partial volume of $A$, is nonzero, we use initial volumes
$\widetilde{V}_{A_{1}}$ and $\widetilde{V}_{B_{1}}$ such that
$\widetilde{V}_{A_{1}}+\widetilde{V}_{B_{1}}=V_{B_{1}}$ (the initial volume is
conserved) and
$\widetilde{V}_{A_{1}}{V}_{B_{2}}=\widetilde{V}_{B_{1}}V_{A_{2}}$ (the voxels
are treated as a single larger voxel). The change in volume is

$\displaystyle\Delta{V}_{A}$ $\displaystyle=$
$\displaystyle{V}_{A_{2}}-\widetilde{V}_{A_{1}}={V}_{A_{2}}-{V}_{A_{2}}\frac{V_{B_{1}}}{{V}_{A_{2}}+{V}_{B_{2}}},$
(12) $\displaystyle\Delta{V}_{B}$ $\displaystyle=$
$\displaystyle{V}_{B_{2}}-\widetilde{V}_{B_{1}}={V}_{B_{2}}-{V}_{B_{2}}\frac{V_{B_{1}}}{{V}_{A_{2}}+{V}_{B_{2}}}.$
(13)

When the original partial volume of $A$, $V_{A_{1}}$, is nonzero and the new
partial volume of $A$, $V_{A_{2}}$, is zero, we force $\widetilde{V}_{A_{1}}$
to be zero and treat $\widetilde{V}_{B_{1}}$ as a single larger voxel. To
treat the two voxels as one, we first average the properties of $A$ and its
neighbor $B$, transferring any lost kinetic energy to internal energy. The
change in volume is then

$\Delta{V}_{B}=V_{B_{2}}-\widetilde{V}_{B_{1}}=V_{B_{2}}-{\rho_{A}{V}_{A_{1}}+\rho_{B}{V}_{B_{1}}\over\rho_{B}}$
(14)

by making sure that
$\rho_{B}\widetilde{V}_{B_{1}}=\rho_{A}V_{A_{1}}+\rho_{B}V_{B_{1}}$ (mass is
conserved).

##  5 Secondary Effects

Figure 8:  Refraction of light from a blast wave. Each frame is 10 ms apart.
The index of refraction is exaggerated tenfold to enhance the effect.

Figure 9:  A fireball after one second of simulation time. Tracer particles
from the fluid simulation determine the position and coloration of the
fireball.

An explosion creates a number of visual secondary effects including the
refraction of light, fireballs, and dust clouds. These secondary effects do
not significantly affect the simulation, so they can be generated and edited
as a post-process.

One of the most stunning, but often ignored, effects of an explosion is the
bending of light from the blast wave. Because the blast wave is substantially
denser than the surrounding air, it has a higher index of refraction, $\eta$.
Light travels at the same velocity between molecules, but near molecules it is
slowed down from interactions with electrons. This concept is expressed
numerically as $\eta-1=k\rho,\;k=2.26\times
10^{-4}\;\mathrm{m}^{3}/\mathrm{kg}$, by the Dale-Gladestone law[11]. We
capture the refraction of light by ray tracing through the fluid volume. (See
Figure 8.) As the ray is traced through the volume, the index of refraction is
continually updated based on the interpolated density of the current position.
For simplicity, we compute the density of each point using a trilinear
interpolation of the densities of the neighboring voxels. When the index of
refraction changes by more than a threshold, the new direction of the ray is
computed via Snell’s law using the density gradient as the surface normal. The
trilinear interpolation results in minor faceting effects that cause small
errors in the reflected direction.

Figure 10:  A fireball coming around a corner. The images are spaced 333 ms
apart.

An advantage of using a full volumetric fluid representation for explosions is
that the simulation can be used to model a fireball in addition to the blast
wave. We assume that the fireball is composed of detonated material from
inside the explosive. To track this material, the system initializes the
fireball by placing particles inside the shape specified by the user for the
explosion. The particles are massless and flow with the fluid, allowing the
fluid dynamics model to capture effects critical for a fireball such as
thermal conductivity and buoyancy. Some fluid simulations[4, 18] model thermal
buoyancy explicitly; in our simulation, thermal buoyancy is a behavior derived
from the governing equations. For rendering, each particle takes on a
temperature that is interpolated based on its position in the volume. The
particles are rendered as Gaussian blobs with values for red, green, blue, and
opacity. The color values are based on blackbody radiation at appropriate
wavelengths given the temperature of the particle[11]. Figure 9 shows a
fireball and corresponding tracer particles after one second of simulation.
Figure 10 shows a fireball coming around a corner; the hallway is illuminated
by the flames.

The tracer particles couple the appearance of the fireball to the motion of
the fluid, and although heat generated by the initial explosion is added to
the fluid model, any additional heat generated by post-detonation combustion
is ignored. Radiative energy released at detonation could also be modeled for
rendering. Much like the difficulties encountered with rendering the sun[15],
the high contrast of this effect may require contrast-reduction techniques
such as LCIS[21].

The blast wave and other secondary waves create dust clouds by disturbing fine
particles resting on surfaces. The creation of dust clouds is difficult to
quantify either experimentally or analytically, so the rate at which the dust
becomes airborne is left as a control for the animator. Once a dust particle
is airborne, its behavior is dictated by its size. The smaller it is, the more
it is influenced by drag forces and the less it is influenced by inertial
forces. Smaller dust particles have lower terminal velocities and exhibit more
Brownian motion. With the exception of coagulated particles, experiments
reveal that most dust particles are approximately ellipsoidal with low
eccentricity, and the particles do not orient themselves to the fluid flow[6].
The difference in dynamics between these particles and spherical particles is
not that significant, so we assume dust particles to be spherical. We
implement dust as metaparticles, each representing a Gaussian density of
homogeneous dust particles. The dust size for each metaparticle is chosen
according to size distributions from experimental data for blasted shale[6].
The metaparticles travel through the fluid as if single particles were located
at their centers. Their variances grow according to the mean Brownian
diffusion per unit time. Figure 11 shows dust clouds in a city scene.

Figure 11:  An explosion among buildings. The images are spaced 667 ms apart.

##  6 Results and Discussion

We ran the system with several scenarios. The physical constants used in the
simulation were constants for air that were taken from an engineering
handbook[7]. Table 1 shows the voxel width, timestep, total simulation time,
initial volume of the explosion (proportional to yield), and initial pressure
and temperature of the explosion. The timesteps $\Delta t$ increase by a
factor of five once the blast wave leaves the volume.

The simulations ran on a single 195 MHz R10K processor and used a $101\times
101\times 101$ volume. The running times per timestep varied considerably from
several seconds to two minutes because of the pruning described in Section
3.3. For coupling with fracture, I/O became a major factor because in each
iteration the entire volume was written to disk; however, using better
compression would reduce this expense. Running times of the simulations varied
from a few hours (Figure 5) to overnight (Figures 2, 6, 7, and 8) to a few
days (Figures 1, 9, 10, and 12).

Example | $h$ | $\Delta t$ | $t_{tot}$ | $V_{0}$ | $P_{0}$ | $T_{0}$
---|---|---|---|---|---|---
(figure) |  ($\mathrm{m}$) |  ($\mathrm{ms}$) |  ($\mathrm{ms}$) |  ($\mathrm{m}^{3}$) |  ($\mathrm{atm}$) |  ($\mathrm{K}$)
projectile (1) | 1.0 | 0.10 | 450 | 73.60 | 1000 | 2900
barrier (2) | 0.2 | 0.01 | 25 | 0.52 | 1000 | 2900
shapes (5) | 1.0 | 0.10 | 30 | 1000.00 | 1000 | 2900
fracture (6,7) | 0.2 | 0.02 | 20 | 0.52 | 1000 | 2900
fireball (8,9) | 1.0 | 0.10 | 1000 | 65.40 | 1000 | 2900
corner (10) | 1.0 | 0.10 | 10000 | 268.08 | 1000 | 2900
city (11) | 1.0 | 0.10 | 5000 | 65.40 | 1000 | 2900
nuclear (12) | 50.0 | 0.50 | 30000 | $\mathsf{9.1}\\!\\!\times\\!\\!\mathsf{10}^{\mathsf{7}}$ | 345 | $\mathsf{1}\\!\\!\times\\!\\!\mathsf{10}^{\mathsf{5}}$
Table 1:  Parameters for simulations: voxel width, timestep, total simulation
time, and initial volume, pressure, and temperature of detonation.

We use an explicit integration technique to compute the motion of the pressure
wave caused by the detonation. Despite its magnitude, the wave does not
transport fluid large distances. Previously, fluid dynamics has been used most
often in computer graphics to capture the effects of macroscale fluid
transport where the fluid does move a significant distance. Implicit
integration techniques with large timesteps are appropriate for these
situations because they achieve stability by damping high frequencies. The
propagation of the pressure wave in our stiff equations, however, is
characterized by these high frequencies and it is essential that they not be
artificially damped. We chose, therefore, to use an explicit integration
technique; however, an implicit integration technique could be used to
simulate the fireball and dust clouds after the blast wave and the secondary
waves have left the volume. Using an implicit integration technique in the
slow flow regime could allow larger timesteps and faster execution times.

We assume that the voxels in the fluid volume are of a size appropriate for
the phenomena that we wish to capture. In particular, if solid objects have a
dimension smaller than a voxel, then they will not create a hard boundary that
prevents fluid flow. For example, a wall that is thinner than a voxel will
permit the blast wave to travel through it because partial volumes do not
maintain any geometric information about the sub-voxel shape of the object.
The difficulty of a two-way coupling with fracturing objects stems from having
to model subvoxel cracks, which should allow flow to go through. If small
objects are required, the voxel size could be decreased or dynamic remeshing
techniques could be used to create smaller voxels in the areas around
boundaries.

There are effects from explosions that we have not investigated. Although
smoke is often a visible feature of an explosion that includes a fireball, we
do not have a physically based model for smoke creation. Incomplete combustion
at lower temperatures results in smoke, and that observation could be used as
a heuristic to determine where smoke should be created in the fireball and how
densely. Stam and Fiume used a similar heuristic model[19]. Textures of
objects could be modified to show soot accumulation and scorching over time.
Dust clouds are created when an object fractures or pulverizes. Dust could be
introduced into our system when the finite-element model produces small
tetrahedra or when cracks form.

We made several assumptions in constructing our model of explosions. Most
discounted effects that did not contribute noticeably to the final rendered
images; however, some could produce a noticeable change in behavior in certain
situations and may warrant further investigation. We only model the blast wave
traveling through air. However, waves travel through other media, including
solid objects, and complex interface effects occur when a wave travels between
two different media. For large-scale explosions, meteorological conditions
such as the change in pressure with respect to altitude or the interface
between atmospheric layers (the tropopause) may need to be considered.

Figure 12:  A large-scale high-temperature explosion resembling a nuclear
explosion: after 3 s, 6 s, 12 s, 24 s.

Our goal in this work has been to create a physically realistic model of
explosions. However, this model should also lend itself to creating less
realistic effects. Even though our model does not incorporate high-temperature
effects such as ionization, we can still obtain interesting results on high-
temperature explosions. The fireball in Figure 12 resulted from an initial
detonation at $10^{5}\,K$. Explosions used in feature films often include far
more dramatic fireballs than would occur in the actual explosions that they
purport to mimic. By using more tracer particles and adjusting the rendering
parameters of the fireballs, we should be able to reproduce this effect. Noise
could be added either to the velocity fields or particle positions post-
process to make the explosion look more turbulent. Similarly, explosions in
space are often portrayed as more colorful and violent than explosions that
occurred outside of the atmosphere should be. Imparting an initial outward
velocity to the explosion, turning off gravity, and increasing the thermal
buoyancy by modifying the state equations might create a similar effect.

##  7 Acknowledgments

This project was supported in part by NSF NYI Grant No. IRI-9457621,
Mitsubishi Electric Research Laboratory, and a Packard Fellowship. The second
author was supported by a Fellowship from the Intel Foundation.

## References

  * [1] J. D. Anderson Jr.  Modern compressible flow: with historical perspective.  McGraw-Hill, Inc., 1990.
  * [2] W. E. Baker.  Explosions in air.  University of Texas Press, 1973.
  * [3] D. Ebert, K. Musgrave, D. Peachy, K. Perlin, and S. Worley.  Texturing and Modeling: A Procedural Approach.  AP Professional, 1994.
  * [4] N. Foster and D. Metaxas.  Modeling the motion of a hot, turbulent gas.  Proceedings of SIGGRAPH 97, pages 181–188, August 1997.
  * [5] N. Foster and D. Metaxas.  Realistic animation of liquids.  Graphics Interface ’96, pages 204–212, May 1996.
  * [6] H. L. Green and W. R. Lane.  Particulate Clouds: Dusts, Smokes and Mists.  D. Van Nostrand Company, Inc., 1964.
  * [7] A. M. Kuethe and C. Chow.  Foundations of aerodynamics: bases of aerodynamic design.  John Wiley and Sons, Inc., 1998.
  * [8] C. L. Madder.  Numerical modeling of detonations.  University of California Press, 1979.
  * [9] K. H. Martin.  Godzilla: The sound and the fury.  Cinefex, pages 82–107, July 1998.
  * [10] O. Mazarak, C. Martins, and J. Amanatides.  Animating exploding objects.  Graphics Interface ’99, pages 211–218, June 1999.
  * [11] J. R. Meyer-Arendt.  Introduction to classical and modern optics.  Prentice-Hall, Inc., 1984.
  * [12] M. Neff and E. Fiume.  A visual model for blast waves and fracture.  Graphics Interface ’99, pages 193–202, June 1999.
  * [13] F.S. Nooruddin and G. Turk.  Simplification and repair of polygonal models using volumetric techniques.  Technical Report GIT-GVU-99-37, Georgia Institute of Technology, 1999\.
  * [14] J. F. O’Brien and J. K. Hodgins.  Graphical modeling and animation of brittle fracture.  Proceedings of SIGGRAPH 99, pages 137–146, August 1999.
  * [15] A. J. Preetham, P. Shirley, and B. E. Smits.  A practical analytic model for daylight.  Proceedings of SIGGRAPH 99, pages 91–100, August 1999.
  * [16] W. T. Reeves.  Particle systems—a technique for modeling a class of fuzzy objects.  ACM Transactions on Graphics, 2(2):91–108, April 1983.
  * [17] K. Sims.  Particle animation and rendering using data parallel computation.  Computer Graphics (Proceedings of SIGGRAPH 90), 24(4):405–413, August 1990.
  * [18] J. Stam.  Stable fluids.  Proceedings of SIGGRAPH 99, pages 121–128, August 1999.
  * [19] J. Stam and E. Fiume.  Depicting fire and other gaseous phenomena using diffusion processes.  Proceedings of SIGGRAPH 95, pages 129–136, August 1995.
  * [20] R. Street.  Volcano: Toasting the coast.  Cinefex, pages 56–84, September 1997.
  * [21] J. Tumblin and G. Turk.  LCIS: A boundary hierarchy for detail-preserving contrast reduction.  Proceedings of SIGGRAPH 99, pages 83–90, August 1999.
  * [22] M. C. Vaz.  Journey to Armageddon.  Cinefex, pages 68–93, October 1998.