This work focuses on the long-standing problem of inertial dynamics of an interface with interfacial mass flux and reports new mechanisms for the interface stabilization and destabilization. The interface is a phase boundary separating fluids of different densities and having interfacial mass flux. To analyze the interface dynamics from a far field, we develop and apply the general matrix method to rigorously solve the boundary value problem involving the governing equations in the fluid bulk and the boundary conditions at the interface and at the outside boundaries of the domain. We find the fundamental solutions for the linearized system of equations and analyze the interplay of interface stability with flow fields’ structure by directly linking rigorous mathematical attributes to physical observables. We find that the interface is stable when the dynamics conserves the fluxes of mass, momentum, and energy; the stabilization is due to an inertial mechanism causing small oscillations of the interface velocity. In the classic Landau’s dynamics, the postulate of perfect constancy of the interface velocity leads to the development of Landau–Darrieus instability. This destabilization is also linked to the imbalance of the perturbed energy at the interface. The classic Landau’s solution is found to have degeneracy; lifting of the degeneracy may lead to singularity and self-similar dynamics. Our results compare well with traditional theories of combustion and propose new experiments to study the dynamics of the interface and the flow fields in combustible systems. We further conduct reactive molecular dynamics simulations to elucidate the complexity of chemical processes, to study the destabilizing effect of energy fluctuations on the interface stability, and to illustrate the chemistry-induced instabilities. In summary, we identify the extreme sensitivity of the interface dynamics to the interfacial boundary conditions, including the formal properties of fundamental solutions and the qualitative and quantitative properties of the flow fields. This provides new opportunities for studies, diagnostics, and control of multiphase flows in a broad range of processes in nature and technology.

Interfacial mixing and transport control a broad range of natural and technological processes in fluids, plasmas, and materials from celestial to atomistic scales (Abarzhi and Goddard, 2019; Abarzhi et al., 2013; Zeldovich and Raiser, 2002; and Landau and Lifshitz, 1987). Examples include supernovae and fusion, planetary convection and reactive fluids, formation of phases in liquefied natural gas, and material transformations under impact (Remington et al., 2019; Zhakhovky et al., 2019; Schlimpert et al., 2016; Arjomandnia et al., 2014; Miglani et al., 2014; Abarzhi, 2010; Azechi et al., 2007; Bell et al., 2004a; 2004b; Stein and Norlund, 2000; Rast, 1998; Arnett, 1996; and Prosperetti and Plesset, 1984). In realistic circumstances, material transport is often characterized by sharply and rapidly changing flow fields and by relatively small effects of dissipation and diffusion. This may lead to the formation of interfaces (discontinuities) separating the flow phases (non-uniformities) at macroscopic (continuous) scales (Abarzhi et al., 2019; 2015; Anisimov et al., 2013; Abarzhi, 2010; Kadau et al., 2010; Meshkov, 2006; and Zeldovich and Raizer, 2002). Here, we consider from a far field the evolution of a hydrodynamic discontinuity that separates ideal incompressible fluids of different densities and has the interfacial mass flux (Landau and Lifshitz, 1987). We develop the general matrix method to rigorously solve the linearized boundary value problem and directly link the interface stability to the flow fields’ structure. We find extreme sensitivity of the flow dynamics to the interfacial boundary conditions, identify new mechanisms of the interface stabilization and destabilization, and thus reveal new opportunities for studies of interfacial dynamics in multiphase flows (Abarzhi et al., 2019).

This paper is organized as follows: We start with introduction in Sec. I and present the method in Sec. II, including the governing equations (Sec. II A), linearized dynamics (Sec. II B), and fundamental solutions (Sec. II C). Results in Sec. III present the investigations of the conservative system (Sec. III A), the classic Landau’s system (Sec. III B), and the dynamic Landau’s system (Sec. III C). The properties of interfacial dynamics are discussed in Sec. IV, including stable and unstable solutions (Sec. IV A), the mechanisms of the interface stabilization and destabilization (Sec. IV B), and the effect of energy fluctuations (Sec. IV C). Comparison of our results with traditional studies of combustible systems is given in Sec. V, including the outlines of traditional studies (Sec. V A) and the new experiments (Sec. V B). Chemically reactive systems are analyzed in Sec. VI, including reactive molecular dynamics (Sec. VI A), the RMD2Kin method (Sec. VI B), and the chemistry-induced instability (Sec. VI C). Interface dynamics influenced by acceleration is briefly compared with our results for inertial dynamics in Sec. VII. Multiphase flows outcomes of our theory are outlined in Sec. VIII. Discussion and conclusion are provided in Secs. IX and X, respectively.

In hydrodynamics, the evolution of an interface with mass flow across it is usually studied in premixed combustion and in fluid vaporization (Abarzhi et al., 2019; Zhakhovsky et al., 2019; Prosperretti, 2017; Klimenko and Class, 2000; Class et al., 2003a; 2003b; Peters, 2000; Mayo et al., 1990; Sivashinsky, 1983; 1980; Williams, 1965; and Zeldovich, 1944). Similar problems are considered in plasmas (stability of ablation fronts in inertial confinement fusion), astrophysics (thermonuclear flashes on the surface of compact stars), materials science (melting and evaporation of materials under high pressures and high strain rates), and industrial applications (atomization of liquid jets in scramjets) (Remington et al., 2019; Zhakhovsky et al., 2019; 2014; Hurricane et al., 2014; Anisimov et al., 2013; Abarzhi, 2010; Kadau et al., 2010; Arnett et al., 2009; Gorokhovski and Herrmann, 2008; Azechi et al., 2007; Marmottant and Villermaux, 2004; Bell et al., 2004a; 2004b; Hillebrandt and Niemeyer, 2000; Cherkov and Yakhot, 1998; Piriz et al., 1997; Wu et al., 1997; and Sanz, 1994). The classic theoretical framework for the problem was developed by Landau (1944). It analyzed the dynamics of a discontinuous interface separating ideal incompressible fluids of different densities. By balancing at the interface the fluxes of mass and momentum and by implementing a special condition for the perturbed mass flux, Landau (1944) found that the interface was unconditionally unstable: small perturbations of the interface grew exponentially with time, leading to the development of Landau–Darrieus Instability (LDI) (Darrieus, 1938).

To connect the classical framework (Landau, 1944) to realistic environments, several theoretical approaches have been developed. In high energy density plasmas, significant deviations have been detected between the growth-rates of ablative Rayleigh–Taylor (RT) and accelerated Landau–Darrieus (LD) instabilities (Hurricane et al., 2014; Sanz 1944; and Kull and Anisimov, 1985). In reactive and super-critical fluids, stabilizing influences have been found for dissipation, diffusion, and finite interface thickness on the interfacial dynamics at small scales (Class et al., 2003a; 2003b; Peters, 2010; 2000; Clavin, 2000; Prosperetti and Plesset, 1984; Pelce and Clavin, 1982; Matkowsky and Sivashinsky, 1979; Williams, 1971; Markstein, 1951; and Zeldovich, 1944). More recently, it has been found that the interface stability is sensitive to the flux of energy fluctuations produced by the perturbed interface (Abarzhi et al., 2015). These approaches successfully expanded the classical framework (Landau, 1944; Zeldovich, 1944) to explain the observations (Remington et al., 2019; Hurricane et al., 2014; Aglitsky et al., 2010; Arnett et al., 2009; Azechi et al., 2007; Bell et al., 2004a; 2004b; Class et al., 2003a; 2003b; Peters, 2000; Sivashinsky, 1980; Williams, 1971; Buckmaster, 1979; Markstein, 1951; and Zeldovich, 1944). Some fundamental challenges remain.

First, the general theoretical framework (Landau, 1944) describes the evolution of a phase boundary and is relevant to phenomena far beyond the processes with gradually changing flow fields (Zhakhovsky et al., 2014; Hurricane et al., 2014; Anisimov et al., 2013; Abarzhi, 2010; Kadau et al., 2010; Arnett et al., 2009; Meshkov, 2006; Class et al., 2003a; and 2003b). We have yet to understand whether the interface is stable when the flow quantities experience sharp changes and the effects of dissipation and diffusion are negligible.

Second, direct diagnostics of various physical effects on the dynamics of a multi-phase flow requires detailed information of the interface structure. Such information is often a challenge to directly obtain in experiments and simulations (Remington et al., 2019; Zhakhovsky et al., 2019; Hurricane et al., 2014; Azechi et al., 2007; Bell et al., 2004a; 2004b; Class et al., 2003a; 2003b; Peters, 2000; Clavin, 2000; Rast, 1998; Wu et al., 1997; and Williams, 1971). Besides, one usually observes the flow evolution from a far field and at time scales and length scales that are substantially greater than characteristic scales induced by, e.g., diffusion, dissipation, and finite interface thickness (Remington et al., 2019; Zhakhovsky et al., 2019; Abarzhi et al., 2015; Kadau et al., 2010; Landau and Lifshitz, 1987; and Sivashinsky, 1983). Thus, there is a need to directly connect the interface stability to the flow fields’ structure and to associate the mechanism of (de)stabilization of the interfacial dynamics with qualitative and quantitative properties of the vector and scalar fields of the flow (Abarzhi et al., 2019).

Third, there are still formal questions to clarify in the classic theory of Landau (1944) for ideal incompressible fluids. For instance, Landau’s framework involves consideration of four independent variables: velocity potential of the heavy fluid, velocity potential of the light fluid, vortical potential of the light fluid, and interface perturbation (Landau and Lifshitz, 1987; Landau, 1944). These variables obey the principles of conservation of mass and normal and tangential components of momentum and a special condition for the perturbed mass flow at the interface. However, only two eigenvalues of the linearized dynamics are usually derived (Peters, 2010; Class et al., 2003a; 2003b; Clavin, 2000; Landau and Lifshitz, 1987; Zeldovich, 1994; Williams, 1971; Landau, 1944). These eigenvalues are the unstable positive real eigenvalue (leading to LDI) and its stable negative counterpart. It is unclear whether the other eigenvalues exist and, for each of the fundamental solutions, what are the flow fields, how these flow fields link to interface stability, and what is the effect of fluctuations on the flow stability and fields.

In this work, we address these questions by systematically studying, in a far field approximation, the evolution of the discontinuous interface separating fluids of different densities and having inertial mass flux. We consider a two-dimensional flow of ideal incompressible fluids. By expanding the classic framework (Landau and Lifshitz, 1987), we implement and apply a general matrix method to solve the boundary value problem involving the governing equations in the fluid bulk and the boundary conditions at the interface and at the outside boundaries of the domain. We find formal fundamental solutions of the dynamics, analyze the interplay of interface stability with flow fields’ structure, and directly link rigorous mathematical attributes to physical observables (Abarzhi et al., 2019; 2015). New mechanisms are of the interface stabilization and destabilization are identified. We find that the interfacial dynamics is stable when it conserves the fluxes of mass, momentum, and energy; the stabilization is due to inertial effects, leading to small oscillations of the interface velocity. This mechanism is absent in the classic Landau’s dynamics due to the postulated constancy of the interface velocity, causing LDI to develop. The interface destabilization is also associated with the imbalance of the perturbed energy at the interface, which is in agreement with the classic results (Landau, 1944). We identify the extreme sensitivity of the interface dynamics to the interfacial boundary conditions, including the formal properties of fundamental solutions and the qualitative and quantitative properties of the flow fields.

In this section, we outline the equations governing the interface dynamics, formulate the boundary value problem, and introduce a general matrix method to solve it.

1. Conservation laws and interfacial boundary conditions

The flow dynamics is governed by the conservation of mass, momentum, and energy. For ideal fluids, these are the Euler equations. In an inertial frame of reference,

(1)

where xi are the spatial coordinates with (x1, x2, x3) = (x, y, z), t is time, (ρ, v, P, E) are the fields of density ρ, velocity v, pressure P, and energy density E = ρ(e + v2/2), with the specific internal energy e in its physics definition. The specific internal energy refers to the energy per unit mass that is contained in the system while excluding the kinetic energy of motion of the system as a whole and the potential energy of the system as a whole due to external force fields. The inertial frame of reference refers to a frame of reference moving with a constant velocity V0 = (0, 0, V0). The Euler equations [Eqs. (1)] are augmented with the closure equation—the equation of state associating the internal energy and the pressure (Landau and Lifshitz, 1987).

We introduce a continuous local scalar function θ(x, y, z, t), whose derivatives θ̇ and ∇θ exist (with dot marking a partial time-derivative), such that θ = 0 at the interface, and the heavy (light) fluid is located in the region θ > 0 (θ < 0) (Abarzhi et al., 2015). In the bulk of the heavy (light) fluid, the flow fields are (ρ, v, P, E) → (ρ, v, P, E)h(l), and they are marked with subscript h(l). We represent the flow fields in the domain as (ρ, v, P, E) = (ρ, v, P, E)hH(θ) + (ρ, v, P, E)lH(−θ) by using the Heaviside step-function H(θ) (Abarzhi et al., 2019; Ilyin et al., 2018; and Abarzhi, 2010).

The balances of fluxes of mass and normal and tangential components of momentum and energy across the interface obey the boundary conditions (Landau and Lifshitz, 1987),

(2)

Here, [⋯] = 0 denotes the jump of functions across the interface, n and τ are the unit normal and tangential vectors at the interface with n = ∇θ/|∇θ| and (n · τ)=0, j̃=ρnθ̇/θ+v is the mass flux across the moving interface, and W is the specific enthalpy in its physics definition (Landau and Lifshitz, 1987). The value of specific enthalpy W includes the enthalpy of formation (Landau and Lifshitz, 1987).

Note that the interfacial boundary conditions [Eqs. (2)] are derived from the Euler equations [Eqs. (1)] independently of the equation of state and hold true for ideal fluids with any equation of state, both compressible and incompressible (Landau and Lifshitz, 1987).

Emphasize that in Eqs. (2), since j̃n=0 and j̃=ρnθ̇/θ+v, the condition of continuity of the tangential component of momentum at the interface j̃nj̃τ/ρτ=0 is equivalent to the condition of continuity of the tangential component of velocity at the interface vτ=0 and the absence of shear at the interface, in agreement with the classic results (Landau and Lifshitz, 1987).

2. Derivation of the interfacial boundary conditions

For readers’ convenience, we provide some details on how the interfacial boundary conditions are derived from the governing equations. As a sample case, we consider first of these conditions—the mass conservation condition (Abarzhi et al., 2019; Ilyin et al., 2018).

In the domain, the density of the fluid is represented as ρ = ρhH(θ) + ρlH(−θ), leading to

Here, we account for ∂H(θ)/θ = δ(θ) and δ(θ) = δ(−θ), with δ(θ) being the Dirac delta-function. The other term is calculated as

We substitute these expressions in the equation for mass conservation, introduce the unit normal and tangential vectors at the interface n and τ as

and find

The expression is further reduced to the conditions in the fluids’ bulk and at the interface as in Eqs. (2),

The other interfacial boundary conditions in Eqs. (2) are derived likewise.

Note that the interfacial function θ is applied in numerical modeling of the multiphase flows, including the level set method and the in-phase field method (Osher and Fedkiw, 2002; Anderson et al., 1998; and Sethian, 1996).

There is an important particular case in Eqs. (2), when the mass flux is conserved at the interface, j̃n=0, and when the conserved mass flux is zero at the interface, j̃nθ=0=0. This special condition transforms the boundary conditions in Eqs. (2) to the continuity of the normal component of velocity [v · n] = 0, the continuity of the pressure [P] = 0, and the arbitrariness of the jumps of the tangential component of velocity [v · τ] = arbitrary and enthalpy[W] = arbitrary (Abarzhi et al., 2019; Ilyin et al., 2018). This case corresponds to the dynamics of contact discontinuity (in, e.g., immiscible fluids), such as Rayliegh–Taylor and Richtmyer–Meshkov instabilities (Meshkov and Abarzhi, 2019; Dell et al., 2017; Abarzhi, 2010; Meshkov, 2006; Davies and Taylor, 1950; and Lord Rayleigh, 1883).

3. Inertial reference frame, interface velocity, and outside boundaries

The governing equations and the boundary conditions [Eqs. (1) and (2)] are derived in an inertial frame of reference that moves with a constant velocity V0 (Fig. 1). Some discussion is required on the relation between the velocity of the inertial frame of reference V0 and the velocity Ṽ of the interface as a whole. For a planar and steady interface, the interface velocity Ṽ is constant (Ṽ=Ṽ0), and it may be equal to the velocity V0 of the inertial frame of reference (Ṽ0=V0). In the general case of a non-steady non-planar interface ṼṼ0, the interface velocity Ṽ and the velocity of the inertial frame of reference are distinct quantities (ṼṼ0), and the interface velocity Ṽ obeys the condition Ṽn=vnθ=0=j̃/ρnθ=0 (Peters, 2010; Landau and Lifshitz, 1987).

FIG. 1.

Schematics of the flow dynamics in a far field approximation (not to scale). Blue dashed line denotes the planar interface, and blue solid line denotes the perturbed interface.

FIG. 1.

Schematics of the flow dynamics in a far field approximation (not to scale). Blue dashed line denotes the planar interface, and blue solid line denotes the perturbed interface.

Close modal

In our general theoretical framework, we consider the velocity V0 as the constant velocity of the inertial frame of reference. To study the effect of the interfacial boundary conditions on the interface velocity, we set that the velocity of the steady planar interface equals the velocity of the inertial frame of reference, as Ṽ0=V0.

In the inertial frame of reference, the boundary conditions at the outside boundaries of the domain are

(3)

with constant Vh(l) (Fig. 1). The initial conditions are the initial flow fields at the interface and in the bulk. They define the characteristic length-scale and time-scale of the dynamics (Zeldovich, 1944).

