Numerical simulation of compressible flows using two-phase observable Euler equations

TL;DR

This paper introduces an observable regularization technique for simulating two-phase compressible flows, effectively handling high wavenumber modes and material interfaces, with demonstrated accuracy and computational efficiency.

Contribution

It proposes a novel inviscid regularization method based on observable regularization for two-phase compressible flows, enabling efficient and accurate simulations across different flow regimes.

Findings

The method accurately captures material interfaces and shock interactions.

Significant computational savings are achieved, especially in 3D simulations.

Results agree well with analytical, numerical, and experimental data.

Abstract

Most fluid flow problems that are vital in engineering applications involve at least one of the following features: turbulence, shocks, and/or material interfaces. While seemingly different phenomena, these flows all share continuous generation of high wavenumber modes, which we term the '$k_\infty$ irregularity.' In this work, an inviscid regularization technique called 'observable regularization' is proposed for the simulation of two-phase compressible flows. The proposed approach regularizes the equations at the level of the partial differential equation and as a result, any numerical method can be used to solve the system of equations. The regularization is accomplished by introducing an 'observability limit' that represents the length scale below which one cannot properly model or continue to resolve flow structures. An observable volume fraction equation is derived for capturing the material interface, which satisfies the pressure equilibrium at the interface. The efficacy of the observable regularization method is demonstrated using several test cases, including a one-dimensional material interface tracking, one-dimensional shock-tube and shock-bubble problems, and two-dimensional simulations of a shock interacting with a cylindrical bubble. The results show favorable agreement, both qualitatively and quantitatively, with available exact solutions or numerical and experimental data from the literature. The computational saving by using the current method is estimated to be about one order of magnitude in two-dimensional computations and significantly higher in three-dimensional computations. Lastly, the effect of the observability limit and best practices to choose its value are discussed.

Tables

Table 1. Table 1: Requirements for the averaging kernel.

Properties	Mathematical Expression
Normalized	$\int g$=1
Nonnegative	$g\left(𝐱\right)>0,\forall 𝐱$
Decreasing	$\left|{𝐱}_1\right|\ge \left|{𝐱}_2\right|\Rightarrow g\left({𝐱}_1\right)\le g\left({𝐱}_2\right)$
Symmetric	$\left|{𝐱}_1\right|=\left|{𝐱}_2\right|\Rightarrow g\left({𝐱}_1\right)=g\left({𝐱}_2\right)$

Table 2. Table 2: Comparison of wave and interface velocities in a Mach 1.22 shock-R22 bubble interaction from experiment [ 41 ] with different numerical results from literature and observable results. The grid resolution for each computational result is shown in the parenthesis. V s subscript 𝑉 𝑠 V_{s} : incident shock velocity, V R subscript 𝑉 𝑅 V_{R} : refracted shock velocity, V T subscript 𝑉 𝑇 V_{T} : transmitted shock velocity, V u ​ i subscript 𝑉 𝑢 𝑖 V_{ui} : initial velocity of upstream edge of bubble, V 𝑢𝑓 subscript 𝑉 𝑢𝑓 V_{\mathit{uf}} : final velocity of upstream edge of bubble, and V d ​ i subscript 𝑉 𝑑 𝑖 V_{di} : initial velocity of downstream edge of bubble. Observable results are calculated by a linear least square fitting during the following time periods, 0 ​ μ ​ s ⩽ t V s ⩽ 172 ​ μ ​ s 0 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑠 172 𝜇 𝑠 0\,\mu s\leqslant t_{V_{s}}\leqslant 172\,\mu s , 5 ​ μ ​ s ⩽ t V R ⩽ 172 ​ μ ​ s 5 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑅 172 𝜇 𝑠 5\,\mu s\leqslant t_{V_{R}}\leqslant 172\,\mu s , 199 ​ μ ​ s ⩽ t V T ⩽ 290 ​ μ ​ s 199 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑇 290 𝜇 𝑠 199\,\mu s\leqslant t_{V_{T}}\leqslant 290\,\mu s , 0 ​ μ ​ s ⩽ t V u ​ i ⩽ 172 ​ μ ​ s 0 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑢 𝑖 172 𝜇 𝑠 0\,\mu s\leqslant t_{V_{ui}}\leqslant 172\,\mu s , 376 ​ μ ​ s ⩽ t V 𝑢𝑓 ⩽ 506 ​ μ ​ s 376 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑢𝑓 506 𝜇 𝑠 376\,\mu s\leqslant t_{V_{\mathit{uf}}}\leqslant 506\,\mu s , and 199 ​ μ ​ s ⩽ t V d ​ i ⩽ 269 ​ μ ​ s 199 𝜇 𝑠 subscript 𝑡 subscript 𝑉 𝑑 𝑖 269 𝜇 𝑠 199\,\mu s\leqslant t_{V_{di}}\leqslant 269\,\mu s .

	$V_s$	$V_R$	$V_T$	$V_{\mathrm{𝑢𝑖}}$	$V_{\mathrm{𝑢𝑓}}$	$V_{di}$
Experiment [41]	415	240	540	73	90	78
Quirk and Karni [13] ($\Delta x=56\mu m$)	420	254	540	74	90	116
Shyue [9] (tracking, $\Delta x=125\mu m$)	411	243	538	64	87	82
Shyue [9] (capturing, $\Delta x=125\mu m$)	411	244	534	65	86	98
Observable (this paper, $\Delta x=222\mu m$)	424	247	530	73	88	78
% Error: Quirk and Karni [13]	1.2	5.8	0.0	1.4	0.0	48.7
% Error: Shyue [9] (tracking)	1.0	1.3	0.4	12.3	3.3	5.1
% Error: Shyue [9] (capturing)	1.0	1.7	1.1	11.0	4.4	25.6
% Error: Observable	2.2	2.9	1.9	0.0	2.2	0.0

Equations

$$\overline{f}=g\ast f$$

$$\underset{\Omega }{\iiint }\mathrm{odiv}\mathbf{F}=\underset{\Omega }{\iiint }\left(\overline{f}\mathrm{div}\mathbf{V}+\overline{\mathbf{V}}\cdot \mathrm{grad}f\right)dV=\underset{S}{\iint }\mathbf{F}\cdot \mathbf{n}dS$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho u=0,$$

$$\frac{\partial \rho u}{\partial t}+\nabla \cdot \rho uu+\nabla p=0,$$

$$\frac{\partial \rho E}{\partial t}+\nabla \cdot \left(\rho Eu+pu\right)=0,$$

$$\rho ={\rho }_1{z}_1+{\rho }_2{z}_2,\rho E={\rho }_1{E}_1{z}_1+{\rho }_2{E}_2{z}_2,$$

$$\text{Rate of generation of f}=Q\text{vspace}-2mm$$

$$\frac{\partial f}{\partial t}+\overline{f}\mathrm{div}\mathbf{V}+\overline{\mathbf{V}}\cdot \mathrm{grad}f=Q.\text{vspace}-2mm$$

$$\frac{\partial \rho }{\partial t}+\overline{\rho }\nabla \cdot \mathbf{u}+\overline{\mathbf{u}}\cdot \nabla \rho =0.$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\overline{\rho \mathbf{u}}\nabla \cdot \mathbf{u}+\overline{\mathbf{u}}\cdot \nabla (\rho \mathbf{u})=-\nabla p.$$

$$\frac{\partial \rho E}{\partial t}+\overline{\rho E}\nabla \cdot \mathbf{u}+\overline{\mathbf{u}}\cdot \nabla (\rho E)=-\left(\overline{p}\nabla \cdot \mathbf{u}+\overline{\mathbf{u}}\cdot \nabla p\right)+\rho \mathbf{u}\cdot \mathbf{\Sigma }+S,$$

$$\rho E-\frac{1}{2}\rho u^2=\Gamma p,\Gamma =\frac{1}{\gamma -1},$$

$$\Gamma =\frac{1}{\gamma -1}=\frac{z_1}{{\gamma }_1-1}+\frac{z_2}{{\gamma }_2-1}.$$

$$\frac{\partial \rho }{\partial t}$$

$$\rho \frac{\partial u}{\partial t}$$

$$\frac{\partial \rho e}{\partial t}$$

$$\frac{\partial \rho }{\partial t}+\overline{u}\frac{\partial \rho }{\partial x}=0,$$

$$\frac{\partial \rho e}{\partial t}+\overline{u}\frac{\partial \rho e}{\partial x}=0.$$

$$\Gamma \left(\frac{\partial p}{\partial t}+\overline{u}\frac{\partial p}{\partial x}\right)+p\left(\frac{\partial \Gamma }{\partial t}+\overline{u}\frac{\partial \Gamma }{\partial x}\right)=0.$$

$$\frac{\partial z_1}{\partial t}+\overline{u}\frac{\partial z_1}{\partial x}=0,$$

$$\frac{\partial z_1}{\partial t}+\overline{\mathbf{u}}\cdot \nabla z_1=0.$$