4. The flow dimensionality and the interfacial function

We consider a spatially extended two-dimensional flow propagating in the z direction, being periodic in the x direction, x+λx, and having no motion in the y direction (Fig. 1). We define the local scalar function θ as

(4)

1. Interfacial boundary conditions to the leading and first orders

In Eqs. (1) and (2), the unperturbed flow fields are j̃,v,P,W=J,V,P0,W0, and the normal and tangential unit vectors of the unperturbed interface are {n, τ} = {n0, τ0}. We slightly perturb the flow fields as j̃=J+j, v=V+u, P = P0 + p, and W = W0 + w, with |j| ≪ |J|, |u| ≪ |V|, |p| ≪ |P0|, and |w| ≪ |W0|, respectively. We slightly perturb the fluid interface as n = n0 + n1 and τ = τ0 + τ1, with |n1| ≪ |n0| and |τ1| ≪ |τ0|, and with θ̇/θV. The fluid density is perturbed as ρ → ρ + δρ with |δρ| ≪ |ρ|. In the laboratory frame of reference, the perturbed velocity of the interface is Ṽ=Ṽ0+ṽ, with ṽṼ0.

To the leading order, the boundary conditions at the interface are

(5a)

To the first order, the boundary conditions at the interface are

(5b)

Small perturbations of the flow fields decay away from the interface,

(5c)

The boundary conditions in Eqs. (1) and (2) and (5a)–(5c) are valid for both compressible and incompressible ideal fluids with any equation of state, for two- and three-dimensional flows, and for arbitrary positioning of the interface relative to the mass flux (Landau and Lifshitz, 1987).

2. Directionality, dimensionality, and density variations

The conditions [Eqs. (5a) and (5b)] can be simplified by applying the conditions of the directionality of the mass flux, the dimensionality of the flow, and the density variations in the fluids.

Indeed, to the leading order, the mass flux is J = ρV, the flow fields are uniform in the bulk, (ρ, v, P, W)h(l) = (ρ, V, P0, W0)h(l), and obey conditions [Eqs. (3)] at the outside boundaries of the domain. The components of mass flux normal and tangential to the interface are Jn = J · n0 and Jτ = J · τ0, respectively. We presume that the leading order mass flux is normal to the planar interface; hence, its tangential component is zero. This simplifies the interfacial boundary conditions to the leading order in Eqs. (5a) as

(5d)

For a two-dimensional flow [Eqs. (4)], to the leading order, the normal and tangential vectors of the interface are n0 = (0, 0, −1) and τ0 = (1, 0, 0), respectively. The first order perturbations of the normal and tangential vectors of the interface are n1 = (∂Z*/∂x, 0, 0) and τ1 = (0, 0, ∂Z*/∂x), respectively. This leads to J · n1 = 0 and further simplifies the interfacial boundary conditions to the first order in Eqs. (5b) as

(5e)

For small density variations, |δρ/ρ| ≪ 2|(j · n0)/(J · n0)| and |δρ/ρ| ≪ 2|(J · j)/(J2)|, with |δρ/ρ| ≪ 2|u|/|V| and |δρ/ρ| ≪ |u|/|V|, the expressions in Eqs. (5e) are estimated as

For |δρ/ρ| ≪ |u|/|V|, the first order interfacial boundary conditions [Eqs. (5e)] are further simplified as

(5f)

3. Incompressible dynamics

In this paper, we consider the limiting case of incompressible dynamics in which the speed(s) of sound in the fluid(s) is substantially greater than other velocity scales, and the effect of density variations is negligible, (δρ/ρ) = 0. For ideal gases, the speed of sound is c=γP/ρ, where γ is the fluid’s adiabatic index. For (c/V)h(l) → ∞, the values approach P0+Jn2/ρhlP0hl and W0+J2/2ρ2hlW0hl. Equations (5a) and (5d) are transformed to

(6)

By introducing temperature (T)h(l) and the specific heat at constant pressure (cP)h(l), we present the physics enthalpy as W0hl=W¯0+cPThl. The enthalpy W¯0 is often used in engineering applications, and it has a jump at the interface, with W¯0=Q¯, where Q¯ is the heat release per unit mass. In reactive fluids, the value Q¯ is considered as the specific heat release of a chemical reaction at absolute zero temperature (Landau and Lifshitz, 1987).

Note the difference between the definitions of enthalpy W0 in physics sense accounting for the enthalpy of formation and the enthalpy W¯0 commonly applied in reactive fluids and engineering. Their difference is discussed in detail in Landau and Lifshitz (1987). In full consistency with this discussion, in the incompressible limit, the physics enthalpy W0 is continuous at the interface, [W0] = 0, whereas the enthalpy W¯0 is discontinuous at the interface and its jump at the interface is W¯0=Q¯.

To the first order, the boundary conditions [Eqs. (5f)] at the interface take the form

(7)

where the first order perturbation of the mass flux is j=ρu+n0θ̇. Its normal component is jn=jn0.

In addition to the boundary conditions at the interface, the governing equations in the bulk and at the outside boundaries of the domain to the first order are

(8)

In the laboratory frame of reference, the interface velocity is Ṽ=Ṽ0+ṽ, and the velocity perturbation is ṽ,

(9)

1. Solution structure

We seek solutions for the boundary value problem [Eqs. (7) and (8)] with the structure in which the velocity field of the heavy fluid is potential and the velocity field of the light fluid has both potential and vortical components (Landau and Lifshitz, 1987),

(10a)

The rationale for this structure is the following: In the sub-domain of the heavy fluid, the flow is potential in accordance with the Kelvin theorem. In the sub-domain of the light fluid, the flow can be a superposition of the potential and vortical components. This structure is applied in the theoretical framework (Landau, 1944) and agrees with observations (Peters, 2000; Zeldovich, 1944). It is expected that the solution structure is established for any initial conditions (Landau, 1944; Zeldovich, 1944).

The scalar and vector potentials of the fluid and the interface perturbation are

(10b)

Here, Ω is the growth-rate, i.e., the characteristic frequency or the eigenvalue of the system of linear differential equations, k = 2π/λ is the wavevector, and λ is the wavelength (spatial period). The pressure perturbations ph(l) and the vortex length-scale λ̃=2π/k̃ are identified from Eqs. (8) as follows:

(10c)

The vortical field does not contribute to pressure perturbations (Landau and Lifshitz, 1987). This leads to

(10d)

In the system of Eqs. (7) and (8), the initial conditions set the characteristic length-scale 1/k and the time-scale 1/kVh. We use dimensionless values for the growth-rate ω = Ω/kVh and for the density ratio R = ρhl with R ≥ 1. This leads to Vl/Vh = R and k̃/k=ω/R. We use dimensionless values for the flow fields and the interface as φ = Φ/(Vh/k), φ̃=Φ̃/Vh/k, ψ = Ψ/(Vh/k), z¯=kZ and for the variables as kxx, kzz, kVhtt. In dimensionless units, the fluid potentials are φh=φeix+z+ωt,φl=φ̃eixz+ωt, and ψl = (0, ψl, 0) with ψl=ψeixk̃/kz+ωt, the fluid velocities are uh = ∇φh, ul = ∇φl + ∇ × ψl, and the interface perturbation is z*=z¯eix+ωt.

2. Eigenvalues and eigenvectors of fundamental solutions

With Eqs. (10), the system of differential equations governing the interface dynamics [Eqs. (7) and (8)] is reduced to the linear system Mr = 0, where the matrix M is defined by the boundary conditions at the interface [Eq. (7)]. In dimensionless form, the elements of the matrix M are the functions of the growth-rate (the frequency, the eigenvalue) ω and the density ratio R, and the vector r is r = (φh, φl, z*, ψl)τ. Traditionally, the linear system Mr = 0 is solved by applying the compatibility condition: its eigenvalues are derived, and the (in)stability of the associated solutions is studied (Peters, 2010; Class et al., 2003a; 2003b; Peters, 2000; Clavin, 2000; Landau and Lifshitz, 1987; Williams, 1971; Landau, 1944; Zeldovich, 1944). In this work, we proceed in a more formal way. Particularly, we develop and apply the general matrix method to analyze the interfacial dynamics [Eqs. (1)–(10)]. We find the fundamental solutions for the system of linear ordinary differential equations, including their eigenvalues and eigenvectors, and analyze, for each of these solutions, the stability and the structure of the flow fields. This formal analysis is required to better understand the interfacial dynamics. To our knowledge, it has not been done before.

We derive the eigenvalues ωi by using the condition detΜ(ωi, R) = 0 and identify the associated eigenvectors ẽi by reducing the matrix Μ(ωi, R) to row-echelon form. The matrix Μ is a 4 × 4 matrix. For a non-degenerate 4 × 4 matrix, the expected rank is 4, and the number of eigenvalues ωi and associated eigenvectors ẽi is 4, with i = 1,..., 4, corresponding to four degrees of freedom and four independent variables obeying four equations [Eqs. (7) and (10)].

Solution r for the system given by Eqs. (7) and (10) is a linear combination of fundamental solutions ri,

(11)

with Ci being the integration constants and ri=riωi,ẽi being the fundamental solutions with ri=ẽieωit and ẽi=φeix+z,φ̃eixz,z¯eix,ψeixk̃/kziT, and the associated vector ei=φ,φ̃,z¯,ψiT.

The linear system Μr = 0 with Μ = Μ(ω, R) results from a linear system Pṙ=Sr, where P, S are 4 × 4 matrices, P = P(R), S = S(R), under the assumption that vector r varies in time as reωt. This leads to Μ = (S − ωP). In a non-degenerate case, the inverse P−1 exists, and the system Pṙ=Sr can be reduced to a standard form ṙ=P1Sr. The eigenvalues ωi of the dynamics can then be found from the condition detP1SωI=0, ω = {ωi}, where the unit matrix is I, and index i marks the independent degrees of freedom. Equations det(P−1S − ωI) = 0 and det M(ω, R) = 0 have the same eigenvalues ωi because a linear combination of differential equations preserves the dynamic properties (Beklemishev, 1971).

By using the general matrix method for solving the boundary value problem [Eqs. (1)–(11)], we directly link the microscopic transport at the interface to macroscopic flow fields in the bulk, and we conduct a systematic study of the interplay of the interface stability with the structure of flow fields.

In this section, we find the fundamental solutions for the boundary value problem in the conservative system and the classic and dynamic Landau’s systems. We identify the stability, the flow fields’ structure, the interface velocity, and the formal properties of these fundamental solutions.

1. Fundamental solutions

We consider the properties of dynamics balances of the fluxes of mass and normal and tangential components of momentum and energy across the interface [Eqs. (6)]. We call this dynamics the conservative dynamics (CD). For the conservative dynamics, the matrix Μ has the form Μ = M,

(12a)

Its rank is 4. Its determinant is detM=iR12/RωRω+Rω2+R, and the eigenvalues ωi and eigenvectors ei (derived as described before) are

(12b)

2. Stability of fundamental solutions

For the conservative system [Eqs. (12)], there are four fundamental solutions. Solutions r12ω12,e12 are stable, with ω1=iR, ω2=ω1*, and Reω12=0. Solution r33, e3) is unstable, with ω3 = R and Re[ω3] > 0. Solution r44, e4) is stable, with ω4 = −R and Re[ω4] < 0. Figure 2 shows the dependence of real and imaginary parts of the eigenvalues on the density ratio.

FIG. 2.

Dependence of the eigenvalues on the density ratio for the conservative system.

FIG. 2.

Dependence of the eigenvalues on the density ratio for the conservative system.

Close modal

3. Flow fields of fundamental solutions

By substituting fundamental solutions ri = rii, ei) from Eqs. (12) to Eqs. (10), we find the flow fields for the conservative dynamics. For each of the fundamental solutions, Figs. 3–6 illustrate the interface perturbation z*, the perturbed velocity fields uh(l), the perturbed velocity streamlines sh(l) defined as dshl/dt×uhl=0, the contour plots of the perturbed pressure ph(l), and (when applicable) the vortical field of the light fluid velocity, including the velocity vortical component ∇ × ψl and the contour plot of vorticity ∇ × ul with ×ul=0,k2k̃2ψl,0. The real parts of the functions are shown at some density ratio and some instance of time in the motion plane. In contour plots, red (blue) marks the positive (negative) values. Each plot has its own range of values.

FIG. 3.

Fundamental solution r1 for the conservative system at some density ratio and some instance of time in the motion plane. The solution is the stable traveling wave. Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure.

FIG. 3.

Fundamental solution r1 for the conservative system at some density ratio and some instance of time in the motion plane. The solution is the stable traveling wave. Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure.

Close modal
FIG. 4.

Fundamental solution rCD = (r1 + r2)/2 for the conservative system at some density ratio and some instance of time in the motion plane. The solution is the stable standing wave, which resulted from the interference of the stable traveling waves. Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure.

FIG. 4.

Fundamental solution rCD = (r1 + r2)/2 for the conservative system at some density ratio and some instance of time in the motion plane. The solution is the stable standing wave, which resulted from the interference of the stable traveling waves. Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure.

Close modal
FIG. 5.

Fundamental solution r3 for the conservative system as well as for the classic Landau’s system, dynamic Landau’s system, and the conservative system with energy fluctuations at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation. The solution has zero perturbed fields of the velocity and pressure in both fluids and zero vorticity of the light fluid.

FIG. 5.

Fundamental solution r3 for the conservative system as well as for the classic Landau’s system, dynamic Landau’s system, and the conservative system with energy fluctuations at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation. The solution has zero perturbed fields of the velocity and pressure in both fluids and zero vorticity of the light fluid.

Close modal
FIG. 6.

Fundamental solution r4 for the conservative system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

FIG. 6.

Fundamental solution r4 for the conservative system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

Close modal

Fundamental solution r11, e1) in Eqs. (12) describes stable oscillations of the flow fields in the vicinity of the interface, which decay away from the interface. For this solution, the velocity fields are potential, and all boundary conditions are satisfied, including balances of mass, momentum, and energy at the interface and the conditions at the outside boundaries of the domain. Note that for solution r11, e1), the velocity field of the light fluid ul has a phase shift π/2 relative to that of the heavy fluid uh. That is, the field ul is shifted relative to uh by λ/4 along the x direction of translation [Eqs. (12)]. The velocity fields are also shifted relative to the interface perturbation. These shifts are needed to balance the perturbed mass flux across the perturbed interface in the incompressible approximation (Fig. 3).

For solution r1, the perturbations of pressure ph(l) oscillate in the vicinity of the interface and decay away from the interface. The pressure fields ph(l) are shifted relative to the interface perturbation. The fields of pressure ph and pl are in anti-phase with one another, thus causing stable oscillations in the z direction. These fields span the same range of values with |pl|max(min) = |ph|max(min). For solution r11, e1), the vortical field is zero, ψl = 0, ∇ × ψl = 0, ∇ × ul = 0 (Fig. 3).

The properties of the fundamental solution r22, e2) in Eqs. (12) are similar to those of r11, e1). This is because these solutions are complex conjugates, with ω2=ω1*,e2=e1*,C2=C1* (Fig. 3).

Note that by their very meaning, the solutions r12 are traveling waves. The interference of these traveling waves results in the appearance of standing waves, whose flow fields experience stable oscillations, particularly solutions rCD = (r1 + r2)/2 and r̃CD=r1r2/2. Figure 4 illustrates the flow fields for solution rCD, including the perturbed velocity, the perturbed velocity streamlines, the pressure perturbations, and the interface perturbation at some density ratio and some instance of time in the motion plane.

Fundamental solution r33, e3) obeys all the boundary conditions and is thus a physical solution. Because Re[ω3(R)] > 0 for any density ratio R ≥ 1, it appears at first glance to describe unstable dynamics. However, a more detailed consideration finds that for this solution, with ω3 = R and e3 = (0, i, 0, 1), the velocity fields of the light and heavy fluids are identically zero, with ul = 0 and uh = 0. Furthermore, the interface perturbation and the pressure fields are also zero, z* = 0 and ph(l) = 0 (Fig. 5). Note that while, for this solution, the vortical velocity component is non-zero, ∇ × ψl ≠ 0, its vorticity is zero, ∇ × ul = 0, because for the vorticity field with ×ul=0,k2k̃2ψl,0, the values are k̃/k2=ω3/R2=1 (Fig. 5). These properties hold true in the entire domain, at any time, and for any integration constant C3. Therefore, solution r3 is unstable because its eigenvalue is ω3 > 0, and it corresponds to the unperturbed fields of the velocity, the pressure, and the interface because its eigenvector is e3 = (0, i, 0, 1) (Fig. 5).

Note that since the solution r3 has zero perturbed fields of the velocity, the vorticity, and the pressure, one may set, for this solution, the integration constant equal to zero. The solution r3 may also reflect the symmetry of the dynamics, associated with the arbitrariness of choice of the velocity of the inertial reference frame. Further investigations are required to clarify this issue.

Fundamental solution r44, e4) appears at first glance to be stable and physical (Fig. 6). A more detailed consideration finds that the required choice of the integration constant C4 is C4 = 0 for this solution. Indeed, according to the expression for ω4, e4, at the outside boundaries of the domain, the velocity field of the heavy fluid decays (uh|z→−∞ = 0) at any time. For the velocity field of the light fluid ul, the potential component is vanishing away from the interfaces, ∇φl|z→+∞ = 0, whereas its vortical component ∇ × ψl, while decaying in time, as ∼eRt, increases away from the interface, as ∼ez. Note that while, for this solution, the vortical component of velocity is non-zero, ∇ × ψl ≠ 0, its vorticity is zero, ∇ × ul = 0, because for the vorticity field with ×ul=0,k2k̃2ψl,0 the value is k̃/k2=1. Finally, note that, for this solution, the pressure fields ph(l) are anti-phase with one another and decay away from the interface. Thus, for solution r4 to satisfy boundary condition ulz+=0, we must have C4 = 0 (Fig. 6). Note that due to the exponential growth of the vortical component of the velocity ∇ × ψl away from the interface, ∼ez, the integration constant must remain zero C4 = 0 even when the boundary conditions away from the interface are slightly perturbed (by, e.g., some external noise).

4. Interface velocity for fundamental solutions

For fundamental solutions r1(2), which describe stable traveling waves, as well as for solutions rCD,r̃CD, which describe stable standing waves resulting from interference of the traveling waves r1(2), the interface velocity is Ṽ=Ṽ0+ṽ in the laboratory frame of reference with ṽn0=uhn0+θ̇θ=0. Because of the oscillatory terms uhn0+θ̇e±iRt, the interface velocity experiences stable oscillations near the steady value Ṽ0. The amplitude of the oscillations is small, ṽṼ0, and is set by the initial conditions and the integration constant(s) C1(2), whereas the period of the oscillations 1/R is substantially smaller than the characteristic time-scale of the dynamics.

For fundamental solution r3 in the laboratory frame of reference, the interface velocity is constant, Ṽ=Ṽ0, because its perturbed velocity fields and perturbed interface are zero for any integration constant C3. For solution r4, in the laboratory frame of reference, the interface velocity is constant Ṽ=Ṽ0 because its integration constant is zero, C4 = 0.

5. Formal properties of the conservative system

To study the formal properties of the conservative system, we find, for matrix M = M, the associated matrices S = SM and P = PM,

(12c)

The solutions of equations det(PM−1SM − ωI) = 0 and detM = 0 yield the same eigenvalues: ω12=±iR and ω3 = R, ω4 = −R. The conservative dynamics has four eigenvalues, four fundamental solutions, and four independent degrees of freedom. It is non-degenerate.

6. Summary of properties of solutions for the conservative system

For the conservative system, there are four fundamental solutions rii, ei), i = 1, 2, 3, 4 [Eqs. (12), Figs. 2–6, and Table I]. Solutions r1 and r2 with integration constant C1,C2=C1* describe the stable oscillations of the flow fields. These solutions are stable traveling waves; their interference results in stably oscillating standing waves rCD and r̃CD. For these solutions, the interface velocity Ṽ experiences small stable oscillations near the steady value Ṽ0.

TABLE I.

Attributes of fundamental solutions for the conservative system (1, 2, 3, 4) with the asterisk marking a quantity.

r1r2r3r4
Re[ω] >0 <0 
Im[ω] >0 <0 
C *(0) 
uh 
ul 
∇×ul 
z* 
r1r2r3r4
Re[ω] >0 <0 
Im[ω] >0 <0 
C *(0) 
uh 
ul 
∇×ul 
z* 

Solution r3, corresponds to the unperturbed fields of the velocity, the pressure, and the interface at any constant C3. For solution r4, the value C4 must be zero to obey all the boundary conditions, as discussed in the foregoing, C4 = 0. For solutions r3(4), the interface velocity is constant, Ṽ=Ṽ0.

1. Fundamental solutions

Landau (1944) analyzed the stability of the interface by considering the boundary conditions balancing the transport of mass and normal and tangential components of momentum at the interface and by omitting the equation for the energy balance at the interface (Peters, 2010; Landau, 1944; Zeldovich, 1944; and Mallard and Le Chatelier, 1883). Since four degrees of freedom and four independent variables require four conditions, Landau employed an additional condition for the continuity of the normal component of the perturbed velocity at the interface un0=0 or jnθ=0=0 at the interface (Landau, 1944).

With this condition, the boundary conditions have the form

(13)

Recall that the perturbed mass flux is j=ρu+n0θ̇ with jn=jn0=ρun0+θ̇. The condition jn|θ=0 = 0 leads to un0=θ̇ for any ρ and thus to un0=0 in Eqs. (13). The condition for the normal component of momentum is transformed to [pn0] = 0, and the field of pressure perturbation is continuous at the interface (Landau, 1944).

For the classic Landau’s system, the matrix Μ is Μ = L,

(14a)

Its rank is 4. Its determinant is detL=iR1/RωRR+1ω2+2RωRR1. Its eigenvalues ωi and eigenvectors ei are

(14b)

2. Stability of fundamental solutions

For the classic Landau’s dynamics, there are three fundamental solutions. Solution r11, e1) is unstable, Re[ω1] > 0. Solution r22, e2) is stable, Re[ω2] < 0. Solution r33, e3) is unstable, Re[ω3] > 0, and is identical to that for the conservative system. Figure 7 illustrates the dependence of eigenvalues on the density ratio.

FIG. 7.

Dependence of the eigenvalues on the density ratio for the classic and the dynamic Landau systems. The classic Landau system has three eigenvalues (1, 2, 3). The dynamic Landau system has four eigenvalues: the first three are identical to the classic system eigenvalues; the fourth behaves as shown.

FIG. 7.

Dependence of the eigenvalues on the density ratio for the classic and the dynamic Landau systems. The classic Landau system has three eigenvalues (1, 2, 3). The dynamic Landau system has four eigenvalues: the first three are identical to the classic system eigenvalues; the fourth behaves as shown.

Close modal

3. Flow fields of fundamental solutions

By substituting solutions ri = rii, ei) in Eqs. (14) to Eqs. (10) for each i, we find the flow fields for fundamental solutions describing the classic Landau’s dynamics.

For each of fundamental solutions rii, ei) for the classic Landau’s dynamics, Figs. 8 and 9 illustrate the interface perturbation z*, the perturbed velocity fields uh(l), the perturbed velocity streamlines sh(l) defined as dshl/dt×uhl=0, the contour plots of the perturbed pressure fields ph(l), and the vortical field of the light fluid velocity, including the vortical velocity component ∇ × ψl and the contour plot of vorticity ∇ × ul with ×ul=0,k2k̃2ψl,0. The real parts of the functions are shown at some density ratio and some instance of time in the motion plane. In contour plots, red (blue) marks the positive (negative) values. Each plot has its own range of values.

FIG. 8.

Fundamental solution rLD = r1 for the classic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

FIG. 8.

Fundamental solution rLD = r1 for the classic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

Close modal
FIG. 9.

Fundamental solution r2 for the classic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

FIG. 9.

Fundamental solution r2 for the classic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

Close modal

Fundamental solution r11, e1)in Eqs. (14) corresponds to the Landau–Darrieus instability (Landau, 1944) (Fig. 8). This solution, although unstable, is a physical solution. It satisfies all boundary conditions, including the assigned boundary conditions at the interface and the conditions at the outside boundaries of the domain [Eqs. (3) and (13)]. This solution describes dynamics in which the interface perturbation and the vortical and potential components are strongly coupled. The perturbed velocity fields uh(l) and the perturbed velocity streamlines sh(l) illustrate the formation of vortical structures near the interface and in the bulk of the light fluid. For solution r11, e1), the vortical component, ∇ × ψl, and the vorticity ∇ × ul of the velocity of the light fluid, while increasing in time, decay away from the interface. The fields of pressure perturbations ph and pl are equal to one another at the interface and decay away from the interface. They are in phase with one another and in anti-phase with the interface perturbation. The fields of pressure ph(l) are symmetric and span the same range of values with |pl|max(min) = |ph|max(min) (Fig. 8).

The vortical field for solution r11, e1) in Eqs. (14) has the following properties (Fig. 8): (1) Consistent with expressions in Eqs. (14), the potential components of uh and ul are in anti-phase with one another, whereas the vortical component of velocity ul is shifted relative to the potential components of velocities uh and ul and the interface perturbation z*. (2) For the length-scale of the vortical field, we have k̃/k=ω1/R, k̃/k=R+R+R2+R3/R1+R. The length-scale of the vortical field is large, with k̃/kR1/20 for R → 1+ and with k̃/kR1/20 for R → ∞. The maximum value of k̃/k (minimum value of λ̃/λ) is achieved at R=2+54.24. (3) For R → 1+, the values are ψ/φ̃R1 and ψ/φ̃R1, whereas for R → ∞, the values are ψ/φ̃R3/2 and ψ/φ̃1. Therefore, the contribution of the vortical field to the dynamics is quantitatively large for fluids with very different densities and is quantitatively small for fluids with similar densities. (4) To quantify the strength of the vorticity field for given ψ, we note that ×ul/ψk2=1k̃/k2. For solution r11, e1), this value is ∼1 for all R ≥ 1.

The other fundamental solution in Eqs. (14) is solution r22, e2) (Fig. 9). Solution r2 describes the dynamics for which the interface perturbation and the vortical and potential components are strongly coupled. It appears stable and physical at first glance. Its fields of pressure ph(l) are equal to one another and are in-phase at the interface and decay away from the interface. For this solution, the field of the perturbed velocity of the heavy fluid vanishes at the outside boundaries of the domain uh|z→−∞ = 0 at all times. The potential component of the perturbed velocity of the light fluid ul also vanishes away from the interface, ∇φl|z→+∞ = 0. Its vortical component, ∇ × ψl, decays exponentially in time but increases far from the interface. The strength of the vorticity field ∇ × ul for this solution also increases far from the interface. For fundamental solution r22, e2) in Eqs. (14) to satisfy the boundary condition ul|z→+∞ = 0 at any time, the value of integration constant C2 should be zero (C2 = 0). This integration constant must remain zero C3 = 0 even when the boundary conditions away from the interface are slightly perturbed (due to the exponential growth of the vortical component of the velocity).

Fundamental solution r33, e3) in Eqs. (14) is identical to that in Eqs. (12) (Fig. 5). Note that solution r3 exists in the classic framework (Landau and Lifshitz, 1987; Landau, 1944) and can be found similarly to Eqs. (13) and (14).

4. Interface velocity for fundamental solutions

For fundamental solution r1, in the laboratory frame of reference, the interface velocity is Ṽ=Ṽ0+ṽ with ṽn0=uhn0+θ̇θ=0. While in this expression, each of the terms uh and θ̇ grows exponentially in time, uh,θ̇eω1t, the exact balance in the boundary conditions un0=0 in Eqs. (13) leads to uhn0+θ̇θ=0=0 and ṽ=0 and thus to the constancy of the interface velocity Ṽ=Ṽ0. The classic Landau’s dynamics is the perfect mathematical match.

For fundamental solutions r23, in the laboratory frame of reference, the interface velocity is constant Ṽ=Ṽ0, similarly to the foregoing.

5. Formal properties of the classic Landau’s system

To study the formal properties of the classic Landau’s system, we find, for matrix Μ = L, the associated matrices S = SL and P = PL,

(14c)

In matrix PL, the fourth row only has null elements and detPL = 0. Thus, the inverse matrix PL1 does not exist. This suggests the degeneracy of the classical Landau’s system and a singular character of the Landau–Darrieus instability.

6. Summary of properties of solutions for the classic Landau’s system

For the classic Landau’s system, there are three fundamental solutions rii, ei), i = 1, 2, 3 [Eqs. (13) and (14), Figs. 5 and 7–9, and Table II]. Solution r11, e1) with C1 ≠ 0 describes the Landau–Darrieus instability. For solution r22, e2), the integration constant must be zero, C2 = 0. Solution r33, e3) corresponds to the unperturbed fields of velocity and pressure at any C3. For each of these solutions, in the laboratory frame of reference, the interface velocity is constant Ṽ=Ṽ0.

TABLE II.

Attributes of fundamental solutions for the classic Landau’s system (1, 2, 3) and the dynamic Landau’s system (1, 2, 3, 4) with the asterisk marking a quantity.

r1r2r3r4o
Re[ω] >0 <0 >0 
Im[ω] 
C *(0) 0(*) 
uh 
ul 
∇×ul 
z* 
r1r2r3r4o
Re[ω] >0 <0 >0 
Im[ω] 
C *(0) 0(*) 
uh 
ul 
∇×ul 
z* 

1. Fundamental solutions

The paradox of the classic Landau’s system (which to our knowledge has not been discussed in detail before) is that it has a smaller than expected number of fundamental solutions. In physics sense, this paradox can be understood from the following considerations: The boundary conditions in the classic Landau’s system [Eqs. (13)] imply dynamic conditions for the conservation of mass and normal and tangential components of momentum at the interface because they contain partial time-derivatives of components of vector r. Landau’s special condition in Eqs. (13) is a “static” condition of continuity of the normal component of the perturbed velocity at the interface [u · n0] = 0. This condition implies that the perturbed mass flux is a constant value at the interface, jn|θ=0 = const., and furthermore, this value is zero, jn|θ=0 = 0. In order to resolve the paradox, we modify the special Landau’s condition jn|θ=0 = 0 to a more general form jn|θ=0 = const. and implement it in a dynamic condition, un0/t=0, at the interface. While other implementations are possible, we note that the condition of constant mass flux at the interface is met in realistic environments (Azechi et al., 2007; Peters, 2000).

For this dynamic implementation of Landau’s system, the modified boundary conditions are

(15)

For the dynamic Landau’s system, the matrix Μ has the form M=L̃,

(16a)

Its rank is 4. Its determinant is detL̃=iR1/RωωRR+1ω2+2RωRR1. It has four eigenvalues ωi and four eigenvectors ei. For i = 1, 2, 3, the eigenvalues ωi and eigenvectors ei are the same as in Eqs. (14). For i = 4, the fundamental solution is r4oω4o,e4o with

(16b)

Symbol “o” in r4oω4o,e4o emphasizes the presence of the null eigenvalue.

2. Stability of fundamental solution

Fundamental solutions rii, ei) with i = 1, 2, 3 of Eqs. (16) are the same as solutions for system [Eqs. (14)] and have the same stability properties. Fundamental solution r4oω4o,e4o is neutrally stable because of the null eigenvalue ω4o=0 (Fig. 7).

3. Flow fields of fundamental solutions

Fundamental solutions rii, ei) with i = 1, 2, 3 in Eqs. (16) are the same as solutions for system [Eqs. (14)] and have the same corresponding fields of the velocity, the pressure, and the interface perturbation (Figs. 5, 8, and 9). For fundamental solution r4oω4o,e4o, Fig. 10 presents the fields of the perturbed velocity uh(l), the perturbed streamlines sh(l), the perturbed pressure ph(l), the perturbed vortical field ∇ × ψl, the vorticity ×ul=0,k2k̃2ψl,0, and the interface perturbation z* at some instance of time for some density ratio in the motion plane. Each plot has its own range of values.

FIG. 10.

Fundamental solution r4o, rDLD=r4o, for the dynamic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

FIG. 10.

Fundamental solution r4o, rDLD=r4o, for the dynamic Landau’s system at some density ratio and some instance of time in the motion plane. (a) Plots of the interface perturbation, the perturbed velocity vector fields, the perturbed velocity streamlines, and the perturbed pressure. (b) Plots of the perturbed velocity vortical component, the vorticity, and the interface perturbation.

Close modal

Solution r4oω4o,e4o is independent from fundamental solutions rii, ei) with i = 1, 2, 3 in Eqs. (14) and has the following properties: (1) The vortical component of velocity ul is shifted relative to the potential components of velocities uh, ul and the interface perturbation z*, whereas the potential components of uh and ul are in anti-phase with one another. (2) The vortical field length-scale is infinitely large, with k̃/k=0 and λ̃/λ=, due to the null eigenvalue ω4o=0. (3) The pressure fields of ph and pl are in phase and decay away from the interface. The pressure fields are no longer symmetric and |pl|max(min) ≪ |ph|max(min). (4) The vortical component of the velocity is ∇ × ψl ≠ 0, and the vorticity is ∇ × ul ≠ 0. They change periodically in the x direction and are uniform in the z direction. (5) The values are |ψ/φ| ∼ R and ψ/φ̃2 for R → ∞ and ψ/φψ/φ̃R1 for R → 1+, and |∇ × ul|/|ψ|k2 = 1 for any R.

For solution r4oω4o,e4o, the velocity uh vanishes at z → −∞, whereas for the velocity ul at z → +∞, the potential component decays and the vortical component changes periodically. The boundary conditions ul|z→+∞ = 0 require us to set, for this solution, the integration constant equal to zero. Note that for solution r4oω4o,e4o, some slight perturbations of the boundary conditions away from the interface z → +∞ (due to some external noise) may lead to non-zero integration constant C4.

4. Interface velocity for fundamental solutions

For fundamental solutions ri with i = 1, 2, 3 in Eqs. (16), the interface velocity in the laboratory frame of reference is constant, Ṽ=Ṽ0, as discussed in the foregoing.

For fundamental solutions r4o in Eqs. (16) in the laboratory frame of reference, the interface velocity is Ṽ=Ṽ0+ṽ with ṽn0=uhn0+θ̇θ=0. Due to the zero eigenvalue, ω4o=0, the value is θ̇=0, whereas the value is uhn0θ=0=const. Thus, for solution r4o with a constant perturbed mass flux across the interface, un0/t=0, for non-zero integration constant C4, the interface velocity Ṽ=Ṽ0+ṽ is slightly shifted from its steady value Ṽ0, ṽṼ0, and the shift is stationary, ṽ̇=0. For zero integration constant C4 = 0, the interface velocity is constant.

5. Formal properties of the dynamic Landau’s system

a. Non-degeneracy of the dynamic Landau’s system.

To study the formal properties of the dynamic Landau’s system, we find for matrix M=L̃ the associated matrices S=SL̃ and P=PL̃,

(16c)

The equations detPL̃1SL̃ωI=0 and detL̃=0 yield the same eigenvalues. This implies that the dynamic Landau’s system is non-degenerate. It has four fundamental solutions for four degrees of freedom, thus resolving the long-standing formal paradox.

b. Singular character of the classic Landau’s system.