$${k_{S}}=\left\{{\begin{array}[]{*{20}{c}}{2\pi k/L}&{0\leqslant k<N/2}\\ {2\pi(k-N)/L}&{N/2\leqslant k<N}\end{array}}\right.$$

$$H(a,b,x)=\left\{{\begin{array}[]{*{20}{l}}0&{x\leqslant a}\\ {6{\chi^{5}}-15{\chi^{4}}+10{\chi^{3}}}&{a<x<b}\\ 1&{b\leqslant x}\end{array}}\right.$$

$$RH{S}_{Buffer}=H\left(a,b,\xi \right)RH{S}_{Physical}+\left[1-H\left(a,b,\xi \right)\right]\left[H\left(c,d,\xi \right)RH{S}_L+\left(1-H\left(c,d,\xi \right)\right)RH{S}_R\right],$$

$$g^{\alpha }(x)=C_{\alpha }{e}^{\left(-|x|/\alpha \right)},$$

$$f=\overline{f}-{\alpha }^2{\nabla }^2\overline{f}.$$

$$\begin{array}{c}{{\left(\rho ,u,p,\gamma \right)}^T}_{x⩽0}=(1,0.5,1/1.4,1.4{)}^T,\text{hfill}\\ {{\left(\rho ,u,p,\gamma \right)}^T}_{x>0}=(10,0.5,1/1.4,1.2{)}^T.\text{hfill}\end{array}$$

$$\left({{\rho_{1}}z,{\rho_{2}}(1-z),u,p,z}\right)=\left\{{\begin{array}[]{*{20}{l}}{\left({{\rm{1.3764}},0,{\rm{0.3947}},{\rm{1.5698}},1}\right)}&{0\leq x<0.25,}\\ {\left({1,0,0,1,1}\right)}&{0.25\leq x<0.4,}\\ {\left({0,0.138,0,1,0}\right)}&{0.4\leq x<0.6,}\\ {\left({1,0,0,1,1}\right)}&{0.6\leq x<1.}\end{array}}\right.$$

$$\left({{\rho_{1}}z,{\rho_{2}}(1-z),u,p,z}\right)=\left\{{\begin{array}[]{*{20}{l}}{\left({0.386,0,26.59,100,1}\right)}&{-1\leq x<-0.8,}\\ {\left({0.1,0,-0.5,1,1}\right)}&{-0.8\leq x<-0.2,}\\ {\left({0,1,-0.5,1,0}\right)}&{-0.2\leq x<1.}\end{array}}\right.$$

$$\rho e=\Gamma p+\Pi ,\Gamma =\frac{1}{\gamma -1},\Pi =\frac{\gamma p^{\infty }}{\gamma -1}$$

$$\Gamma \left(\frac{\partial p}{\partial t}+\bar{u}\frac{\partial p}{\partial x}\right)+p\left(\frac{\partial \Gamma }{\partial t}+\bar{u}\frac{\partial \Gamma }{\partial x}\right)+\left(\frac{\partial \Pi }{\partial t}+\bar{u}\frac{\partial \Pi }{\partial x}\right)=0.$$

$$\begin{array}{c}(\rho \text{,}u\text{,}P\text{,}\gamma \text{,}{P}_{\infty }{)}_{\text{L}}^{\text{T}}=(1.241\text{,}0\text{,}2.753\text{,}1.4\text{,}0{)}^{\text{T}}\text{,}\text{hfill}\\ (\rho \text{,}u\text{,}P\text{,}\gamma \text{,}{P}_{\infty }{)}_{\text{R}}^{\text{T}}=(0.991\text{,}0\text{,}3.059\times 10^{-4}\text{,}5.5\text{,}1.505{)}^{\text{T}}\text{.}\text{hfill}\end{array}$$

Full text

1 Introduction

Material interfaces and/or shockwaves are the source of complex flow behavior in many fluid engineering problems. Scientific understanding of these problems is the base for technological advances needed in engineering and daily life. For example, shock-induced collapse is a damaging phenomena leading to propeller erosion in naval engineering [1]; similar mechanisms may lead to tissue damage and internal bleeding in shockwave lithotripsy [2, 3, 4], a non-invasive medical treatment of kidney stones. Understanding such phenomena can help us to develop designs that specifically avoid cavitation or control it to our advantage [5]. In the past couple of decades, with the improvement of computational power, numerical study of such problems has become mainstream. While solving the conservation equations, discontinuity in flow variables creates challenges that need specific attention.

In the past three decades, significant effort has been given to numerical treatment of material interfaces and shocks. Many approaches are developed for tracking/capturing interfaces. In Front Tracking [6], a separate unstructured grid is used to follow the location of the interface, and although it has many advantages when it comes to control over the topology change and sharpness of the interface, it has more complexity compared to other methods and also is challenging for parallel processing. Eulerian interface tracking approaches have lower complexities in maintaining and parallel processing of the interface. However, they still need reconstruction of the interface, which adds to computational cost of method. Examples of such methods are ghost-fluid method (GFM) and volume-of-fluid methods. In GFM, a level set function is introduced with a zero value representing the interface. The equation for each fluid is solved on each side of the interface with the help of “ghost-fluid” points that are reconstructed across the interface at each time step [7, 8]. Volume of fluid methods keep the interface sharp by reconstructing the interface from the advected volume fraction [9]. Front capturing methods are the simplest type for multiphase flow. These methods solve an additional equation for evolving the interface without the need for interface reconstruction. The pressure non-equilibrium effect at the interface was initially a downside of these methods as a result of numerical diffusion. Karni [10] first addressed this issue by solving the governing equations in non-conservative form and solving an advection equation for the mass fraction of one of the phases. Many other methods were later developed to address this issue and each solves an additional advection equation for either a material property [11, 12], mass fraction [13], or volume fraction [14]. In addition to development of methods for tracking/capturing the interface, several Riemann solver type approaches are developed for physical approximation of fluxes at discontinuities, specifically for problems including shock waves; Roe solver [15], HLL [16], HLLC [17], and flux vector splitting [18] are a small sample of such methods. All the efforts in solving fluid equations are based on the assumption that any set of conservation equations (e.g. Euler) are valid when solved with finite resolution. This assumption holds in regions of smooth variations where the flow variables and their derivatives can be represented using the limited resolution with minimal error. However, this may not be the case in regions with sharp changes.

When there are sharp variations in a field quantity, one needs to have a large number of wave modes in the Fourier space to correctly represent the data, and as the thickness of the jump decreases the required number of modes approaches infinity. Considering that discontinuous solutions are part of admissible solutions of Euler equations; this infinite resolution requirement is inherent in the Euler equations before any discretization scheme imposed. Consequently, it should not be a surprise that naive discretization of Euler equations could lead to numerical instability in simulation of shocks, interfaces, and even turbulence. In order to avoid this conflict, Mohseni [19, 20] developed the concept of observable divergence in which the finite resolution limit (observability limit) is introduced into the derivation of conservation equations; for a derivation of the observable divergence theorem see [20]. Zhang and Shen [21] demonstrated that observable divergence can be used to efficiently regularize the relativistic Burgers equation. They showed that the regularized equation preserves the shock speed and conservation properties of the original equation. It is believed that observable Euler equations do not allow discontinuous solutions; as a result of this regularity, these equations are expected to have better numerical behavior while handling shocks and interfaces. Using observable divergence for conservation laws requires a reevaluation of the governing equations for the existing problems in order to systematically resolve the challenges in the simulations of fluid flows, including but not limited to multiphase flows, shocks, and turbulence.

In the observable approach the regularization is done at the level of the partial differential equation. The aim of this paper is to numerically demonstrate that the governing equations derived using this concept regularize shocks and interfaces with no need for special numerical interface treatment. To this end, we use pseudo-spectral method which has zero numerical viscosity. Note that the selection of pseudo-spectral method is just to demonstrate that even for a numerical scheme with no numerical dissipation the observable equations are regularized. One can easily use any other numerical scheme. We have also used finite difference and high order Padé schemes in the past. The regularizing capability of the observable set of equations is demonstrated with some examples. We first develop the observable volume fraction equation for a two-phase flow problem. Then, present the numerical solution of 1D and 2D shock-interface interaction problems using the observable set of equations and compare them to available exact solutions or experimental/numerical results from literature. The effect of observability limit on the results is investigated at the end and a guideline for selecting the observability limit is presented.

2 Motivation and Observability Concept

In physics, a field is a physical quantity that can be assigned a value at any and each mathematical point in space and time; in other words a field is assumed to be a continuum. In this section, a critical review of some of the basic concepts and axiomatic assumptions in classical field theory is presented. We further show the intimate connection between these assumptions and many of the most well known, intractable challenges in classical physics such as turbulence and high dimensionality problems in shocks and multiphase flows. The main difficulties in these problems stem from the requirement of resolving fields to infinitely-small scales while in reality no experiments or simulations are capable of such resolution. A resulting consequence is what we refer to as $k_{\infty}$ (wavenumber infinity) irregularity; see the discussion below. Our proposed solution to these issues is to introduce an observable field quantity. Equipped with this new observability concept, we have systematically revisited some of the generally accepted concepts from the calculus and show that one can derive governing observable dynamical equations which are a more realistic representation of the observable physical world. To explain our ideas, in the following we consider three seemingly different phenomenon from fluid mechanics; namely shocks, turbulence, and sharp interfaces in two-phase flows. We will argue that these different challenges are, in fact, manifestations of the infinite observability assumption in the definition of a field quantity. We will then offer a remedy.

*$k_{\infty}$ Irregularity. * To introduce the problem of assuming infinite observability, consider the simplest and standard example of shock formation in the inviscid Burgers’ equation, $u_{t}+uu_{x}=0$. It is well known that a smooth initial condition for the Burgers’ equation can generate a finite-time discontinuity. Formation of a discontinuity is characterized by a rather fast generation of increasingly higher wavemodes as the solution profile steepens toward a shock; see Figure 1. For the one-dimensional Burgers equation the cascade of energy to smaller scales has a slope of $-2$. The conventional wisdom to remedy nonphysical discontinuities has been the addition of a viscous term, e.g. in the form of a Laplacian of the velocity field, in order to control and reduce the creation of smaller scales. The Laplacian term changes the nature of the partial differential equation (PDE) from a hyperbolic PDE to a regularized parabolic one. Note that the term responsible for this high-wavenumber generation is the $\mathbf{u}\cdot\nabla\mathbf{u}$ term which is the difference between an Eulerian and Lagrangian observer; we have termed this phenomenon the high wavenumber ($k_{\infty}$) irregularity.

The $k_{\infty}$ irregularity is a common feature of almost all challenging problems in fluid dynamics including turbulence, shocks and sharp interface two-phase problems where one requires comparing observations between Eulerian and Lagrangian observers ($\mathbf{u}\cdot\nabla\mathbf{u}$ term). The symptoms of the problem in each case (shocks, turbulence, two-phase flows) are identical; namely a never ending tail of the Fourier spectrum with the signal in the physical domain showing formation of a discontinuity or excessive filamentation. Again, the ad hoc remedy in each case has been the introduction of a parabolic term in the form of a Laplacian or its higher order versions. To this end, the $k_{\infty}$ irregularity is responsible for excessive computational cost of high Reynolds number flows, high Mach number shocks, and high density ratio two-phase flows. Interestingly, this $k_{\infty}$ irregularity has not been recognized among these seemingly different problems and as a result different remedies were considered in different disciplines.

Now, a fundamental question is: at which step in deriving the governing PDEs in classical field theory were nonphysical solutions allowed to creep in as a part of the domain of admissible solutions to a physical problem? We have recently argued that the answer to this question is in the infinite observability assumption that is automatically accepted when a physical quantity is represented as a field. In other words, the values of the quantity at any two points can be distinguished from each other no matter how small the distance between the two are! Notice that the wavenumber $k\rightarrow\infty$ is the same as the length scale $L\rightarrow 0$. Basically, by merely defining a field quantity we have axiomatically accepted this hidden assumption. As a result, we have made the assumption that any length scale, no matter how small, is observable by a field quantity. It is clear then that this is not a physically correct assumption. No field quantity can be experimentally or numerically observed to an infinitely-small length scale. This physical limit in our observation is termed the observability scale, $\alpha$. Anything below this scale is virtually unmeasurable or unobservable from a physical stand point. As a result, what we often observe is an averaged quantity rather than the mathematical limit where the observation volume approaches zero. This has significant ramifications when one calculates the flux of a quantity at the face of a volume surrounding a point of interest. This limit is not reflected in the classical field theory equations and we believe that it is the source of many sustained theoretical and numerical difficulties in solving some of the challenging problems in classical field theory, including fluid mechanics.

Upon asking this question a few years ago, we were led to ask a series of others: What would happen if we start with a field quantity with a finite observability limit? Do we obtain the same governing dynamical equations? Do we still obtain the same integral equations in calculus (e.g., divergence and Stokes’ theorems) if the derivations are performed with an observable field quantity and not, as is currently done, with an infinity observable quantity? The following provides answers to some of these questions; see [19, 20] for details.

Observable Field Quantity. An observability limit on a field quantity, $f$, could be imposed by a convolution kernel, $g$, with a length scale $\alpha$ to obtain

[TABLE]

where $\ast$ is the convolution operator. Here, $f$ is the infinitely observable quantity as defined in all classical field theory and $\overline{f}$ is the observable field quantity. The following theorem addresses how the classical divergence theorem is altered for observable field quantities:

Observable Divergence Theorem [19, 20]. Let $\Omega$ be a region with surface boundary $S$ oriented outward. Let $\mathbf{F}(\mathbf{x})=f\mathbf{V}(\mathbf{x})=(u(\mathbf{x}),v(\mathbf{x}),w(\mathbf{x}))f$ be a vector function that is continuous and has continuous first partial derivatives in a domain containing $\Omega$. Then

[TABLE]

where the operator $\bar{(\cdot)}$ is defined by $\bar{(\cdot)}=g\ast(\cdot)$, $\ast$ is the convolution operator, $g$ is the kernel of the convolution with an observability length scale $\alpha$. $g^{\alpha}$ could be any kernel with properties stated in Table 1.

See [19, 20] for a proof.

Observable method was systematically evaluated numerically and analytically for multidimensional Burgers equations [22, 23]. Norgard and Mohseni [22] proved the existence and uniqueness of the solution to observable Burgers’ equation and its convergence to the weak solution of Burgers’ equation as $\alpha\to 0$ for a subset of kernels with the form $\hat{g}(k)=1/\left({1+\sum\nolimits_{j=1}^{N}{{C_{j}}{k^{2j}}}}\right)$; a hat represents the Fourier transform of a quantity and the $C_{j}$ are the kernel constants. Furthermore, the observable Burgers’ equation is a special case of a more general equation, which is shown to be Hamiltonian by Holm and Staley [23].

Equipped with the observable divergence theorem one can systematically derive the observable equations for any conservation laws in classical field theory. In the following section, the conservation equations for a compressible flow are presented first and then an observable volume fraction equation is derived for capturing the material interface and defining material properties inside each phase.

3 Governing Equations for Two-phase Flow

Infinitely observable Euler equations. Conservation of mass, momentum, and energy for the mixture fluid properties are inferred in different ways [14, 24] from Baer-Nunziato model [25] or Drew derivation of two phase conservation equations [26, 27]:

[TABLE]

where $\rho$, $\mathbb{u}$, $p$, and $E$ are density, velocity vector, pressure and specific total energy of fluid mixture, respectively. In the derivation of Eqs. 3, 4 and 5, an instantaneous mechanical relaxation between phases is assumed, $u=u_{1}=u_{2}$ and $p=p_{1}=p_{2}$ where subscript 1 and 2 refer to phase 1 and 2. The mixture quantities are defined using the following relations:

[TABLE]

here, $z_{1}$ and $z_{2}$ are the volume fractions of phase 1 and 2 and $z_{1}+z_{2}=1$. In the rest of this text $z$ with no subscript refers to $z_{1}$.

Observable Euler equations. Mohseni [19, 20] derived the observable mass, momentum, and energy equations. Here, we repeat the derivation for the sake of completeness. We develop the differential equations that must be satisfied by a fluid with an observability limit of $\alpha$. The equations are expected to resolve all relevant dynamical quantities up to this observable scale $\alpha$. As mentioned before, this scale is dictated by the resolution of a particular numerical simulation or the length scale resolution of a particular experimental equipment. For scales below $\alpha$ we assume that the medium still acts as a continuum and one can take the mathematical limit of $\Delta V\rightarrow 0$, where $\Delta V$ is the control volume of interest. However, our observable scale remains at $\alpha$ even as this mathematical limit is taken to zero.

Before the differential form of the observable conservation laws for mass, momentum, and energy of a continuum are presented, we will formulate the general form of the observable differential conservation law for a quantity per unit volume, $f$. We assume $f$ satisfies a conservation law of the form

[TABLE]

where $Q$ is the sum of all sources of $f$ inside the volume. For instance $Q$ is the total external force exerted on a control volume if one considers the conservation of momentum inside the volume. The rate of generation of $f$ is basically the outflow of $f$ minus the inflow at the boundaries of the volume plus the storage inside the volume. The rate of outflow minus the inflow at the boundaries can be conveniently calculated by the operator ${\rm odiv}\,f{\bf V}$ where ${\bf V}$ is the continuum velocity, see Equation 1. Therefore, the observable conservation of $f$ reduces to the following differential equation

[TABLE]

Now, the observable compressible Euler equations can be easily derived by using the conservation of mass ($\rho$), momentum ($\rho\mathbf{u}$), and energy ($\rho E=\rho(e+\frac{1}{2}u^{2})$) to obtain

[TABLE]

where $\mathbf{\Sigma}$ is the sum of all external forces acting on the control volume or its surface and $S$ is the sum of all energy sources inside the control volume. Note that in the limit of $\alpha$ approaching zero, that is when the observable scale approaches a mathematical zero, one recovers the classical Euler equations. For the problems presented in this work $\mathbf{\Sigma}$ and $S$ is assumed to be zero.

In order to close this system of equations, an equation of state (EOS) needs to be used. Here, we use the ideal gas equation of state which is used by many authors for the gas-gas problems [13, 11]:

[TABLE]

where $\gamma$ is the ratio of specific heat of material at constant pressure to its specific heat at constant volume. The mixture $\Gamma$ is defined using:

[TABLE]

3.1 Observable Interface Capturing

To apply the observable Equations 9a, 9b and 9c to two-phase flows, we need to also obtain an equation for capturing the interface, which is needed to define material properties in each phase. In this section we present a derivation of an observable equation for evolving the volume fraction. For the purpose of capturing the interface while keeping pressure equilibrium at the interface, we consider the case of an interface-only problem, following Shyue [12], in which velocity and pressure are constant in the domain while density and material properties change across the interface. For simplicity, we present the one-dimensional equations. Eqs. 9a, 9b and 9c can be rewritten in the following non-conservative form:

[TABLE]

where $e$ is the specific internal energy with $E=e+\frac{1}{2}{u^{2}}$. For the given interface-only problem, the equations for the evolution of interface can be inferred to be:

[TABLE]

By inserting the equation of state into Eq. 16 and rearranging we have:

[TABLE]

With the pressure equilibrium requirement, the first parenthesis is zero. Since this equation needs to be satisfied for any pressure, the term in the second parentheses should be zero. Writing second parenthesis using Eq. 11, the simplified equation is

[TABLE]

which is the observable equation for capturing a material interface in one-dimension and in general form of the equation is:

[TABLE]

In order to validate and verify the applicability of Eqs. 9a, 9b and 9c and Eq. 18 we perform several numerical experiments. In the next section, the numerical method used for this work is explained. Here, we used ideal gas equation of state to derive the volume fraction equation as other authors have done in their work [13, 11]. Using stiffened equation of state, as reported in [28], which is suitable for liquids in addition to gases, the same equation for volume fraction will be achieved as demonstrated in appendix A.

3.2 Numerical Method

As mentioned before, the observable regularization is done at the level of the PDE and not at the discretization level. As a result, one can utilize any numerical technique in solving the observable equations. In the past we have used many different numerical methods including high order finite difference compact Padé in solving observable Burgers’ equation or 1D Euler equations in our group [29]. However, to show the robustness of the regularization in the observable technique, here we used pseudo-spectral method, which provides a numerical technique with no/minimal numerical dissipation as compared with other methods. Our objective is to show that a spectrally high order technique can be used without any numerical dissipation.

The pseudo-spectral method is a standard method [30, 31, 32], which is used here with no modifications. Basically, all spatial derivatives (e.g. $u_{x}$, $u_{xx}$, $u_{yy}$) are calculated in the Fourier domain while all multiplications (e.g. $u\cdot\nabla u$) are done in the physical domain. An explicit method is used for the time integration and a Helmholtz kernel is used for computing the observed quantities.

Computing derivatives. Using pseudo-spectral for taking derivative and applying an averaging kernel such as Helmholtz kernel is very straightforward [30, 31, 32]. Having a function $f$ on $N$ discrete points, ${f^{\prime}_{n}}=f^{\prime}(nL/N)$, where $n$ is the index of the point and $L$ is the length of the domain. Assuming $N$ is even, the scaled wave number is defined as

[TABLE]

and to take the first derivative of a function $f$, first, Fourier coefficients $F_{k}$ for $0\leqslant k<N$ are calculated by taking FFT of $f$ using $f_{n}$ for $0\leqslant n<N$. Then, $F_{k}$ is multiplied by $ik_{S}$ for $k\neq N/2$ (where $i$ is the imaginary unit) and $F_{N/2}$ is multiplied by zero to obtain $F^{\prime}_{k}$ [32]. Finally, ${f^{\prime}_{n}}$ is calculated by taking the inverse FFT of $F^{\prime}_{k}$.

Buffer zone. Since a pseudo-spectral method is used here and a Fast Fourier Transform inherently assumes that the variable is periodic, a buffer zone is added to the sides of the physical domain to satisfy the following conditions: (I) periodicity of all the variables in the computational domain, (II) smooth transition from the condition on the right to the left of the physical domain, (III) prevent reflection from boundaries. To achieve the second condition while satisfying the third condition we use a smooth $5^{th}$ order polynomial as a weighting function to transition from the governing equation in the physical domain to a linear advection equation with a damping term in the buffer zone. The governing equation in the buffer zone is applied using similar ideas proposed by Freund [33]. For the weighting function, the following smooth step function is used:

[TABLE]

where $\chi=\left({x-a}\right)/\left({b-a}\right)$. For any governed quantity, here represented by $f$, the governing equation in the right and left buffer zones are $\frac{{\partial f}}{{\partial t}}=-{U_{R}}\frac{{\partial f}}{{\partial x}}+\sigma\left({{f_{R}}-f}\right)=RHS_{R}$ and $\frac{{\partial f}}{{\partial t}}=-{U_{L}}\frac{{\partial f}}{{\partial x}}+\sigma\left({{f_{L}}-f}\right)=RHS_{L}$, respectively. The first term in these equations is a linear advection with a prescribed constant velocity $U_{R}$ and $U_{L}$ which is set equal to the initial horizontal velocity at the right and left edge of physical domain, respectively. The second term in the buffer zone equations is a sponge term that brings the quantity $f$ towards the prescribed values $f_{R}$ and $f_{L}$, which are again set based on initial value of $f$ at the right and left edge of physical domain. $\sigma$ is a damping coefficient which is here defined to be $1/\Delta x$. In order to find the right hand side of the equation in the buffer zones, the $x$ location of the right boundary is mapped linearly to $\xi\in[0,0.5)$ and the $x$ location of the left buffer zone mapped to $\xi\in[0.5,1)$. Then the right hand side is defined as

[TABLE]

where $RHS_{Physical}$ is the right hand side of the governing equation in the physical domain. $a$, $b$, $c$, and $d$ are 0.1, 0.5, 0.5, and 0.9, respectively, which specify the start and end of the smooth transitions. The values of these constants will be mentioned for each problem in the results section. Note that the buffer zone setup presented here is not necessarily the best choice, but it satisfies the three major requirements mentioned above for the problems solved in this manuscript.

Computation of observed quantities. It is demonstrated for the regularized Burgers’ equation that using different averaging kernels, including Gaussian and Helmholtz, has negligible effect on the large scale features, features above the observability scale $\alpha$, [34, 22]. While any kernel can be used to implement the method proposed in this work, since applying Helmholtz kernel is equivalent to simple algebraic operations in the Fourier space, without loss of generality (see [22]) and for the sake of simplicity, we use this kernel as the averaging kernel defined by,

[TABLE]

where $C_{\alpha}$ is a normalization constant. Alternatively, the Helmholtz kernel can also be applied by solving the Helmholtz equation:

[TABLE]

Computing $\overline{f}$ using a Helmholtz operator, Equation 23, is easily done by multiplying $F_{k}$ by $1/(1+\alpha^{2}k_{S}^{2})$ and taking an inverse FFT.

Time integration and dealiasing. For the time integration, a $3^{rd}$ order TVD Runge-Kutta method is used [35]. For dealiasing, a low pass exponential filter, $G(k)=\exp({-36{{\left|k/k_{\max}\right|}^{36}}})$, is used [36]. In a similar manner to applying Helmholtz kernel, this low pass filter is computed by multiplying $F_{k}$ by $G(k_{S})$ and computing the inverse FFT of the result. The dealiasing filter is applied to all conservative variables in addition to the volume fraction after each Runge-Kutta step. For all initial conditions double filtering is used as explained in [37].

4 Results and discussion

Here, we present several 1D and 2D test cases to show the performance of the proposed method compared to other available methods, experiments and exact solutions.

4.1 One-dimensional advection of isolated interface

First, an isolated interface is examined for a one-dimensional domain with periodic boundary condition and the initial condition defined by

[TABLE]

Since the domain is periodic and velocity is constant across the doma in, after one time period, at $t=4$, all the flow variables should be exactly equal to initial condition as shown in Fig. 2. Observed density and ratio of heat capacities are shown using blue open circles which are in good agreement with the exact solution of the Euler equation (black solid lines). In addition the results of a WENO5 finite volume method calculated by Johnsen and Colonius [38] (orange open triangles) is included. The magnitude of differences between the final and initial observed velocity and pressure are demonstrated to be zero to the order of machine precision, as expected. This shows that our method preserves the pressure equilibrium at material interfaces. The grid resolution is the same as [38], $\Delta x=0.01$, and non-dimensional observability limit, $\alpha/\Delta x$, is set to $0.75$.

4.2 One-dimensional interaction of a shock wave with a Helium bubble

Next a shock-bubble interaction, proposed by Quirk and Karni [13], is studied. In this problem a Mach 1.22 shock wave traveling in air interacts with a Helium bubble. This interaction leads to different wave structures that show the performance of the method in handling all possible waves in an inviscid compressible flow, rarefaction wave, shock wave, and contact discontnuity. The initial condition for this problem is

[TABLE]

The ratio of specific heat, $\gamma$, is set to 1.67 and 1.4, inside and outside the bubble, respectively. Figure 3 shows the results for this problem at $t=0.35$. The solid line is the exact Euler solution. The red dotted line and blue dash-dot line are simulations using 1024 grid points with $\alpha/\Delta x$ equal to 1.25 and 3.75, respectively.

When the incident shock hits the interface, it reflects a left going rarefaction wave while a shock wave passes into the bubble. This shock hits the other end of the bubble and creates a transmitted shock and a reflected shock. These interactions consequently result in two left going and three right going shock waves outside the bubble and a very weak (almost invisible) left going shock wave inside the bubble until $t=0.35$. As shown in the figure, for this problem the method can regularize the equation with a non-dimensional observability limit as small as 1.25. Obviously, as $\alpha/\Delta x$ increases, the equations become more regularized but at the same time less details will be captured for any quantities of interest.

4.3 One dimensional interaction of a strong shock wave with an air bubble

As stated by Coralic and Colonius [39], the interaction of a shock with a material interface can challenge numerical methods that do not solve the equations in conservative form. The main known artifact is likely error in the resulting wave speeds, which deteriorate as time passes and the error in mass conservation increases. To investigate the performance of our method in this regard, the interaction of a strong shock wave with a material interface is studied here.

This problem is originally proposed by Liu et al. [40], and then a modified form of it, as used here, is solved by Johnsen [4] and Coralic and Colonius [39]. In this problem a Mach 8.96 shockwave moving in Helium interacts with an air bubble. Initially, the unshocked region has a small velocity towards the shock. The same resolution as [39] is used, $\Delta x=0.005$, with a non-dimensional observability limit $\alpha/\Delta x=1$. A buffer zone of length 2 is used on both sides of the physical domain. The initial condition is

[TABLE]

Figure 4 shows good agreement between our results and the exact solution and demonstrate that all the wave velocities are captured correctly. Since we are using pseudo-spectral method to calculate derivatives, any sharp changes need to be taken place over several grid points in order to avoid the Gibbs phenomena. As a result of the strong shock, using a non-dimensional observability limit of less than unity creates spurious oscillations in the result, which would be expected as it means the computation is under-observed. Using the Helmholtz kernel, the observability limit imposes the shock thickness of about $4.6\alpha$ (which here translates to about 5 cells or 6 grid points), as shown by Norgard and Mohseni [37]. While, the WENO5 computation of Coralic and Colonius [39] is demonstrating sharper interface and shockwave results using the same resolution, it also creates small amplitude non-physical oscillations noticeable in the pressure and velocity plots for $x\in(-0.4,0.5)$ while this spurious oscillation is not seen in observable results.

4.4 Two-dimensional interaction of a shock wave with a R22 bubble

In this section our method is validated with available experimental results. This problem is a two-dimensional counterpart of the one that is solved in Section 4.2; a Mach 1.22 shockwave in air interacts with a cylindrical R22 bubble with $\gamma=1.249$ and $\rho_{R22}=3.712\,kg/m^{3}$. This problem was first studied experimentally by Haas and Sturtevant [41] and soon became a benchmark for numerical methods, [42, 43, 44, 45, 46, 47, 48], since it has complex wave interactions and also contains a Richtmyer-Meshkov instability that poses a severe challenge on computational methods. The current simulation is done using a grid resolution of $\Delta x=222\,\mu m$ and a non-dimensional observability limit $\alpha/\Delta x=1$. For plotting numerical results, the visualization technique introduced by Quirk and Karni [13] is used here. Figure 5 shows the schematic of the problem setup. At the left and right side of the domain, a buffer zone is used to make the numerical domain periodic (for the purpose of pseudo-spectral computations). Additionally, the buffer zone damps the outgoing waves and minimizes reflection from boundary. The top and bottom boundaries assume a periodic boundary condition which, with the symmetry of the problem, act like a reflective boundary condition. In order to initialize the curved interface on a rectangular grid and minimize the stairs-like initialization of the interface, which could lead to additional waves and interface instability in compressible flow, here a $tanh$ function is used for the definition of initial volume fraction, $z(x,y)=0.5+0.5*tanh((D(x,y)-r_{bubble})/\delta)$, where $r_{bubble}$ is the bubble radius and $D(x,y)$ is the distance from the center of the bubble equal to $\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$, and $\delta$ is approximately one-third of the interface thickness which is set to 1.5 here.

Figure 6 shows experimental schlieren images of shock-cylindrical bubble interaction from Haas and Sturtevant [41] in addition to numerical schlieren from our results. The incident shock hits the interface at $t=0\,\mu s$; a shock wave reflects from the interface and another one refracts into bubble. Since the speed of sound is lower inside the bubble, the refracted shock moves with a lower speed compared to the incident shock outside the bubble, and as a result the shock converges towards the end of bubble; see snapshots (a), (b), and (c). The shock initially reflected from the interface reflects from top and bottom wall, as shown in snapshot (c). The incident shock diffracted around the bubble arrives at the downstream edge of bubble faster than the refracted shock (snapshot d). Consequently, the pressure outside the bubble becomes larger than the pressure of unshocked gas inside the bubble for a short period of time until the refracted shock reaches the downstream edge of the bubble. During this time period the downstream edge starts deflecting towards the upstream edge; this can be noticed by looking closely at red solid line in Fig. 7 between $t=161\,\mu s$ and $t=194\,\mu s$. This deflection is very small and when the refracted shock, which is focused at the downstream edge, hits the interface it pushes back the interface towards the downstream and creates a jet at that point, see snapshots (e) and (f). It is interesting to note that depending on bubble size, shock strength, and the material property of the bubble and surroundings, a reverse jet may run towards upstream, and even pierce the front edge of bubble in spherical bubble case, instead of a jet towards downstream; for such examples see [46] and [49]. If the initial defracted shock has enough time to make a dent into the downstream of the bubble before the refracted shock reaches downstream edge of the interface, the high pressure, created as the refracted shock transmits out of the bubble, pushes the trapped air in the dented area towards upstream and creates a reverse jet as opposed to what is seen here. After the refracted shock is focused at the downstream edge of the bubble, it transmits in an almost circular shape, trailing the defracted wave, snapshots (e) and (f). Interaction of the shock with the interface creates rollup and instability of the interface which become visible in snashots (g) and (h).

The results in Fig. 6 show good agreement between experimental and simulation snapshots. All the main waves and interfaces are captured correctly. To show the accuracy of the results with respect to the experiment, several wave and interface locations over time are extracted from the simulation and compared in Fig. 7 with what is reported in experiment [41]. In this figure the lines show our results while symbols show the result from the experiment. A small schematic of this shock interaction problem shows the location of different waves and interfaces where data are extracted. As shown in Fig. 7, the observable results are in good agreement with the reported experimental results. A quantitative comparison of different wave and interface speeds is shown in Table 2, including experimental results, numerical results from Quirk and Karni [13] and Shyue [9], in addition to the observable results of the current simulation. The reported observable wave and interface speeds are calculated by linear least square fitting the data shown in Fig. 7 for the time intervals mentioned in the caption of the table. As shown in the table, all the reported speeds from the observable method are in good agreement with experimental results and are well within the 10% error reported for the experimental results. Since the resolution of our simulation is at least about twice as coarse as the other reported simulations, while producing similar or better results (e.g. $V_{di}$) compared to the experiment, we estimate that the current method reduces the computational cost by about an order of magnitude in two-dimensions (2 times less grid points in each direction and twice larger time step size based on CFL criterion). Obviously, the savings will be much higher in three-dimensional computations if the current method is used.

4.5 Two-dimensional interaction of a shock wave with a Helium bubble

Haas & Sturtevant [41] also conducted their experiment for a Mach 1.22 shock wave in the air interacting with a cylindrical Helium bubble contaminated with 28% air by mass. In this case, the material properties inside the bubble are $\gamma=1.648$ and $\rho_{He+28\%Air}=0.214\,kg/m^{3}$. The problem setup is exactly the same as the shock-R22 bubble problem explained in the previous section. Since the speed of sound is now higher inside the bubble compared to its surrounding, the refracted shock diverges as opposed to the converging case for the R22 bubble. Because of a higher speed of sound inside the bubble and a higher variation of flow variables across the interface and shocks compared to the R22 bubble case, this case is numerically more challenging.

Figure 8 shows a comparison of the numerical schlieren snapshots of our simulation with the schlieren images from the experiment [41]. In this problem the shock meets the interface at a zero angle at $t=0\,\mu s$, which results in a regular refraction that initially consists of a reflected rarefaction wave and a transmitted shock wave, see snapshot (a). As the incident shock moves downstream, the angle between the shock and interface increases and the reflected wave first turns into a Mach wave and then a shock wave before the refraction transforms into an irregular one. In the irregular refraction, the refracted shock runs ahead of the incident shock and a side shock wave appears, which creates a Mach stem that can be seen clearly in snapshots (b) and (c). For more details on the regular and irregular shock refraction in a slow-fast configuration, the reader is referred to [50, 51, 52]. The refracted shock wave also gets reflected internally, where the shock touches the interface, see snapshots (b) and (c). The numerical schlieren results shown in snapshots (a), (b), and (c) are processed at instants, which more closely resemble the experimental snapshots. As mentioned in the caption of the figure, the timing of these snapshots are slightly different from the corresponding times reported in the experiment. This mismatch is also reported in the work of other authors [39, 48, 53], which might be associated with the timing error in the experiment or the non-uniform contamination which is not accounted for in the simulations. At $t=56\,\mu s$ the refracted shock retransmits into the air, see snapshot (d). The internally reflected shock wave partly transmits and runs behind the doubly transmitted shock wave and partly reflects upstream (circular shape inside the bubble), see snapshot (e). At $t=102\,\mu s$, snapshot (f), the internally reflected shock wave is transmitted into the air and surrounds the bubble. In snapshot (g), the beginning of interface roll up is visible at the upstream side of the bubble; at this time the reflected wave from the top and bottom wall have already interacted with the bubble and can be seen in the front and back of the bubble. The vorticity generated at the interface by the passage of the shockwave creates a jet on the upstream side of the bubble, which evolves the bubble into a kidney-like shape, see snapshot (h). The comparisons in this figure show a good agreement between the interface shape and the wave locations in the observable results with the reported experiment results.

In the set of two-phase Euler equations, the thickness of the interface is theoretically zero. Unless the interface is tracked exactly, as opposed to captured, it is not possible to simulate a zero thickness with any numerical method, this means that it is the numerical method and computational resolution that dictates the smallest scales of the flow close to interface. As a result, most numerical methods cannot have a grid-independent solution, especially when there are small scale features growing next to the interface. Fig. 9 shows the numerical schlieren for the shock-bubble problem described above at $t=245\,\mu s$. Snapshot (a) is the experimental schlieren while snapshot (d) is the result from our method. Snapshot (b) is the numerical schlieren from Quirk and Karni [13], which is achieved by solving the two-phase Euler equation in non-conservative form using an advection equation for the mass fraction, with a computational resolution $\Delta x=56\,\mu m$. Snapshot (d) is the result computed by Coralic and Colonius [39] using a fifth order finite volume WENO method accompanied with a volume fraction equation for capturing the interface. Snapshot (e) is a slightly different simulation in that the bubble is not contaminated with air (pure Helium bubble). It is computed by Hejazialhosseini et al. [54] using a fifth order finite volume WENO method with an advection equation for a smoothed Heaviside function for capturing the interface, included here because of its high spatial resolution. The observable result, (c), more closely resembles the experiment. The three other numerical results show very different interfacial instability depending on the numerical method and the computational resolution. Snapshot (e) with the highest resolution shows that as the resolution is increased the interface instability starts growing much earlier and consequently results in different outcome at a later time. Obviously, by adding viscous terms, and solving the Navier-Stokes equations, some of these issues can be subdued. However, it will add considerable computational cost that might not be nessessary.

4.6 Effect of observability limit and its proper choice

The observability limit specifies how much of high wave numbers or details of the flow will be captured by the set of observable equations. Intuitively, the observability limit needs to be of the order of grid resolution. For moderate shock-tube problems, the minimum non-dimensional observability limit from our experience can be as low as 0.5 ($\alpha=0.5\Delta x$); note that the length scales below $\Delta x$ are not resolved by the limit of computational resolution. However, as the strength of the shock increases and more complexities like interfaces are added to the problem, the minimum acceptable non-dimensional observability limit increases. For the Helmholtz kernel used here, usually a non-dimensional observability limit of slightly higher than one is recommended. Obviously, as the observability limit increases, the equations become more regularized and less susceptible to instabilities with the caveat that some of the small flow structures and details will not be captured. Figure 10 shows an example of this trend demonstrated by solving the same shock-R22 bubble interaction problem for different observability limits. As shown in part (a), (b), and (c) of the figure, by increasing the observability limit, but keeping the same numerical resolution some of the details of the flow are lost. Part (d) of the figure shows the same simulation on a grid twice as coarse as the ones shown in parts (a), (b) and (c) with an observability limit equal to the one used in (b). It can be seen that the results shown in (b) and (d) are almost identical, which demonstrates that while we can increase the observability limit as much as we want, any non-dimensional observability limit that is not as much as possible close to unity will waste computational resources.

5 Conclusion

In this work we used the concepts of observability and observable divergence [20] to develop a method for simulation of two-phase compressible flows. This method uses an observability limit that represents the length scale below which flow structures cannot be resolved. The introduction of the limit is applied at the level of the governing differential equation and therefore is not a filter on an already-computed solution.

An observable equation for tracking the material interface is developed that preserves pressure equilibrium at the interface. A pseudo-spectral method is used to avoid numerical dissipation. Using no special numerical treatment for performing high order calculation, e.g. WENO/ENO, and having no numerical dissipation in the computations shows that it is the governing equations that regularize the evident discontinuities in the problem. A few 1D test cases are studied and the results are compared with some of the best shock and interface capturing methods in the literature and the exact Euler solution. In all cases the current method produce comparable or better results. Additionally, two 2D shock-bubble interactions are studied and the results are compared with available experimental data. The method captured the major wave and interface speeds correctly and the qualitative comparison of numerical schlieren snapshots show good agreement with the experiment. We also compared different results from the literature for the shock-Helium bubble interaction and discussed that because of the infinitely thin interface assumption in the two-phase Euler equations, it is the numerical method and the computational resolution which dictate the thickness of the interface. As a result, depending on the numerical method, different results may be achieved. The observable approach, by introducing the observability limit, defines the smallest length scales in the flow This method reduces the computational cost by approximately one order of magnitude in 2D computations in comparison to other computations and we expect that the savings would be considerably more in 3D. Finally, we also looked at the effect of the observability limit and demonstrated the best practices for selecting the optimum observability limit parameter. Here, an interface-capturing approach is chosen for recognizing the location of the interface since the simplicity of its implementation attracted attention in fluid dynamics community. However, it is necessary to note that the concept of observability can be implemented using other interface-tracking approaches. Since the observable method is regularizing the equations at the level of the PDE, one can use any numerical method to solve the equations and consequently numerical methods for unstructured grids can be used to deal with complex geometries.

Acknowledgment

We thank Dr. Doug Lipinski for his early work on similar problems which helped development of current work. We would also like to thank Dr. Adam DeVoria for his helpful comments in preparation of this manuscript. This work was partially supported in part by the Air Force Office of Scientific Research as well as the National Science Foundation.

Appendix A

Since the equation of state is used in the derivation of the governing equation for the volume fraction, changing the equation of state might result in slight changes in the volume fraction equation. Here, the derivation of observable volume fraction equation is presented for a more general equation of state, the stiffened equation of state. In this case the derived equation is the same as the one derived for ideal gas equation of state.

The stiffened equation of state can be written as

[TABLE]

in which $P^{\infty}$ is usually zero for gases and non-zero for liquids. To derive the observable volume fraction equation, we insert Eq. 27 into Eq. 16 and by rearranging we have:

[TABLE]

With the pressure equilibrium requirement, the first parenthesis is zero. Since this equation needs to be satisfied for any pressure, the terms in the second and third parentheses should be zero. By defining the mixture property $\Pi$ similar to Eq. 11, $\Pi=\frac{{\gamma{p^{\infty}}}}{{\gamma-1}}=\frac{{{\gamma_{1}}p_{1}^{\infty}{z_{1}}}}{{{\gamma_{1}}-1}}+\frac{{{\gamma_{2}}p_{2}^{\infty}{z_{2}}}}{{{\gamma_{2}}-1}}$, and additionally using Eq. 11, the simplified observable equation for the volume fraction becomes Eq. 18, $\frac{{\partial{z_{1}}}}{{\partial t}}+\overline{\bf{u}}\cdot\nabla{z_{1}}=0.$.

Here, the observable equations are tested using a gas-liquid Riemann problem used to model underwater explosions [55, 12, 38]. The location of the discontinuity is at the origin at time equal zero and the initial condition is:

[TABLE]

Figure 11 shows the result of this problem by solving the two-phase observable Euler equations with a grid resolution of $\Delta x=1/20$ or 200 points across physical domain and an observability limit $\alpha=\Delta x$. As shown in the figure, the observable method produce good results and captures all the wave speeds correctly.