An important consequence of the zero eigenvalue ω4o=0 is the singular nature of the LDI and the classic Landau’s solution (Kadanoff et al., 1967; Barenblatt et al., 1962). Indeed, differential equations [Eqs. (13)–(16)] govern the dynamics of slight perturbations r near the equilibrium state of the uniform flow fields in Eqs. (6) and (7)–(10) for which r = 0. It is presumed that these deviations are exponential in time, with r varying as reωt, and with characteristic time-scale being ∼|ω−1|. For ω = 0, formal substitution leads to retconst. This implies that slight deviations from the equilibrium dynamics are actually power-law functions of time, rta, rather than exponential functions, reωt (Kadanoff et al., 1967; Barenblatt et al., 1962). The requirement of the asymptotic stability of fundamental solution r4oω4o,e4o leads to condition Re[a] < 0 (Landau and Lifshitz, 1987). For Re[a] < 0, slight initial deviations from the equilibrium, rta, vanish asymptotically with time, r → 0 for t → ∞. To derive the exponent a, one has to account for nonlinear terms. We address this problem to future research.

Note that the power-law functions are scale-invariant and describe processes with no characteristic scales. The power-law dependence rta may also imply that the system may have a singularity and experience a phase transition (Kadanoff, 2000; Kadanoff et al., 1967). Note also that with ω4o=0, solution r4oω4o,e4o in Eqs. (16) can be viewed as a steady vortical field that is imposed at some instance of time (by some external noise) and is neutrally stable (Figs. 8 and 10). This vortical field can “seed” the vortical field in Landau’s solution r11, e1) in Eqs. (14) and trigger the unstable dynamics.

6. Summary of properties of fundamental solutions for the dynamic Landau’s system

For the dynamic Landau’s system in Eqs. (15) and (16), there are four fundamental solutions rii, ei) with i = 1, ..., 4 [Eqs. (14) and (16), Table II, and Figs. 5 and 7–10]. Solution r11, e1) with C1 ≠ 0 describes the Landau–Darrieus instability. For solution r22, e2), the integration constant is C2 = 0. Solution r33, e3) corresponds to the unperturbed fields of the velocity, the pressure, and the interface at any C3. Solution r4oω4o,e4o is neutrally stable and may lead to a power-law dynamics near the equilibrium seeding LDI.

For solutions ri with i = 1, 2, 3, the interface velocity in the laboratory frame of reference is constant, Ṽ=Ṽ0, as discussed in the foregoing. For solution r4owith C4 ≠ 0, the interface velocity Ṽ=Ṽ0+ṽ is slightly shifted from its steady value Ṽ0, ṽṼ0, and the shift is stationary, ṽ̇=0. For zero integration constant C4 = 0, the interface velocity is constant.

In this section, we focus on the physical properties of the interfacial dynamics that are prescribed by fundamental solutions for the conservative system and the classic and dynamic Landau’s systems. We further find the new mechanisms of the interface stabilization due to the inertial effects, and we model the influence of energy fluctuations on the interface dynamics.

1. Stable, neutrally stable, and unstable standing waves

We compare the properties of fundamental solution rCD for the conservative system, fundamental solution rLD for the classic Landau’s system, and fundamental solution rDLD for the dynamic Landau’s system, where rCD = (r1 + r2)/2 in Eqs. (12), rLD = r1 in Eqs. (14), and rDLD=r4o in Eqs. (16). Solutions rCD, rLD, and rDLD are standing waves; they have distinct properties (Tables III–V).

TABLE III.

Eigenvalues of standing waves.

rCD ωCD=±iR 
rLD ωLD=R+R3+R2RR+1 
rDLD ωDLD=0 
rCD ωCD=±iR 
rLD ωLD=R+R3+R2RR+1 
rDLD ωDLD=0 
TABLE IV.

Properties of the conservative dynamics (CD) and the classic Landau–Darrieus (LD) dynamics.

CDLD
Conservation properties Conserves mass, momentum, Conserves mass and momentum 
and energy at the interface and has zero perturbed mass flux at the interface 
Interface velocity Slight stable oscillations near Constant value (postulate) 
the constant value 
Flow field Potential velocity fields Vortical structures in the light fluid 
Interfacial shear Shear-free Shear-free 
Formal properties Non-degenerate; four fundamental Degenerate; three fundamental 
solutions and four solutions and four degrees of 
degrees of freedom freedom 
Stability Stable; stabilized by inertial Unstable 
effects 
CDLD
Conservation properties Conserves mass, momentum, Conserves mass and momentum 
and energy at the interface and has zero perturbed mass flux at the interface 
Interface velocity Slight stable oscillations near Constant value (postulate) 
the constant value 
Flow field Potential velocity fields Vortical structures in the light fluid 
Interfacial shear Shear-free Shear-free 
Formal properties Non-degenerate; four fundamental Degenerate; three fundamental 
solutions and four solutions and four degrees of 
degrees of freedom freedom 
Stability Stable; stabilized by inertial Unstable 
effects 
TABLE V.

Properties of the classic Landau–Darrieus (LD) and dynamic Landau–Darrieus (DLD) systems.

LDDLD
Conservation properties Conserves mass and momentum Conserves mass and momentum 
and has the null perturbed and has a constant perturbed 
interfacial mass flux mass flux at the interface 
Interface velocity Constant value (postulate) May be slightly shifted near the 
constant value 
Flow field The vortical field length-scale The vortical field length-scale is 
depends on the density ratio and independent of the density ratio 
is large and finite when compared and is infinite when compared to 
to the perturbation wavelength the perturbation wavelength 
Interfacial shear Shear-free Shear-free 
Formal properties Degenerate; three fundamental Non-degenerate; four fundamental 
solutions and four degrees of solutions and four degrees 
freedom of freedom 
Stability Unstable Neutrally stable 
LDDLD
Conservation properties Conserves mass and momentum Conserves mass and momentum 
and has the null perturbed and has a constant perturbed 
interfacial mass flux mass flux at the interface 
Interface velocity Constant value (postulate) May be slightly shifted near the 
constant value 
Flow field The vortical field length-scale The vortical field length-scale is 
depends on the density ratio and independent of the density ratio 
is large and finite when compared and is infinite when compared to 
to the perturbation wavelength the perturbation wavelength 
Interfacial shear Shear-free Shear-free 
Formal properties Degenerate; three fundamental Non-degenerate; four fundamental 
solutions and four degrees of solutions and four degrees 
freedom of freedom 
Stability Unstable Neutrally stable 

The dynamics rCD is stable, the dynamics rLD is unstable, and the dynamics rDLD is neutrally stable. For solution rCD, the velocity fields are potential; for solutions rLD and rDLD, the potential and vortical components of the velocities and the interface perturbation are strongly coupled. There is zero shear at the interface for each of the dynamics rCD, rLD, and rDLD. Solution rCD conserves mass, momentum, and energy at the interface. Solution rLDG conserves mass and momentum and has zero perturbed mass flux at the interface. Solution rDLD conserves mass and momentum and has a constant perturbed mass flux at the interface. The dynamics rCD is non-degenerate (four solutions and four degrees of freedom). The dynamics rLD is degenerate (three solutions and four degrees of freedom). The dynamics rDLD is non-degenerate (four solutions and four degrees of freedom) (Tables III–V).

For solutions rDLD and rLD, the potential and vortical components of the velocities are strongly coupled with the interface perturbations. Yet, for solution rLD, the length-scale of the vortical field λ̃=2π/k̃ depends on the density ratio, whereas for solution rDLD, the length-scale of the vortical field is independent of the density ratio and is infinitely large (λ̃=) (Tables III–V).

2. Other unstable solutions

The conservative system, the classic Landau system, and the dynamic Landau system have a fundamental solution in common, r33, e3) in Eqs. (12), (14), and (16), which can be implemented for any integration constant C3. This solution is the fastest unstable solution, as suggested by its eigenvalue ω3. This solution has the null perturbation fields of the velocity, pressure, and interface, as suggested by its eigenvector e3 (Fig. 5).

Note that this solution is in agreement with the classic works (Landau and Lifshitz, 1987) and can be derived from the governing equations in the framework (Landau, 1944). The existence of this solution may imply that the velocity of the inertial frame of reference V0 is a free quantity, and it may differ from the velocity of the planar steady interface Ṽ0. It may also imply that the interfacial dynamics may have flow fields that are uniform and are indistinguishable from those of the equilibrium state. While the solution is formally unstable, its dynamics may appear as being “fully stabilized.” This apparent stabilization is achieved in the case of ideal incompressible fluids and at any density ratio. In realistic fluids, various physical effects may influence the null perturbation fields of the velocity, pressure, and interface associated with this solution, to be studied in the future research.

3. Other stable solutions

Our analysis identifies some new properties of stable fundamental solutions, which are the stable solutions that decay exponentially with time due to their real negative eigenvalues. These include solution r4 in Eqs. (12) for the conservative system and solution r2 in Eqs. (14) and (16) for the classic and dynamic Landau’s systems. As discussed in the foregoing, these stable solutions require a null integration constant in order to obey the boundary conditions at the outside boundaries of the domain at any times.

A somewhat similar requirement was discussed in the study of the LDI by means of the Laplace transform (Zeldovich, 1944). The study interpreted a non-zero vortical component of the velocity field away from the interface as vortical structures that are produced at the interface and are moved to the bulk. The study required a null integration constant for the corresponding solution because these vortical structures have no influence on the development of LDI (Zeldovich, 1944). In our analysis, we require the null integration constant for the exponentially decaying stable solutions to obey the boundary conditions at the outside boundaries of the domain at any time for any density ratio. We further note that for stable solutions decaying exponentially with time, the requirement of zero integration constant directly follows from the solution structure. Thus, it should hold in realistic fluids, especially in the case of fluids with contrasting densities, when the contribution of the vortical field to the dynamics is large.

1. Interface stabilization for the conservative dynamics

To better understand the mechanism of the interface stabilization of the conservative dynamics, we consider the interface velocity. In the laboratory reference frame, the interface velocity is Ṽ=Ṽ0+ṽ, with ṽn0=un0+θ̇θ=0.

For the conservative dynamics rCD, the values are un0θ=0e±iRt,θ̇θ=0e±iRt, leading to un0+θ̇θ=0e±iRt and ṽn0e±iRt. Thus, the interface velocity experiences small stable oscillations near the steady value Ṽ0,ṼṼ0n0e±iRt, as discussed in the foregoing.

This suggests the inertial effect as the stabilization mechanism of the conservative dynamics rCD. When the interface is slightly perturbed, the parcels of the heavy fluid and the light fluid follow the interface perturbation causing the change of momentum and energy of the fluid system. Yet, the dynamics is inertial. To conserve the momentum and energy, the interface should slightly change its velocity. This causes the reactive force occurs and to stabilize the dynamics. The resulting flow in the conservative dynamics is a superposition of two motions: the background motion of the fluids following the interface, whose velocity slightly oscillates near the steady value, and the stable oscillations of the interface perturbations (Fig. 11).

FIG. 11.

Schematics of the flow dynamics in a far field approximation (not to scale) for the conservative dynamics (left) and the classic Landau’s dynamics (right) in the inertial reference frame. Blue dashed line denotes the planar interface, and blue solid line denotes the perturbed interface (solid). For the conservative dynamics, the blue double arrows denote the oscillations of the interface perturbations (solid) and the interface velocity as a whole (dashed) with the latter occurring due to inertial effects and causing the reactive force to occur. For the classic Landau’s dynamics, the single blue arrows denote the growth of the interface perturbations; the velocity of the interface is the postulated constant.

FIG. 11.

Schematics of the flow dynamics in a far field approximation (not to scale) for the conservative dynamics (left) and the classic Landau’s dynamics (right) in the inertial reference frame. Blue dashed line denotes the planar interface, and blue solid line denotes the perturbed interface (solid). For the conservative dynamics, the blue double arrows denote the oscillations of the interface perturbations (solid) and the interface velocity as a whole (dashed) with the latter occurring due to inertial effects and causing the reactive force to occur. For the classic Landau’s dynamics, the single blue arrows denote the growth of the interface perturbations; the velocity of the interface is the postulated constant.

Close modal

2. Interface destabilization for the classic Landau’s dynamics

For the classic Landau–Darrieus dynamics rLD, the interface velocity is constant, Ṽ=Ṽ0, as it is postulated by the boundary condition un0=0, leading to un0+θ̇θ=0=0,ṽ=0. The constancy of the interface velocity prevents the reactive force to occur. The dynamics is destabilized. The LD unstable flow is a superposition of two motions—the background motion of the fluids following the interface with the constant velocity and the growth of the interface perturbations (Fig. 11).

Since, in the classic LD dynamics, mass and momentum are conserved, for capturing the destabilization mechanism, we consider how the condition un0=0 may influence the energy transport. Remarkably, for ideal incompressible fluids, the solution rLD is incompatible with the condition for energy balance at the perturbed interface (Abarzhi et al., 2019). Indeed, by substituting un=0jn=0 in the condition for the perturbed energy balance Jnw+Jj/ρ2=0, one obtains Jnw+Jj/ρ2=Jnw=0 and, with [Jn] = 0, reduces it further to [w] = 0 (Abarzhi et al., 2015).

The enthalpy perturbations are w = p/ρ + δe, where δe are the fluctuations of internal energy (in physics sense). In ideal incompressible fluids, free from energy sources, these fluctuations are zero, δe = 0, because ė=0,e=0. Thus, with w = p/ρ and with ρh ≠ ρl, the condition for energy balance, [w] = 0, contradicts the condition for momentum balance, [p] = 0. We see that for ideal incompressible fluids, the classic Landau’s dynamics, while conserving the energy (and the physics enthalpy) to the zero order, [W0] = 0, requires energy imbalance at the interface to the first order. This imbalance is induced by the first order enthalpy perturbations, [w] = [p/ρ]. This is the work done by the fluid when a parcel of fluid of a unit mass expands its volume from ρh1 to ρl1 under pressure p.

In realistic fluids, this imbalance of the perturbed energy can be induced by energy fluctuations. The effect can be self-consistently derived from entropy conditions accounting for chemical reactions. In ideal fluids, to evaluate the effect of the imbalance of perturbed energy on the interface stability, we may introduce an artificial energy flux and study a transition from stable to unstable dynamics with an increase in its strength (Abarzhi et al., 2015).

The theoretical framework of Landau (1944) describes the evolution of a phase boundary: On each side of the interface, the fluids are ideal and incompressible. They are transformed into one another at the interface. The process is accompanied by a mass flow across the interface. At the planar interface, the fluxes of mass, momentum, and energy are fully balanced. For the perturbed interface, in some circumstances, the interface becomes unstable. The instability is accompanied by the appearance of a vortical field and requires the constancy of the interface velocity, the complete balance of the perturbed momentum and, for ideal incompressible fluids, an imbalance of perturbed energy, or energy fluctuations by the perturbed interface. Here, we model the effect of energy fluctuations on the interface stability and the flow fields’ structure (Abarzhi et al., 2015).

Mathematically, it is straightforward to verify that adding a stationary energy source at the moving interface in the governing equations, such as j̃nW+j̃2/2ρ2Q̃=0, keeps the boundary conditions and the interfacial dynamics the same, except for the modification of the zero-order enthalpy as W0W̃0, assuming that the fluids are ideal and incompressible. To influence the interface stability, the energy imbalance, or energy fluctuations, should be produced by the perturbed interface.

For realistic fluids, the flux of energy fluctuations occurs due to the difference in the thermodynamic properties of the fluids (Kadanoff, 2000; Landau and Lifshitz, 1987). It should be derived self-consistently from the condition for entropy conservation (Landau and Lifshitz, 1987). For ideal incompressible fluids, we model the effect of energy fluctuation by introducing an artificial energy flux produced by the perturbed interface. We modify the boundary conditions as (Abarzhi et al., 2015)

(17)

This model [Eqs. (17)] can be used only for estimates of the effect of the perturbed energy imbalance, or the energy fluctuations, on the interface stability. A self-consistent study is addressed to the future.

1. Fundamental solutions

For the conservative system with energy fluctuations, with Qhl=qhlVh2 in Eqs. (17), the matrix Μ is M=M̃,

(18a)

Its rank is 4. It has four eigenvalues ωi and four corresponding eigenvectors ei. With qhql = q, its determinant is detM̃=iR12ωRω+Rω2+R+iqωRR+1ω2+2RωRR1. The parameter q, q > 0, describes the strength of energy fluctuations (the perturbed energy imbalance).

With detM̃=RdetMqR/R1detL and detM̃=RdetMq/ωR/R1detL, the condition detM̃=0 leads to (R − 1)detM = qdetL and R1detM=q/ωdetL̃. Similarly to Abarzhi et al. (2015), we scale the fluctuations strength as q=fR12R/R+1, where f ≥ 0 is a constant, and reduce equation detM̃=0 to

(18b)

The dependence of fundamental solutions ri = rii, ei) with i = 1, 2, 4 on the density ratio R and the fluctuation strength f is cumbersome and not presented here. Solution r3 = r33, e3) with ω3 = R and e3 = (0, i, 0, 1) is the same as for the conservative system, the classic Landau’s system, and the dynamic Landau’s system. Figure 12 illustrates the dependence of the eigenvalues on the fluctuation strength and the dependence of the eigenvalues on the density ratio for weak and strong fluctuations.

FIG. 12.

Dependence of the eigenvalues on the fluctuations strength for the conservative system with energy fluctuations (a) at some density ratio; (b) top—for a broad range of fluctuations’ strength at some density ratio; bottom—for (left) weak fluctuations and (right) strong fluctuations as a function of density ratio.

FIG. 12.

Dependence of the eigenvalues on the fluctuations strength for the conservative system with energy fluctuations (a) at some density ratio; (b) top—for a broad range of fluctuations’ strength at some density ratio; bottom—for (left) weak fluctuations and (right) strong fluctuations as a function of density ratio.

Close modal

2. Stability of fundamental solutions

In the limiting cases of weak (f → 0) and strong (f → ∞) fluctuations for fluids with very different densities (R → ∞) and with very similar densities (R → 1+) in Eqs. (18), the properties of the dynamics are the following (Table VI and Fig. 12).

TABLE VI.

Attributes of the fundamental solution for the conservative system with energy fluctuations with the asterisk marking a quantity.

f ≡ 0f > 00 < f < fcrfcr < f
r1 
Re[ω] <0 <0 >0 
C 
r2 
Re[ω] <0 <0 <0 
C 
r3 
Re[ω] >0 >0 >0 >0 
C 
r4 
Re[ω] <0 <0 <0 <0 
C 
f ≡ 0f > 00 < f < fcrfcr < f
r1 
Re[ω] <0 <0 >0 
C 
r2 
Re[ω] <0 <0 <0 
C 
r3 
Re[ω] >0 >0 >0 >0 
C 
r4 
Re[ω] <0 <0 <0 <0 
C 
a. Weak fluctuations.

For weak fluctuations and very different densities,

(18c)

In Eqs. (18c), solutions ri with i = 1, 2 are stable and describe the decaying oscillations of the flow fields. Solution r3 is the same as in the foregoing cases. Solution r4 is stable (Fig. 12).

For weak fluctuations and similar densities,

(18d)

In Eqs. (18d), the stability of fundamental solutions is similar to the corresponding solutions in Eqs. (18c), except for the difference in the density ratio (Fig. 12).

b. Strong fluctuations.

For strong fluctuations and very different densities,

(18e)

In Eqs. (18e), the fundamental solutions are similar to the corresponding solutions in the classic Landau’s system. Solution r1 is unstable and corresponds to the Landau–Darrieus instability. Solution r2 is formally stable. Solution r3 is the same as before. Solution r4 is stable (Fig. 12).

For strong fluctuations and similar densities,

(18f)

In Eqs. (18f), the stability of fundamental solutions is similar to the corresponding solutions in Eqs. (18e), as expected for the difference in the density ratio (Fig. 12).

3. Flow fields of fundamental solutions

a. Weak fluctuations.

In Eqs. (18), in the case of weak fluctuations and very different densities, solutions rii, ei) with i = 1, 2 describe slightly decaying oscillations of the flow fields (Fig. 12). At the same time, the condition Re [ω1(2)] < 0 leads to an increase in the vortical velocity component far from the interface; thus, the values of integration constants C12,C1=C2* should be zero C1(2) = 0 in order to satisfy the boundary conditions at the outside boundaries at all time. Solution r33, e3) corresponds to the unperturbed flow fields in the entire domain at any time for any C3 as before. For solution r44, e4), due to Re[ω4] < 0, the integration constant must be zero, C4 = 0, in order to satisfy the boundary conditions at the outside boundaries. In Eqs. (18), in the case of weak fluctuations and very similar densities, for solutions rii, ei) with i = 1, 2, 4, the corresponding integration constants Ci should be zero, Ci = 0, in order to satisfy all the boundary conditions. As before, solution r33, e3) corresponds to the unperturbed flow fields at any time for any C3.

b. Strong fluctuations.

In Eqs. (18), in the case of strong fluctuations and very different densities, fundamental solutions rii, ei) with i = 1, 2, 3 are similar to the corresponding solutions in the classic Landau’s system. Specifically, (1) for solution r11, e1), the interface perturbations are strongly coupled with the vortical and potential components of the velocity fields. This solution is unstable. (2) For solution r22, e2), the interface perturbation and vortical and potential velocity components are also strongly coupled. While this solution is formally stable, its velocity field has a vortical structure that decays exponentially in time yet grows away from the interface because of Re[ω2] < 0. For solution r22, e2) to satisfy the boundary condition ul|z→+∞ = 0 at all times, the integration constant must be zero C2 = 0. (3) Solution r33, e3) corresponds to the unperturbed flow fields in the entire domain at any time and for any C3, as before. (4) For the formally stable solution r44, e4) due to Re[ω4] < 0, the integration constant must be C4 = 0 in order to satisfy the boundary condition at the outside boundaries of the domain at any time. In Eqs. (18), in the case of strong fluctuations and similar densities, the solutions are similar to those with very different densities. These similar features include the development of LD-type of instability for solution r11, e1), the choice of C2 = 0 for solution r22, e2) with Re[ω2] < 0, the unperturbed flow fields for solution r33, e3), and the choice of C4 = 0 for solution r44, e4) with Re[ω4] < 0.

4. Interface velocity for fundamental solutions

The properties of fundamental solutions for the conservative dynamics with the effect of energy fluctuations indicate that for f > 0, the interface velocity in the laboratory frame of reference is constant, Ṽ=Ṽ0, for these solutions in Eqs. (18) in the case of either strong f → ∞ or weak f → 0+ fluctuations. Indeed, solution r3 has the null perturbed fields of the velocity and the interface at any integration constant C3. For stable solutions r2(4), the integration constants should be zero to satisfy the boundary conditions at the outside boundaries of the domain, C2(4) = 0. For solution r1 the integration constant should be zero in the case of weak fluctuations, f → 0+, f ≠ 0, in order to satisfy the boundary conditions at the outside boundaries of the domain; in the case of strong fluctuations f → ∞, this solution approaches the classic LD solution r1rLD and has a constant interface velocity Ṽ=Ṽ0, as discussed before. This result holds even when the boundary conditions away from the interface are slightly perturbed, since, for these solutions, the vortical component of the velocity grows exponentially far from the interface.

In the limit of zero fluctuations, f ≡ 0, the interface velocity experiences stable oscillations near the steady value Ṽ0. These oscillations enable the inertial stabilization mechanism of the conservative dynamics to occur, as found in the foregoing.

5. Formal properties of the conservative system with effect of energy fluctuations

To study the formal properties of the dynamic Landau’s system, we find for matrix M=M̃ the associated matrices S=SM̃ and P=PM̃,

(18g)

Equations detPM̃1SM̃ωI=0 and detM̃=0 yield the same eigenvalues ω=ωi,i=1,...,4. The conservative system with the effect of energy fluctuations is non-degenerate and has four fundamental solutions for four degrees of freedom.

6. Summary of properties of fundamental solutions for the conservative system with effect of energy fluctuations

a. Outline of properties of fundamental solutions.

For strong and weak fluctuations and for fluids with similar and contrasting densities, there are four fundamental solutions rii, ei) with i = 1, 2, 3, 4 [Eqs. (18), Table VI, and Fig. 12]. Solution r33, e3) corresponds to the unperturbed flow fields for any C3 and any f. In the case of strong fluctuations, f → ∞, solution r11, e1) describes (for C1 ≠ 0) unstable dynamics of LD-type, whereas for solutions r22, e2) and r44, e4), the constants are C2(4) = 0. In the case of weak fluctuations, f → 0, for solutions r11, e1), r22, e2), and r44, e4), the constants are C1 = C2 = C4 = 0. For f > 0, the interface velocity in the laboratory frame of reference is constant, Ṽ=Ṽ0, for these solutions.

At f ≡ 0, solutions r11, e1) and r22, e2) with non-zero integration constants C1, C2 describe stable oscillations of the flow fields for the conservative dynamics, accompanied by stable oscillations of the interface velocity near the steady value Ṽ0.

b. Transition from the stable to the unstable dynamics.

Transition from the stable to unstable dynamics is illustrated in Fig. 12 showing eigenvalues ωi(R, f) with i = 1, 2, 3, 4 in Eqs. (17) and (18) for a fixed R as a function of f as well as in cases of small and large f. For very small f, the values Im[ω1(2)] and ω3(4) are indistinguishable from the corresponding values in conservative dynamics. For very large f, eigenvalues ωi with i = 1, 2, 3 are indistinguishable from the corresponding values in the classic Landau’s system.

The properties of the flow fields of fundamental solutions for the conservative system with energy fluctuations indicate that when the value f changes from 0 to ∞ in the system [Eqs. (17) and (18)], an abrupt transition can occur to the unstable dynamics. The critical value f = fcr, at which the transition is expected, can be estimated from condition ω = 0 in Eq. (18), leading to fcr = (R + 1)/(R − 1), q = qcr = R(R − 1), and QhQlcr=RR1Vh2.

The flow dynamics can behave as follows: At f ≡ 0, the flow fields experience stable oscillations near the uniform. For 0 < f < fcr, the flow fields are uniform. For f > fcr, the dynamics becomes unstable. For large fluctuations, f → ∞, the unstable dynamics fully recreates the properties of the Landau–Darrieus instability.

In realistic fluids, other mechanisms are possible that account for the influences of entropy flux and the thermodynamic properties of the fluids on energy fluctuations produced by the perturbed interface (Kadanoff, 2000; Landau and Lifshitz, 1987) as well as for nonlinear effects (Sivashinsky, 1983).

In this section, we briefly discuss the connection of our work to the studies in the theory, simulations, and experiments of combustion. We propose new experiments and experimental diagnostics on the basis of our analysis.

While the theoretical framework (Landau, 1994) describes the dynamics of a phase boundary broadly defined, LDI is traditionally associated with the theory of combustion. Landau (1944) considered the evolution of a flame in premixed combustion by treating the unburned gas as a heavy fluid, the burned gas as a light fluid, both fluids as ideal and incompressible, and the interface as a discontinuity.

A premixed flame is a self-sustained wave of an exothermic reaction in a gaseous mixture (Sivashinsky, 1983; Williams, 1965; and Zeldovich, 1944). The wave is initiated by a spark, propagates through the mixture, and produces combustion products. A freely propagating flame is a multi-scale, multi-parameter, multi-physics, and multi-chemistry system. Its evolution is a long-standing problem (Peters, 2000). For planar and steady flames, the interface dynamics is well captured by rigorous theories (Mayo et al., 1990; Frank-Kamenetsky, 1969; Zeldovich and Frank-Kamenetsky, 1938; and Mallard and Le Chatelier, 1883). For unsteady and curved flames, the flame interaction with the hydrodynamic flow that it produces is a challenging problem (Class et al., 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Liberman et al., 1994; Sivashinsky, 1983; Matalon and Matkowsky, 1982; and Williams, 1965).

Landau’s result on the unconditional instability of the flame front contradicts some experiments because stable flames exist and are observed in a laboratory (Clavin, 2000; Clanet and Searby, 1998; Buckmaster, 1979; and Markstein, 1951). Several approaches were developed to explain these observations. Theoretical research was focused on bridging the gap between the idealistic framework and realistic conditions and on augmenting the theory (Landau, 1944) with physical effects typical for flames in premixed combustion (Clavin and Searby, 2016; Class et al., 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Liberman et al., 1994; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Williams, 1971; Sivashinsky, 1983; Matkowsky and Sivashinsky, 1979; Istratov and Librovich, 1966; and Zeldovich, 1944). These effects included the finite thickness of the flame and the finite curvature of the interface (Markstein, 1951), the molecular diffusion and the thermal expansion (Williams, 1971), the dynamics of a reaction zone inside a thin flame (Matalon and Matkowsky, 1982), and the kinematic viscosity (Frankel and Sivashinsky, 1982) as well as the assumptions of weak compressibility and a gradual change in the flow quantities across the interface (Clavin, 2000; Williams, 1971). These fundamental works have found that various physical effects can stabilize the development of LDI at small scales (Sivashinsky, 1983). Peters (2000) and Clavin (2000) summarized these studies, reconciled the theory and the experiment, and found that in premixed combustion, the flame stability depends on the interplay of compressibility, mass diffusion, and heat diffusion. Class et al. (2003a; 2003b), Klimenko and Class (2000), and Sivashinsky (1983) analyzed the stability of a uniformly propagating planar flame as a solution of the unified model that accurately captured Landau–Darrieus, cellular, and pulsating instabilities including a high-wavenumber cutoff for each.

Our results agree with Landau’s results (Landau, 1944) and thus substantiate that the theories of premixed combustion are accurate and complete, within the range of their applicability (Class et al., 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Sivashinsky, 1983; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Sivashinsky, 1980; Matkowsky and Sivashinsky, 1979; Williams, 1971; Markstein, 1951). Particularly, by applying proper transformations of variables, one can directly match our results and the classical results (Landau and Lifshitz, 1987). This includes the governing equations, the boundary conditions, the structure of the solution, the characteristic length scale of the vortical field, and the instability growth-rate. Furthermore, the flow fields in Figs. 8 and 9 can be obtained with some effort from classical works (Peters, 2000; Landau and Lifshitz, 1987; Sivashinsky, 1983; Landau, 1944; Zeldovich, 1944).

In experiments, the implementation and diagnostics of Landau–Darrieus instability is an extreme challenge. It is remarkable that the first experimental study of LDI (Clanet and Searby, 1998) was reported nearly 60 years after the theoretical work was published (Landau, 1944). These experiments have measured the LDI growth-rate in the compressible reactive fluids (gases) for a relative thick interface (Markstein parameter 4.0) in the presence of acceleration (gravity) under the influence of acoustic parametric forcing. More recently, the observations of Landau–Darrieus instability have also been reported for a Bunsen burner and for an inverted V-flame (Clavin, 2000). These experiments have observed that under the influence of Landau–Darrieus instability, the amplitude of the flame wrinkles increases; the instability growth-rate agrees with the theory accounting for the finite flame thickness and non-ideal effects (Markstein parameters 4.3, 4.6, 5.0) (Clavin, 2000). A direct quantitative comparison of these experimental results and our theory is a challenge since the experimental conditions strongly deviate from our theoretical approximation (Clavin, 2000; Clanet and Searby, 1998).

For numerical modeling of Landau–Darrieus instability of combustion fronts in premixed gases, the simulations apply the so-called low Mach number approximation and consider non-ideal fluids with simplified chemical kinetics and gradually changing flow fields (Schlimpert et al., 2016; Miglani et al., 2014; Clavin, 2000; and Denet and Haldenwang, 1995). Particularly, simulations (Denet and Haldenwang, 1995) adopt the isobaric approximation valid for low Mach number flows, apply momentum–pressure formulation for variable density flows, and solve numerically the equations for a premixed flame in a continuous non-ideal fluid. The simulations obtain curved flames with large amplitudes and provide measurements of the growth-rate of LDI (Denet and Haldenwang, 1995). This dispersion relation is consistent with theory (Clavin, 2000; Pelce and Clavin, 1982) and detects the stabilization of short wavelengths.

In addition to the quantitative results in a broad parameter regime, the simulations (Denet and Haldenwang, 1995) also present the flow fields, including temperature, velocity streamlines, vorticity lines, and the interfacial shear. The presence of the vortical field in the bulk of the light fluid is consistent with traditional theories of combustion and with our results for LDI (Fig. 8) (Clavin, 2000; Williams, 1971; Landau, 1944). The simulations also report the presence of shear at the interface with a finite thickness (Denet and Haldenwang, 1995), whereas for Landau’s solution, the velocity field is shear free at the discontinuous interface [Eqs. (2) and (14) and Fig. 8] (Landau and Lifshitz, 1987). The absence of interfacial shear at the discontinuous interface with interfacial mass flux is a rigorous theoretical result for ideal fluids [Eqs. (2)] (Abarzhi et al., 2019; Ilyin et al., 2018; and Landau and Lifshitz, 1987). Further investigations are required to better understand the processes of formation of shear at the interface and the vortical field in the bulk in real fluids.

Numerical simulations of premixed combustion, while employing realistic conditions of weakly compressible non-ideal fluids, usually presume the constancy of the interface velocity (Schlimpert et al., 2016; Clavin, 2000; Denet and Haldenwang, 1995). The direct observation of the inertial stabilization mechanism, which leads to the stable oscillations of the interface velocity near the constant value, and which our analysis finds for ideal incompressible fluids, may be a challenge to detect in these simulations. Theories of premixed combustion also often presume the constancy of the interface velocity (Class et al., 2003a; 2003b; Clavin, 2000; and Peters, 2000). For instance, this assumption is applied in the studies of the influence of temperature dependence of diffusivities on the dynamics of flame fronts. This analysis has derived the 4 × 4 matrix, has found the eigenvalues for the matrix, and has identified the instability growth-rate (Clavin, 2000). A derivation of eigenvectors for the 4 × 4 matrix may further enable the study of structure of flow fields, as in Eqs. (1)–(18).

The existing experiments and simulations of the interfacial dynamics are focused on the evolution of the perturbation amplitude (Class et al., 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Clanet and Searby, 1998; Sivashinsky, 1983; Matalon and Matkowsky, 1982; and Matkowsky and Sivashinsky, 1979). Our theoretical analysis, in addition to identifying the evolution of the perturbation amplitude, finds that the flow is highly sensitive to the interfacial boundary conditions and that the properties of interfacial transports at microscopic scales can be linked to the flow fields at macroscopic scales. Note that Figs. 3–6 and 8–10 present the perturbations of the flow fields. In order to obtain the flow fields, the corresponding zero-order component should be added to the perturbations.

Figure 13 illustrates the velocity streamlines Sh(l) of the heavy and light fluids for the conservative dynamics and the classic Landau’s dynamics at some density ratio and some instance of time. Each plot has its own range of values. The streamlines are defined as dShl/dt×vhl, where the heavy (light) fluid velocity is vhl=Vhl+uhl. The integration constants for solutions rCD and rLD are chosen such that the amplitude of the interface perturbation is the same for the conservative dynamics and the classic Landau’s dynamics and equals to z¯=1/3. Figure 13 illustrates slight deviations of the velocity streamlines from the straight lines for both conservative dynamics and classic Landau’s dynamics. For the conservative dynamics, the deviations of the velocity streamlines are clearly seen in the heavy fluid; they are also present in the light fluid though they are somewhat less pronounced due to the large value of the uniform velocity of the light fluid. For the classic Landau’s dynamics, the deviations of the velocity streamlines are clearly seen in both fluids. For a given value of the perturbation amplitude, the fluid motion is more intense in the classic Landau’s dynamics (due to the presence of the vortical components of the velocity of the light fluid) when compared to that in the conservative dynamics. The results in Figs. 4, 8, and 13 are in consistency with one another.

FIG. 13.

Velocity streamlines for the conservative (left) and the classic Landau’s (right) dynamics.

FIG. 13.

Velocity streamlines for the conservative (left) and the classic Landau’s (right) dynamics.

Close modal

Note that due to the small values of the perturbed velocities, uhlVhl, the properties of the flow fields in the fluid bulk may be a challenge to precisely diagnose. For addressing these challenges, advanced experimental technologies may be applied (Orlov et al., 2010). New experiments with new diagnostics can be conducted on the basis of our results to better understand the interfacial dynamics. Here, we outline some parameters of a new experiment.

According to our results, in the experiment, the interface velocity is sub-sonic, the dynamics is observed from a far field, and its characteristic length scale 1/k = λ/2π and time scale τ = 1/kVh = λ/2πVh are greater than those that are induced by the processes of dissipation and diffusion and by a finite thickness of the interface. These scales are set by the initial conditions. To select the density ratio for the experimental fluid system, we note that, in the classical Landau’s solution for a given λ, the length-scale λ̃ of the vortical field is a non-monotone function of the density ratio R with λ̃/λ=k/k̃=R1+R/R+R+R2+R3. The minimum value of λ̃/λ is achieved at R=2+54.24. Thus, the vortical field can be easier to diagnose when the fluids density ratio is R=2+54.24 and the ratio of uniform velocities is Vl/Vh=2+54.24.

In this case, the Landau solution has the following characteristics:

(19)

or, numerically, R = 4.24, λ̃/λ=4.24, ω = 1.00, and e1 = (−0.62 i, 1.6 i, 0.62 i, 1). The choice of the integration constant for the solution is suggested by the requirement that the perturbations of the flow fields are small. By evaluating the amplitude of the interface perturbation for which linear theory is valid as z¯=1/3, we estimate the integration constant as C=5.4×101. By implementing in experiments the fluid density R=2+54.24 and the corresponding initial conditions, one can further diagnose the flow fields in the vicinity of the interface and in the bulk and compare with the theoretical solutions.

The experimental outcome may include the following scenarios: (1) The interface is stable; the perturbed flow fields are the same as those of the conservative solution; they slightly oscillate near the equilibrium. This can be interpreted as the absence of energy fluctuations and/or as their weak effect. (2) The interface is unstable; the flow fields are the same as in the classical Landau’s solution. This can be interpreted as the existence of LDI and the applicability of Landau’s framework for the experimental system. (3) The interface is unstable; the flow fields deviate qualitatively from those of the classical LDI and from those of the conservative dynamics with strong energy fluctuations. This can be interpreted as that the approximation of ideal incompressible fluids is not a valid approximation for the fluid system. It may also suggest that the fluid system experiences a hydrodynamic instability that is different from LDI. If the growth-rate of this instability agrees with the results of the seminal theories of combustion, it may give grounds to name those instabilities in honor of those who studied its properties (Class et al., 2003a; 2003b; Klimenko and Class, 2000; Peters, 2000; Clavin, 2000; Clanet and Searby, 1998; Liberman et al., 1994, Sivashinsky, 1983; Matalon and Matkowsky, 1982; Pelce and Clavin, 1982; Matkowsky and Sivashinsky, 1979; Williams, 1971; and Markstein, 1951).

To even more accurately outline the regime of experimental parameters, our analysis can be expanded to include compressibility and non-ideal effects. We address these studies to the future.

In reactive fluids, it is generally well understood that chemically reactive systems can be hydro-dynamically unstable (Ilyin et al., 2019; Zhakhovsky et al., 2014; Mayo et al., 1990). Yet, it is a challenge to construct a model experiment or a simulation that clearly displays a significant chemical reaction instability at a simple interface. We can use atomistic simulations to study the energy transport in a reactive system and the effect of chemical reaction on the interface stability. Here, we briefly outline the properties of a simple chemically reactive system. By applying reactive molecular dynamics simulations, we find the following: (i) Simple chemical reactions can be accompanied by complex processes with transient states (Ilyin et al., 2019). (ii) Energy fluctuations due to chemical reactions can destabilize the interface, leading to chemistry-induced instability. (iii) Energy scatter (via rotational and vibrational degrees of freedom rather than translational degrees of freedom) may hinder the instability development. Our results are illustrated in Figs. 14–17.

FIG. 14.

Reactive molecular dynamics of hydrogen peroxide decomposition identifies five dominant species and seven dominant reactions.

FIG. 14.

Reactive molecular dynamics of hydrogen peroxide decomposition identifies five dominant species and seven dominant reactions.

Close modal
FIG. 15.

Reactive molecular dynamics of the formation of the interface with the radicals introduced initially at the interface. (a) Simulation system for computational investigation of instability showing dimensions, bulk HOOH, location of the radical-initiated flat front, and location of the inset. (b) Propagation and disappearance of a perturbation in the inset. Energetic radicals are propelled forward and initiate the formation of a plume at 10 ps; however, the plume fails to propagate at 12 ps due to energy scatter into the bulk; by 24 ps, the plume has effectively disappeared into the rest of the flame front.

FIG. 15.

Reactive molecular dynamics of the formation of the interface with the radicals introduced initially at the interface. (a) Simulation system for computational investigation of instability showing dimensions, bulk HOOH, location of the radical-initiated flat front, and location of the inset. (b) Propagation and disappearance of a perturbation in the inset. Energetic radicals are propelled forward and initiate the formation of a plume at 10 ps; however, the plume fails to propagate at 12 ps due to energy scatter into the bulk; by 24 ps, the plume has effectively disappeared into the rest of the flame front.

Close modal
FIG. 16.

Reactive molecular dynamics simulations of the dynamics of the interface with the radicals introduced initially at the interface. A detailed look at the propagation of a radical-initiated front at the same location as Fig. 15(b). A perturbation is formed around 1 ps and grows over the next 16 ps. However, it does not propagate forward, and by 24 ps, the rest of the front “catches up” and the perturbation disappears.

FIG. 16.

Reactive molecular dynamics simulations of the dynamics of the interface with the radicals introduced initially at the interface. A detailed look at the propagation of a radical-initiated front at the same location as Fig. 15(b). A perturbation is formed around 1 ps and grows over the next 16 ps. However, it does not propagate forward, and by 24 ps, the rest of the front “catches up” and the perturbation disappears.

Close modal
FIG. 17.

A detailed look at the propagation of a thermal shock at the same location as in Fig. 15(b). Instead of radical initiation, the system is thermally shocked in the region outlined in white by rescaling atom velocities. The resulting temperature map shows the location of high-energy molecules. The temperature front also propagates in a stable manner.

FIG. 17.

A detailed look at the propagation of a thermal shock at the same location as in Fig. 15(b). Instead of radical initiation, the system is thermally shocked in the region outlined in white by rescaling atom velocities. The resulting temperature map shows the location of high-energy molecules. The temperature front also propagates in a stable manner.

Close modal

In our reactive molecular dynamics simulations, we find that reactions that are simple at first glance may exhibit complex transient behavior (Ilyin et al., 2019). As a particular example, we consider radical initiated hydrogen peroxide decomposition to water and oxygen. We start with the reactants hydrogen peroxide (HOOH) and hydroxyl radical (OH) and allow them to react under various temperatures and pressures to form a series of intermediate species (OH and HOO in steady-state as well as occasional OH–OH and HOO–HOOH complexes) to products (HOH and O2). Here, we use standard notations O for oxygen and H for hydrogen.

We analyze this system by applying the RMD2Kin method (Ilyin et al., 2019). The molecular dynamics is run with the reactive force field ReaxFF, which allows electric charges to flow freely and permits continuous bond breaking and formation (van Duin et al., 2001). Molecules are detected as they form and break bonds, and the chemical reactions are analyzed directly with no preconceived ideas about the chemistry. We follow the products back in time to find the specific reaction that formed them. We then statistically analyze all reactions at all time steps and at various temperatures to extract reaction rates and barriers using the Arrhenius relation (Ilyin et al., 2019).

The computational setup is as follows: The system is enclosed in a flat periodic box 20 × 2 × 100 nm3 in size. At the start of the simulation, the box is filled with 100 000 peroxide molecules. After equilibration and heating, reactions are initiated by introducing radicals in a specific region [Fig. 15(a)] or by thermal shock in a specific region (Fig. 17). The system is then allowed to react for tens to hundreds of picoseconds.

From the reaction trajectories, the observed dynamics are quite complex. Many reactions can be described as two bodies colliding to form two products (OH + HOOH → H2O + HOO) or three products (HOO + HOOH → O2 + H2O + OH). For uni-molecular reactions (HOOH → OH + OH), third bodies play the role in energizing the reacting molecules without being modified. However, we also find additional interactions with no net reaction that are essential to accurately describe the reactive process. For instance, two OH radicals may end up in a rotationally bound state that temporarily reduces their concentration and affects the overall reaction rate. As another example, HOO and HOOH can transfer the H back and forth, which effectively removes some HOO radicals from interactions and also has to be accounted for.

We find five dominant species and seven dominant reactions. The five dominant species are HOOH, H2O, O2, HOO, and OH. The seven dominant reactions are (1) HOOH + OH → HOO + H2O, (2) HOO + HOOH → O2 + H2O + OH, (3) HOO + OH → O2 + H2O, (4) HOO + HOO → O2 + HOOH, (5) HOO + H2O → OH + HOOH, (6) HOOH + M (a third body) → OH + OH, and (7) OH + OH <=> HO–OH (see Fig. 14). We find that each reaction transfers the translational motion of the reactants into rotational and vibrational energy of the products. Thus, the kinetic energy of the reactants is effectively absorbed into molecular degrees of freedom.

To study the effect of energy fluctuations on the interface dynamics, we carry out atomistic simulations aimed specifically at searching for chemically induced instabilities such as LDI. We use two ways of initiating the reactions: introduction of radicals (Figs. 15 and 16) and thermal shock (Fig. 17).

When the radicals are introduced, as in Figs. 15 and 16, then, at the interface between a region with HOOH seeded with OH and region of pure HOOH, we observe the initiation of a local plume developing into HOOH. While the plume starts to grow, it, however, is hindered before developing a clear instability (Figs. 15 and 16). When the thermal shock is introduced, as in Fig. 17, then, after some growth, small perturbations along the thermal shock interface are hindered as well (Fig. 17). We observe that at late times, both methods produce a stable front propagation.

In Fig. 15, energetic radicals are propelled forward and initiate the formation of a plume at 10 ps; however, the plume fails to propagate at 12 ps due to energy scatter into the bulk, and by 24 ps, the plume has effectively disappeared into the rest of the flame front. Similarly, in Fig. 16, a perturbation is formed around 1 ps and grows over the next 16 ps. However, it does not propagate forward, and by 24 ps, the rest of the front “catches up” and the perturbation disappears (Fig. 16). A thermal shock in the same region also does not lead to instability. The temperature map (Fig. 17) shows the location of high-energy molecules, and the temperature front propagates in a stable manner.

In all systems studied, small local perturbations appear but then quickly disappear into the rest of the reaction front. As a result, the front is not perfectly flat but propagates in a stable way. Any initial forward momentum of the front is rapidly scattered, and any initial perturbations that may have led to an LDI-like instability completely disappear.

We hypothesize two main reasons for the hindrance of perturbations. In order to develop into an instability, a perturbation needs (1) forward momentum and (2) exponential growth. First, the momentum that could potentially cause formation and propagation of a plume comes from the individual radicals. These radicals are smaller than the bulk HOOH molecules, and thus, any initial forward momentum is rapidly scattered after a few collisions. Second, the reactions in this particular chemistry are such that each reaction consumes and produces one radical. There is no multiplication of the number of radicals or chain reaction. In fact, studies of similar chemistry have shown that the number of radicals remains constant throughout (Ilyin et al., 2019). Therefore, the perturbation does not experience exponential growth, and the evolution of the interface is stabilized.

These issues may be solved by (i) introducing outside forces that may provide the necessary forward momentum, (ii) selecting chemistry that yields more easily a radical chain reaction, (iii) including fast propagation of high energy atoms and radicals, and (iv) selecting the proper initial conditions for the flow fields. For instance, we might see more radical multiplication that could lead to the exponential growth of the plume by changing the reactants or parameters of the system to an extreme, e.g., by replacing HOOH with analogs that produce two or more radicals per each consumed radical, such as HOOOH or HOOOOH, or by using a different chemistry altogether. To further disturb the flat interface, the reactive system may be augmented with the fast propagation of high energy atoms and radicals, which, in turn, may initiate the reaction and the formation of hotspot ahead of the interface.

Emphasize that in our reactive molecular dynamics simulations, the energy is effectively scattered via rotational or vibrational degrees of freedom (rather than diffused). This process is substantially quicker than the energy dissipation (diffusion) via translational degrees of freedom, which require multiple collisions. Our theoretical analysis at continuous scales suggests the importance of the large-scale vortical field in the bulk of the reacted material (in order for the instability to occur). Its implementation in reactive molecular dynamics simulations is challenging, and we address the resolving of this challenge to the future research.

Hence, atomistic simulations of combustible systems find that even a simple chemical reaction can be accompanied by complex processes with transient states strongly influencing energy transport. Energy fluctuations by chemical reactions can potentially lead to chemistry-induced instabilities. The hindering of the instability can occur via energy scatter locally into molecular degrees of freedom (rather than the longer-range energy diffusion due to multiple collisions). Further investigations are required to better understand the properties of chemistry-induced instabilities at atomic and continuous scales. Our analytical and numerical methods can serve for accurate modeling of reactive systems and their hydrodynamic stability and for expanding the theory of combustion beyond the diffusive approximation.

The inertial stabilization mechanism of conservative dynamics suggests that for ideal fluids, the interface velocity experiences small stable oscillations near the constant value. In realistic fluids, these free oscillations can be damped by dissipation, diffusion, and other effects and can lead to the constancy of the interface velocity. A question thus appears whether the inertial stabilization mechanism can be directly observed in experiments and simulations. To answer this question, we analyze the interface dynamics influenced by acceleration. Here, we briefly outline the properties of accelerated conservative dynamics and compare them with those of inertial conservative dynamics as well as with accelerated Landau–Darrieus and Rayleigh–Taylor dynamics. Details can be found in Ilyin et al. (2018).

Under the influence of acceleration g directed from the heavy to the light fluid, g = (0, 0, 1), pressure is modified as PP + ρ g z. Accounting for the structure of the solution [Eqs. (10) and (11)], the pressure perturbations are transformed to phl=ρhlΦ̇hl+VhlΦhl/zgz. By introducing the dimensionless acceleration value G=g/kVh2 with G > 0, one can further find, for accelerated conservative dynamics, the matrix MG and the corresponding fundamental solutions. While matrices MG and M appear similar [except for element 33, which is G(R − 1) in matrix MG and is 0 in matrix M], a number of properties of accelerated conservative dynamics deviate substantially from those of inertial conservative dynamics. Particularly, accelerated conservative dynamics exhibits a new fluid instability (Ilyin et al., 2018).

Per condition detMG = 0, the fundamental solutions for accelerated conservative dynamics are

(20)

where Gcr = R(R − 1)/(R + 1) and asterisks mark functions of R, G. Fundamental solutions r1(2) have potential velocity fields in the fluid bulk for any G ≥ 0 and are transformed to the corresponding solutions in Eqs. (12) at G = 0. Fundamental solutions r3(4) are the same as the corresponding solutions in Eqs. (12) for any acceleration value G (Figs. 5 and 6) (Ilyin et al., 2018).

For small acceleration values G < Gcr, solutions r1, r2 describe stable traveling waves whose interference results in stably oscillating standing waves, similarly to Figs. 3 and 4. For large acceleration values G > Gcr, solutions r1(2) describe a standing wave with an exponentially growing (decaying) amplitude. The stable and unstable regions in the (G, R) plane are separated by the curve GR(R − 1)/(R + 1) = 0 and the curve gkVh2ρh/ρlρhρl/ρh+ρl=0 in dimensional form. This suggests that the stability of accelerated conservative dynamics is defined by a balance of inertial and acceleration effects, that is, by reactive force and gravity. For G < Gcr, the inertial effect ∼kVh2 dominates and solution rCDG = (r1 + r2)/2 is stable. In this regime, the interface is stable and the interface velocity experiences stable oscillations near the constant value, similarly to the case of inertial conservative dynamics [Eqs. (12)]. For G > Gcr, acceleration effect ∼g dominates and solution rCDG = r1 is unstable. In this regime, the interface is unstable and the interface velocity increases with time (Ilyin et al., 2018).

The properties of this new fluid instability deviate substantially from those of the accelerated Landau–Darrieus instability and Rayleigh–Taylor instability (Ilyin et al., 2018). Particularly, the new fluid instability develops for acceleration values greater than a threshold, G > Gcr, whereas the accelerated Landau–Darrieus and Rayleigh–Taylor instabilities develop for any acceleration value G > 0. For accelerated conservative dynamics, the velocity field is potential in the bulk, similarly to Rayleigh–Taylor unstable dynamics and unlike accelerated Landau–Darrieus dynamics having the coupled potential and vortical components of the velocity field. The accelerated conservative dynamics is shear-free at the interface, similarly accelerated Landau–Darrieus dynamics and unlike Rayleigh–Taylor unstable dynamics exhibiting the interfacial shear. The interface velocity is time-dependent for accelerated conservative dynamics, is constant for accelerated Landau–Darrieus dynamics, and is zero for Rayleigh–Taylor unstable dynamics. At a given density ratio R > 1 and at large enough acceleration values G > G*, where G* = (R2 − 1)/4 with G* > Gcr, the new fluid instability has the largest growth-rate when compared to the accelerated Landau–Darrieus instability and Rayleigh–Taylor instability (Ilyin et al., 2018).

We emphasize that the velocity fields have the same qualitative features for the accelerated and inertial fluid systems. For conservative dynamics, the velocity field is potential in the bulk and is shear-free at the interface. For Landau’s dynamics, the velocity field has the coupled potential and vortical components in the bulk and is shear-free at the interface. Pressure fields are, however, distinct in the inertial and accelerated cases (Figs. 4 and 8) (Ilyin et al., 2018). In the presence of acceleration, the perturbed pressure fields become asymmetric in the light and heavy fluids. This anisotropy leads to the formation of the structure of bubbles and spikes at the interface and is present in the accelerated conservative dynamics and accelerated Landau–Darrieus and Rayleigh–Taylor dynamics (Ilyin et al., 2018).

Our theoretical framework is developed to study interface dynamics in a multiphase flow. The interface is broadly defined as a phase boundary. It can be an interface between two distinct fluids, an interface between distinct phases of the same fluid, or an interface at which fluids experience a phase transition and change their thermodynamic properties (Abarzhi et al., 2019). In addition to combustible systems, some examples of multiphase flows with interfacial mass flux may include, e.g., explosion of type-Ia supernova, convection in stellar interiors, indirect-drive inertial confinement fusion, and transportation of liquefied natural gas (Abarzhi et al., 2019; Karthik et al., 2018; Arjomandnia et al., 2014; Hurricane et al., 2014; Nepomnyashchy et al., 2012; Bo et al., 2011; Gorokhovski and Herrmann, 2008; Prosperetti, 2017; Bell et al., 2004a; 2003b; Stein and Norlund, 2000; Clavin, 2000; Peters, 2000; and Sivashinsky, 1977).

Our theoretical framework can be employed to study multiphase flows with sharply changing flow fields in nature and technology (Abarzhi et al., 2019; Karthik et al., 2018; Arjomandnia et al., 2014; Hurricane et al., 2014; Nepomnyashchy et al., 2012; Bo et al., 2011; Prosperetti, 2017, Bell et al., 2004a; 2003b; and Peters, 2000). Our analysis finds strong sensitivity of the interface dynamics to the interfacial boundary conditions, directly links microscopic processes at the interface to macroscopic properties of the flow fields in the bulk, and elaborates extensive theoretical benchmarks for experiments and simulations Eqs. (1)–(18) and Figs. 1–12. To the best of the authors’ knowledge, these benchmarks were not discussed before. They can serve for advancement of methods of experimental diagnostics and numerical modeling of multi-phase flows in a broad parameter regime and can expand the capabilities of accurate methods of numerical modeling of the interface dynamics, including the level set and the in-phase fluid method (Karthik et al., 2018; Bo et al., 2011; Bell et al., 2004a; 2003b; Badalassi et al., 2003; Osher and Fedkiw, 2002; Peters, 2000; and Anderson et al., 1998).

Our theoretical approach accurately describes the formal properties of the interface dynamics, finds the degeneracy of the classic Landau’s solution, and connects the lifting of this degeneracy to a singular, possibly—self-similar, character of the dynamics. Our theory finds that for the conservative dynamics, the interface between two fluids, which may mix at small scales due to the interfacial mass flux (i.e., are miscible at molecular scales), can be stable at global scale. One possible application of our results in multiphase geophysical flows can be a qualitative explanation of confluence of rivers—when waters in two rivers meet and not mix at the global scale (Fig. 18).

FIG. 18.

Possible application of our analysis in multiphase geophysical flows. Two rivers meet and do not mix at the confluence of the Green and the Colorado Rivers in Canyonlands National Park, UT, USA (https://www.pinterest.com.au/pin/321092648426240027/). The interface between miscible (red and green) fluids is stable at the global scale.

FIG. 18.

Possible application of our analysis in multiphase geophysical flows. Two rivers meet and do not mix at the confluence of the Green and the Colorado Rivers in Canyonlands National Park, UT, USA (https://www.pinterest.com.au/pin/321092648426240027/). The interface between miscible (red and green) fluids is stable at the global scale.

Close modal

In order to study realistic multiphase flows, including premixed combustion, accounting for compressibility and non-ideal effects, our analysis can be augmented with the conditions for entropy conservation, the corresponding boundary conditions, and the equation of state. The formalism developed in this work provides the principal opportunity to systematically study the influence of realistic properties of the fluids on the quantitative and qualitative properties of interface dynamics and the structure of flow fields. We address these studies to the future.

We have studied the dynamics of the interface separating ideal fluids with different densities and having interfacial mass flux [Eqs. (1)–(3)]. We have developed and applied the general matrix method to solve the boundary value problem for the system of governing equations linearized near the equilibrium state of the uniform flow fields [Eqs. (1)–(11)]. We have formally derived the fundamental solutions, have found their eigenvalues and eigenvectors, and have directly linked the interface stability to the structure of the flow fields (Figs. 1–13 and Tables I–VI). Our theoretical framework is consistent with and generalizes the classic approach (Landau and Lifshitz, 1987) and assumes sharp changes in the flow quantities at the interface and the negligible effects of dissipation, diffusion, compressibility, surface tension, and the interface thickness. We have identified the extreme sensitivity of the interfacial dynamics to the interfacial boundary conditions [Eqs. (12)–(18)]. We have found the new properties of the interface dynamics for the conservative and the Landau systems and the new mechanisms of the interface stabilization and destabilization [Eqs. (12)–(18) and Figs. 1–13].

Particularly, the conservative dynamics is stable, and it is stabilized by the reactive force and the inertial effects. The stable solutions are traveling waves; their interference results in the stably oscillating standing waves. The stable flow is a superposition of two motions—the background motion of the fluids following the interface, whose velocity slightly oscillates near the steady value, and the stable oscillations of the interface perturbations [Eqs. (12) and Figs. 1, 2, 4, and 11]. For the classic Landau’s dynamics, the interface velocity is constant, the reactive force is absent, and the dynamics is unstable. The unstable solution is the standing wave with the increasing amplitude. The LD unstable flow is a superposition of two motions—the background motion of the fluids following the interface with the constant velocity and the growth of the interface perturbations [Eqs. (14) and Figs. 1, 7, 8, and 11]. For the stable conservative dynamics, the velocity fields are potential in the fluids’ bulk (Fig. 4). For the classic unstable LD dynamics, the potential and vortical components of the velocity fields are strongly coupled with the interface perturbation (Fig. 8). In either case, there is the null velocity shear at the interface. In either case, the velocity streamlines slightly deviate from straight lines, though for the unstable Landau’s dynamics, the deviations are stronger when compared to those in the conservative dynamics (Fig. 13).

Our results clearly indicate that the problem of the Landau–Darrieus instability remains a source of inspiration for theory research. Landau (1944) analyzed the interface stability by applying compatibility conditions and identifying two eigenvalues of the dynamics (Landau, 1944). Over the past 70 years, the theoretical framework of Landau (1944) has been extended to incorporate the properties of realistic environments (Class et al., 2003a; 2003b; Peters, 2000; Clavin, 2000; Sivashinsky, 1983; Matalon and Matkowsky, 1982; Williams, 1971; and Markstein, 1951). Our work focuses on general conditions of the stability of the interface with the interfacial mass flux, considers incompressible ideal fluids from a far field, and treats the interface as a discontinuity [Eqs. (1)–(16)]. We have developed and employed the general matrix method and have found the formal fundamental solutions for the linearized dynamics, which directly link the interface stability to the structure of the flow fields’ (Figs. 1–13). To our knowledge, this rigorous consideration has not been performed before.

By applying our general matrix method, we have found that the conservative system is non-degenerate and has four fundamental solutions associated with four independent degrees of freedom [Eqs. (12), Tables I and III, and Figs. 2–6]. These include the two stable traveling waves, whose interference results in two stably oscillating standing waves, and the two other solutions—the unstable solution with the null perturbation fields of the velocity and pressure and the stable solution with the null integration constant. We have further found that the classic Landau’s system is degenerate and has only three fundamental solutions [Eqs. (14), Tables II–IV, and Figs. 5 and 7–9]. These include the unstable standing wave corresponding to LDI, the stable solution with the null integration constant, and the unstable solution with the null perturbation fields of the velocity and pressure, which is the same as for the conservative dynamics. The degeneracy is associated with the static nature of the special postulated condition at the interface (Landau, 1944) and can be lifted with the use of a more general condition.

When the degeneracy is lifted by modifying the postulated condition of zero perturbed mass flux by a constant mass flux across the interface [Eqs. (15)], the fourth solution appears [Eqs. (16), Tables II and V, and Fig. 10]. This solution is neutrally stable. It has a null eigenvalue, and its vortical field can seed the LDI. The null eigenvalue implies that the perturbations of the flow fields can be scale-invariant power-law functions of time rather than scale-dependent exponential functions. Furthermore, by noting that the interface is a phase boundary at which one fluid is transformed into another, this result is in agreement with the linear response theory of phase transitions, where thermodynamic fluctuations are known to be power-laws rather than exponentials (Landau and Lifshitz, 1987; Kadanoff et al., 1967; Kadanoff, 2000).

We have identified the new stabilization mechanism of the conservative dynamics. It is due to the inertial effects that cause the reactive force to occur and, by influencing the interface velocity, to balance the change of momentum produced by the interface perturbation. Note that this inertial stabilization mechanism leads to the new hydrodynamic instability of the accelerated interface that develops when the acceleration is directed from the heavy fluid to the light fluid and when the acceleration magnitude exceeds a threshold (Abarzhi et al., 2019; Ilyin et al., 2018). Both in the inertial and accelerated conservative dynamics, the flow has potential velocity fields in the fluid bulk and is shear free at the interface (Abarzhi et al., 2019; Ilyin et al., 2018).

In LDI, the inertial stabilization mechanism is absent due to the perfect constancy of the interface velocity. Furthermore, in ideal incompressible fluids, the development of the LDI is associated with an imbalance of perturbed energy due to the energy fluctuations by the perturbed front. To model the effect of energy of fluctuations (i.e., the imbalance of the perturbed energy), we have introduced an artificial energy flux and have found that the LDI may develop when the fluctuations are strong. The structure of the flow fields of the solutions for the conservative dynamics with the effect of energy fluctuations suggests a possibility of an abrupt transition to unstable dynamics with the increase in fluctuation strength over a critical value (Fig. 12). In realistic systems, the energy imbalance can be induced by chemical reactions and by difference in the thermodynamic properties of the fluids (Abarzhi et al., 2019; Ilyin et al., 2019).

To our knowledge, these properties of the conservative dynamics and the Landau dynamics have not been discussed in detail in other studies (Abarzhi et al., 2019). Traditionally, it is believed that an interface separating nearly ideal incompressible fluids and having an interfacial mass flux is a subject to the LDI at large scales and that in a realistic environment, the LDI is a challenge to implement because the effects of dissipation, diffusion, and finite interface thickness stabilize the small scales (Peters, 2000). Our analysis is in agreement with these results: For fluids with similar densities, the contribution of vortical fields to the dynamics is small, and the fluxes induced in realistic fluids by the stabilizing effects may dominate the dynamics and define the interface stability (Abarzhi et al., 2015). For fluids with very different densities, a more careful consideration is required. Our analysis yields qualitative and quantitative characteristics of the dynamics that have not been measured before and that require new diagnostics (Figs. 1–13) (Abarzhi et al., 2019; 2015).

One such experiment can be a study of the dynamics of fluids with very different densities, with diagnostics of the flow fields near the interface and in the fluids’ bulk, and with the measurements of the interface evolution, including the interface velocity as a whole and the interface perturbation growth-rate. By comparing the observations with our benchmarks, one can further identify the fundamentals of the interfacial dynamics in realistic environments and elaborate new approaches for the flow control (Figs. 1–13). Several questions may frame these studies: Can the LDI unconditionally develop at the large scales? How can the dynamics be stabilized–by inertial effect, by dissipation and diffusion, or by their combination? How strong should energy fluctuations be to destabilize the flow? Can these fluctuations be induced by chemical reactions? How are the properties of chemistry-induced instabilities compared to those of the LDI?

Our studies of chemically reactive flows by means of reactive molecular dynamics find the following: Even a simple chemical reaction can be accompanied by complex transient processes with strong energy fluctuations. Energy fluctuations by chemical reactions can lead to chemistry-induced instabilities. The hindering of the instability can occur via the energy scatter. Further investigations are required to better understand the properties of chemistry-induced instabilities at atomic and continuous scales.

The existing experimental and numerical studies of the interface stability are focused on the growth of the perturbation amplitude (Abarzhi et al., 2019). Our analysis derives the amplitude growth-rate and finds that the dynamics is highly sensitive to interfacial boundary conditions. According to our theory, by measuring, at macroscopic scales, the flow fields in the bulk and at the interface, one can capture the properties of the interfacial transport at microscopic scales (Figs. 1–13). This information is particularly important in systems where experimental parameters can be difficult to control and where data are a challenge to obtain. They include, for instance, supernovae in astrophysics, laser ablation in plasma physics, plumes in the solar atmosphere, as well as convection in planetary interiors, and industrial applications (Abarzhi et al., 2013). By directly connecting the macroscopic observations far from the interface to the microscopic properties at the interface, one can further get a better understanding of connection between the Lagrangian and Eulerian dynamics away from equilibrium and potentially achieve a better control of complex transport processes involving phase transitions (Abarzhi et al., 2019).

Our approach can be employed in a broad range of problems with interfacial dynamics. Particularly, our theoretical framework can be expanded to other multiphase flows, such as ablative Rayleigh–Taylor instabilities, detonation to deflagration transition, convection in stellar and planetary interiors, and stability of the liquid–liquid and liquid–vapor interface (Abarzhi and Goddard, 2019; Abarzhi et al., 2019; Ilyin et al., 2019; 2018; Remington et al., 2019; Zhakhovsky et al., 2019; 2014; Schlimpert et al., 2016; Hurricane et al., 2014; Anisimov et al., 2013; Abarzhi, 2010; Gorokhovski and Herrmann, 2008; Meshkov, 2006; Bell et al., 2004a; 2003b; Class et al., 2003a; 2003b; Stein and Norlund, 2000; Peters, 2000; Clavin, 2000; Klimenko and Class, 2000; Piriz et al., 1997; Sanz, 1994; Prosperetti and Plesset, 1984; Sivashinsky, 1983; and Williams, 1965). When augmented with linear response theory, it can also serve to analyze dynamic processes in fluids experiencing phase transitions and super-critical fluids (Kadanoff, 2000; Landau and Lifshitz, 1987; and Kadanoff et al., 1967). Our approach can self-consistently account for the effects of compressibility, viscosity, heat conduction, mass ablation, and mass diffusion and thus enable a unified treatment of the interfacial dynamics in realistic environments (Landau and Lifshitz, 1987). We address these studies to the future.

To conclude, we have analyzed from a far field the dynamics of the interface separating incompressible ideal fluids with different densities and with mass flux across it. We have found fundamental solutions of the dynamics and directly linked the interface stability to the flow fields’ structure, associating rigorous mathematical attributes to physical observables. Our approach is consistent with the spirit of Landau’s approach and identifies the new properties of the interface evolution and the new mechanisms of the interface stabilization and destabilization. Our molecular dynamics simulations illustrate the complexity of chemically reactive systems and suggest possible chemistry-induced instabilities.

The authors contributions to this work are as follows: D.V.I.—analysis, methodology, and investigations for the general theoretical framework and the reactive molecular dynamics simulations, writing—original draft, and figures’ preparation; W.A.G.—analysis, methodology, and investigations for the reactive molecular dynamics simulations and chemistry-induced instabilities, student’s supervision, and writing—original draft; S.I.A.—conceptualization, analysis, methodology, and investigations for the general theoretical framework, writing—original draft, and writing—review and editing.

The methods, the results, and the data presented in this work are freely available to the readers and on request from the authors.

This work is dedicated to the memory of Leo P. Kadanoff.

The authors thank for support the University of Western Australia (Australia), Award/Contract Number 10101047 [S.I.A.]; the National Science Foundation (USA), Award/Contract Number 1404449 [S.I.A.]; the Japan Society for the Promotion of Science (Japan), Award/Contract Number S11140 [S.I.A.]; the Summer Undergraduate Research Fellowship Program at California Institute of Technology (USA), the Toni and Bob Perpall Fellowship [D.V.I.]; the Department of Energy (USA), Award/Contract Number DE-SC0017710 [W.A.G.]. S.I.A. expresses her gratitude to the late Professor Leo P. Kadanoff and the late Professor Sergei I. Anisimov for discussions on singularities and stabilization mechanisms and thanks Professor Yasuhide Fukumoto for discussions on fluid instabilities. The authors express their gratitude to Professor Sivashinsky for comments and for discussions of self-similar dynamics.

1.
Abarzhi
,
S. I.
, “
Review of theoretical modelling approaches of Rayleigh-Taylor instabilities and turbulent mixing
,”
Philos. Trans. R. Soc. A
368
,
1809
(
2010
).
2.
Abarzhi
,
S. I.
,
Fukumoto
,
Y.
, and
Kadanoff
,
L. P.
, “
Stability of a hydrodynamic discontinuity
,”
Phys. Scr.
90
,
018002
(
2015
).
3.
Abarzhi
,
S. I.
,
Gauthier
,
S.
, and
Sreenivasan
,
K. R.
, “
Turbulent mixing and beyond: Non-equilibrium processes from atomistic to astrophysical scales II
,”
Philos. Trans. R. Soc. A
371
,
20130268
(
2013
);
Turbulent mixing and beyond: Non-equilibrium processes from atomistic to astrophysical scales I
,”
371
, 20120436 (
2013
);
Turbulent mixing and beyond: Non-equilibrium processes from atomistic to astrophysical scales
,”
371
, 20120435 (
2013
).
4.
Abarzhi
,
S. I.
and
Goddard
,
W. A.
, “
Interfaces and mixing: Nonequilibrium transport across the scales
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
18171
(
2019
).
5.
Abarzhi
,
S. I.
,
Ilyin
,
D. V.
,
Goddard
,
W. A.
, and
Anisimov
,
S. I.
, “
Interface dynamics: Mechanisms of stabilization and destabilization and structure of flow fields
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
18218
(
2019
).
6.
Aglitskiy
,
Y.
 et al., “
Basic hydrodynamics of Richtmyer–Meshkov-type growth and oscillations in the inertial confinement fusion-relevant conditions
,”
Philos. Trans. R. Soc. A
368
,
1739
(
2010
).
7.
Anderson
,
D. M.
,
McFadden
,
G. B.
, and
Wheeler
,
A. A.
, “
Diffuse-interface methods in fluid mechanics
,”
Annu. Rev. Fluid Mech.
30
,
139
165
(
1998
).
8.
Anisimov
,
S. I.
,
Drake
,
R. P.
,
Gauthier
,
S.
,
Meshkov
,
E. E.
, and
Abarzhi
,
S. I.
, “
What is certain and what is not so certain in our knowledge of Rayleigh–Taylor mixing?
,”
Philos. Trans. R. Soc. A
371
,
20130266
(
2013
).
9.
Arjomandnia
,
P.
,
Tade
,
M. O.
,
Pareek
,
V.
, and
May
,
E. F.
, “
Analysis of available data from liquefied natural gas rollover incidents to determine critical stability ratios
,”
AIChE J.
60
,
362
374
(
2014
).
10.
Arnett
,
D.
,
Supernovae and Nucleosynthesis
(
Princeton University Press
,
Princeton
,
1996
).
11.
Arnett
,
D.
,
Meakin
,
C.
, and
Young
,
P. A.
, “
Turbulent convection in stellar interiors. II. The velocity field
,”
Astrophys. J.
690
,
1715
(
2009
).
12.
Azechi
,
H.
 et al., “
Comprehensive diagnosis of growth rates of the Ablative Rayleigh–Taylor instability
,”
Phys. Rev. Lett.
98
,
045002
(
2007
).
13.
Badalassi
,
V. E.
,
Ceniceros
,
H. D.
, and
Banerjee
,
S.
, “
Computation of multiphase systems with phase field models
,”
J. Comput. Phys.
190
,
371
397
(
2003
).
14.
Barenblatt
,
G. I.
,
Zeldovich
,
Y. B.
, and
Istratov
,
A. G.
, “
On diffusion thermal stability of laminar fame
,”
Z. Prikl. Mekh. Tekh. Fiz.
4
,
21
(
1962
).
15.
Beklemishev
,
D. V.
,
Course of Analytic Geometry and Linear Algebra
(
Nauka
,
Moscow
,
1971
) (in Russian).
16.
Bell
,
J. B.
,
Day
,
M. S.
,
Rendleman
,
C. A.
,
Woosley
,
S. E.
, and
Zingale
,
M.
, “
Direct numerical simulations of type Ia supernovae flames. I. The Landau–Darrieus instability
,”
Astrophys. J.
606
,
1029
(
2004a
).
17.
Bell
,
J. B.
,
Day
,
M. S.
,
Rendleman
,
C. A.
,
Woosley
,
S. E.
, and
Zingale
,
M.
, “
Direct numerical simulations of type Ia supernovae flames. II. The Rayleigh–Taylor instability
,”
Astrophys. J.
608
,
883
(
2004b
).
18.
Bo
,
W.
,
Liu
,
X.
,
Glimm
,
J.
, and
Li
,
X.
, “
A robust front Tracking method: Verification and application to simulation of the primary breakup of a liquid jet
,”
SIAM J. Sci. Comput.
33
,
1505
1524
(
2011
).
19.
Buckmaster
,
J. D.
, “
The quenching of two-dimensional premixed flames
,”
Acta Astronaut.
6
,
741
(
1979
).
20.
Chertkov
,
M.
and
Yakhot
,
V.
, “
Propagation of a huygens front through turbulent medium
,”
Phys. Rev. Lett.
80
,
2837
(
1998
).
21.
Clanet
,
C.
and
Searby
,
G.
, “
First experimental study of the Darrieus–Landau instability
,”
Phys. Rev. Lett.
80
,
3867
(
1998
).
22.
Class
,
A. G.
,
Matkowsky
,
B. J.
, and
Klimenko
,
A. Y.
, “
Stability of planar flames as gasdynamic discontinuities
,”
J. Fluid Mech.
491
,
51
63
(
2003a
).
23.
Class
,
A. G.
,
Matkowsky
,
B. J.
, and
Klimenko
,
A. Y.
, “
A unified model of flames as gasdynamic discontinuities
,”
J. Fluid Mech.
491
,
11
49
(
2003b
).
24.
Clavin
,
P.
, “
Dynamics of combustion fronts in premixed gases: From flames to detonations
,”
Proc. Combust. Inst.
28
,
569
585
(
2000
).
25.
Clavin
,
P.
and
Searby
,
G.
,
Combustion Waves and Fronts in Flows
(
Cambridge University Press
,
2016
).
26.
Darrieus
,
G.
, “
Propagation d’un front de flame
,” presented at La Technique Moderne, Paris, France (unpublished) (
1938
).
27.
Davies
,
R. M.
and
Taylor
,
G. I.
, “
The mechanics of large bubbles rising through extended liquids and through liquids in tubes
,”
Proc. R. Soc. London, Ser. A
200
,
375
(
1950
).
28.
Dell
,
Z. R.
,
Pandian
,
A.
,
Bhowmick
,
A. K.
,
Swisher
,
N. C.
,
Stanic
,
M.
etal., “
Maximum initial growth-rate of strong-shock-driven Richtmyer-Meshkov instability
,”
Phys. Plasmas
24
,
090702
(
2017
).
29.
Denet
,
B.
and
Haldenwang
,
P.
, “
A numerical study of premixed flames Darrieus–Landau instability
,”
Combust. Sci. Technol.
104
,
143
167
(
1995
).
30.
Frankel
,
M. L.
and
Sivashinsky
,
G.
, “
The effect of viscosity on hydrodynamic stability of a plane flame front
,”
Combust. Sci. Technol.
29
,
207
224
(
1982
).
31.
Frank-Kamenetsky
,
D. A.
,
Diffusion and Heat Transfer in Chemical Kinetics
(
New York, Plenum
,
1969
).
32.
Gorokhovski
,
M.
and
Herrmann
,
M.
, “
Modeling primary atomization
,”
Annu. Rev. Fluid Mech.
40
,
343
366
(
2008
).
33.
Hillebrandt
,
W.
and
Niemeyer
,
J. C.
, “
Type Ia supernova explosion models
,”
Annu. Rev. Astron. Astrophys.
38
,
191
(
2000
).
34.
Hurricane
,
O. A.
 et al., “
The high-foot implosion campaign on the National Ignition Facility
,”
Phys. Plasmas
21
,
056314
(
2014
).
35.
Ilyin
,
D. V.
,
Fukumoto
,
Y.
,
Goddard
,
W. A.
, and
Abarzhi
,
S. I.
, “
Analysis of dynamics, stability, and flow fields’ structure of an accelerated hydrodynamic discontinuity with interfacial mass flux by a general matrix method
,”
Phys. Plasmas
25
,
112518
(
2018
).
36.
Ilyin
,
D. V.
,
Goddard
,
W. A.
 et al., “
First-principles-based reaction kinetics from reactive molecular dynamics simulations: Application to hydrogen peroxide decomposition
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
18202
(
2019
).
37.
Istratov
,
A. G.
and
Librovich
,
V. B.
, “
The effect of transport processes on the stability of a plane flame front
,”
Prikl. Mat. Mekh.
30
,
541
(
1966
).
38.
Kadanoff
,
L. P.
,
Statistical Physics: Statistics, Dynamics and Renormalization
(
World Scientific
,
2000
).
39.
Kadanoff
,
L. P.
etal., “
Static phenomena near critical points: Theory and experiment
,”
Rev. Mod. Phys.
39
,
395
(
1967
).
40.
Kadau
,
K.
,
Barber
,
J. L.
,
Germann
,
T. C.
,
Holian
,
B. L.
, and
Alder
,
B. J.
, “
Atomistic methods in fluid simulation
,”
Philos. Trans. R. Soc.
368
,
1547
(
2010
).
41.
Karthik
,
K.
,
Dominic
,
K.
, and
Herrmann
,
M.
, “
An in-cell reconstruction finite volume method for flows of compressible immiscible fluids
,”
J. Comput. Phys.
373
,
784
810
(
2018
).
42.
Klimenko
,
A. Y.
and
Class
,
A. G.
, “
On premixed flames as gasdynamic discontinuities: A simple approach to derive their propagation speed
,”
Combust. Sci. Tech.
160
,
23
37
(
2000
).
43.
Kull
,
H. J.
and
Anisimov
,
S. I.
, “
Ablative stabilization in the incompressible Rayleigh-Taylor instability
,”
Phys. Fluids
29
,
2067
(
1986
).
44.
Landau
,
L. D.
, “
On the theory of slow combustion
,”
Acta Phys. U.R.S.S.
19
,
77
(
1944
).
45.
Landau
,
L. D.
and
Lifshitz
,
E. M.
,
Course of Theoretical Physics VI, Fluid Mechanics
(
Pergamon Press
,
New York
,
1987
).
46.
Liberman
,
M. A.
etal., “
Stability of a planar flame front in the slow-combustion regime
,”
Phys. Rev. E
49
,
445
(
1994
).
47.
Lord Rayleigh
, “
Investigations of the character of the equilibrium of an incompressible heavy fluid of variable density
,”
Proc. London Math. Soc.
14
,
170
(
1883
).
48.
Mallard
,
E.
and
Le Chatelier
,
H. L.
, “
Combustion des mélanges gazeux explosifs
,”
Ann. Mines. Ser.
8
(
3
),
274
(
1883
).
49.
Markstein
,
G. H.
, “
Interaction of flow pulsations and flame propagation
,”
J. Aeronaut. Sci.
18
,
428
(
1951
).
50.
Marmottant
,
P.
and
Villermaux
,
E.
, “
On spray formation
,”
J. Fluid Mech.
498
,
73
(
2004
).
51.
Matalon
,
M.
and
Matkowsky
,
B. J.
, “
Flames as gasdynamic discontinuities
,”
J. Fluid Mech.
124
,
239
(
1982
).
52.
Matkowsky
,
B. J.
and
Sivashinsky
,
G. I.
, “
An asymptotic derivation of two models in flame theory associated with the constant density approximation
,”
SIAM J. Appl. Math.
37
,
686
(
1979
).
53.
Mayo
,
S. L.
,
Olafson
,
B. D.
, and
Goddard
,
W. A.
, “
Dreiding: A generic force field for molecular simulations
,”
J. Phys. Chem.
94
,
8897
8909
(
1990
).
54.
Meshkov
,
E. E.
,
Studies of Hydrodynamic Instabilities in Laboratory Experiments
(
FGYC-VNIIEF
,
Sarov
,
2006
) (in Russian), ISBN: 5-9515-0069-9.
55.
Meshkov
,
E. E.
and
Abarzhi
,
S. I.
, “
Group theory and jelly’s experiment of Rayleigh–Taylor instability and Rayleigh–Taylor interfacial mixing
,”
Fluid Dyn. Res.
51
,
065502
(
2019
).
56.
Miglani
,
A.
,
Basu
,
S.
, and
Kumar
,
R.
, “
Insight into instabilities in burning droplets
,”
Phys. Fluids
26
,
032101
(
2014
).
57.
Nepomnyashchy
,
A.
,
Simanovskii
,
I.
, and
Legros
,
J. C.
,
Interfacial Convection in Multilayer Systems
, 2nd ed. (
Springer
,
New York
,
2012
).
58.
Orlov
,
S. S.
,
Abarzhi
,
S. I.
,
Oh
,
S. B.
,
Barbastathis
,
G.
, and
Sreenivasan
,
K. R.
, “
High performance holographic technologies for fluid-dynamics experiments
,”
Philos. Trans. R. Soc., A
368
,
1705
1737
(
2010
).
59.
Osher
,
S.
and
Fedkiw
,
R.
,
Level Set Methods and Dynamic Implicit Surfaces
(
Springer
,
2002
).
60.
Pelce
,
P.
and
Clavin
,
P.
, “
Influence of hydrodynamics and diffusion upon the stability limits of laminar premixed flames
,”
J. Fluid Mech.
124
,
219
(
1982
).
61.
Peters
,
N.
,
Turbulent Combustion
(
Cambridge University Press
,
2000
).
63.
Piriz
,
A. R.
,
Sanz
,
J.
, and
Ibañez
,
L. F.
, “
Rayleigh–Taylor instability of steady ablation fronts: The discontinuity model revisited
,”
Phys. Plasmas
4
,
1117
(
1997
).
64.
Prosperetti
,
A.
, “
Vapor bubbles
,”
Annu. Rev. Fluid. Mech.
49
,
221
248
(
2017
).
65.
Prosperetti
,
A.
and
Plesset
,
M. S.
, “
The stability of an evaporating liquid surface
,”
Phys. Fluids
27
,
1590
(
1984
).
66.
Rast
,
M. P.
, “
Compressible plume dynamics and stability
,”
J. Fluid Mech.
369
,
125
(
1998
).
67.
Remington
,
B. A.
 et al., “
Rayleigh-Taylor Instabilities in High-Energy Density Settings on the National Ignition Facility
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
18233
(
2019
).
68.
Sanz
,
J.
, “
Self-consistent analytical model of the Rayleigh–Taylor instability in inertial confinement fusion
,”
Phys. Rev. Lett.
73
,
2700
(
1994
).
69.
Schlimpert
,
S.
,
Feldhusen
,
A.
,
Grimmen
,
J. H.
etal., “
Hydrodynamic instability and shear layer effects in turbulent premixed combustion
,”
Phys. Fluids
28
,
017104
(
2016
).
70.
Sethian
,
J. A.
,
Level Set Methods
(
Cambridge University Press
,
1996
).
71.
Shigemori
,
K.
 et al., “
Measurements of Rayleigh–Taylor growth rate of planar targets irradiated directly by partially Coherent light
,”
Phys. Rev. Lett.
78
,
250
(
1997
).
72.
Sivashinsky
,
G. I.
, “
Nonlinear analysis of hydrodynamic instability in laminar flames. I. Derivation of basic equations
,”
Acta Astronaut.
4
(
11-12
),
1177
1206
(
1977
).
73.
Sivashinsky
,
G. I.
, “
On flame propagation under conditions of stoichiometry
,”
SIAM J. Appl. Math.
39
,
67
82
(
1980
).
74.
Sivashinsky
,
G. I.
, “
Instabilities, pattern formation, and turbulence in flames
,”
Annu. Rev. Fluid Mech.
15
,
179
(
1983
).
75.
Stein
,
R. F.
and
Nordlund
,
Å. Å.
, “
Realistic solar convection simulations
,”
Solar Phys.
192
,
91
(
2000
).
76.
van Duin
,
A. C. T.
,
Dasgupta
,
S.
,
Lorant
,
F.
, and
Goddard
,
W. A.
, “
ReaxFF: A reactive force field for hydrocarbons
,”
J. Phys. Chem. A
105
,
9396
9409
(
2001
).
77.
Williams
,
F. A.
,
Combustion Theory
, 2nd ed. (
Addison-Wesley
,
Reading, MA
,
1965
).
78.
Williams
,
F. A.
, “
Theory of combustion in laminar flows
,”
Ann. Rev. Fluid Mech.
3
,
171
(
1971
).
79.
Wu
,
P.-K.
,
Kirkendall
,
K. A.
,
Fuller
,
R. P.
, and
Nejad
,
A. S.
, “
Breakup processes of liquid jets in subsonic crossflows
,”
J. Propul. Power
13
,
64
(
1997
).
80.
Zeldovich
,
Y. B.
,
The Mathematical Theory of Combustion and Explosions
, Academy of Sciences of the Soviet Union (
Springer
,
1944
), ISBN: 978-1461294399.
81.
Zeldovich
,
Y. B.
and
Frank-Kamenetsky
,
D. A.
, “
A theory of thermal propagation of flame
,”
Acta Physicochim. URSS
9
,
341
350
(
1938
).
82.
Zeldovich
,
Y. B.
and
Raizer
,
Y. P.
,
Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena
, 2nd English edtion (
Dover
,
New York
,
2002
).
83.
Zhakhovsky
,
V. V.
,
Budzevich
,
M. M.
,
Landerville
,
A. C.
,
Oleynik
,
I. I.
, and
White
,
C. T.
, “
Laminar, cellular, transverse, and multiheaded pulsating detonations in condensed phase energetic materials from molecular dynamics simulations
,”
Phys. Rev. E
90
,
033312
(
2014
).
84.
Zhakhovsky
,
V. V.
etal., “
Mass and heat transfer between evaporation and condensation surfaces: Atomistic simulation and solution of Boltzmann kinetic equation
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
18209
(
2019
).