*4.2.4 Endoscopic retrograde cholangiography*

The rupture of a hydatid cyst in the intrahepatic bile duct can initiate some complications, which might become serious mainly due to the development of

**Figure 5.** *Unilocular cyst (CE1, white arrow). Detached membrane. Cyst (CE2, white star).*

**Figure 6.** *Contrast-enhanced CT. Septated cyst (white arrow).*

**Figure 7.** *Daughter vesicles in mother cyst (white arrow).*

**Figure 8.** *Coronal contrast-enhanced CT. CE3 A, water-lily sign (white arrow).*

**157**

**Figure 12.**

treatment [40].

cholangitis and septicemia. In these cases, the ERC provides the diagnosis and performs the removal of hydatid material, with the objective of improving the general conditions of the patient before carrying out the definitive surgical

*MRI T2 axial. Detached membranes. Hypointense pericyst (white arrow). Water-lily sign.*

*(A) US content, heterogeneous with a solid appearance (CE4). (B) In the same patient, the contrast-enhanced*

*Surgical Treatment of Hepatic Hydatidosis DOI: http://dx.doi.org/10.5772/intechopen.86319*

*Dome location with annular calcification (CE 5, white arrow).*

**Figure 10.**

**Figure 11.**

*CT diagnoses an unilocular hydatid cyst (CE1).*

**Figure 9.** *Central matrix. Daughter vesicles. Cartwheel sign (white arrow).*

*Surgical Treatment of Hepatic Hydatidosis DOI: http://dx.doi.org/10.5772/intechopen.86319*

*Liver Disease and Surgery*

**Figure 7.**

**Figure 8.**

*Daughter vesicles in mother cyst (white arrow).*

*Coronal contrast-enhanced CT. CE3 A, water-lily sign (white arrow).*

*Central matrix. Daughter vesicles. Cartwheel sign (white arrow).*

**156**

**Figure 9.**

#### **Figure 10.** *Dome location with annular calcification (CE 5, white arrow).*

#### **Figure 11.**

*(A) US content, heterogeneous with a solid appearance (CE4). (B) In the same patient, the contrast-enhanced CT diagnoses an unilocular hydatid cyst (CE1).*

cholangitis and septicemia. In these cases, the ERC provides the diagnosis and performs the removal of hydatid material, with the objective of improving the general conditions of the patient before carrying out the definitive surgical treatment [40].

**Figure 13.**

*MRI T2 axial. Detached membranes (white star). Hydatid membranes in bile duct (white arrow).*

**Figure 14.** *MRI coronal. Detached membranes. Bile duct with membranes (white arrow).*

**159**

*Surgical Treatment of Hepatic Hydatidosis DOI: http://dx.doi.org/10.5772/intechopen.86319*

Currently, diagnosis and follow-up of patients with cystic echinococcosis are achieved especially through imaging. Serology is used for the same purpose, using the detection of IgG-specific antigens. However, low sensitivity and specificity rates have been reported. In addition, false positives appear during follow-up due to the persistence of antibodies over time. There is a lot of research (recombined proteins, isotopic antibodies, subisotopic IgG, synthetic peptides), which seek to develop new antibodies by means of the molecular technique, allowing a better diagnosis of

Frequently, there are complications secondary to the hepatic location of the cyst or by involvement of adjacent organs, in which case symptoms and signs of greater

• More intense pain might appear when the Glisson capsule is stressed by larger cysts. In these cases, an abdominal mass mobilizing together with the respiratory movements is visualized during the physical examination. These giant cysts are easily visualized with the ultrasound. CT allows to define in more detail the elevated right hemidiaphragm and the development of secondary

• Another complication of large cysts occurs when there is an intrahepatic breach or rupture of the cyst to the peritoneum either spontaneously or by trauma. The discharge of fertile hydatid fluid in the liver or peritoneal cavity causes anaphylaxis of a different magnitude and new hydatid implants [43].

• Cyst rupture and emptying of hydatid fluid or membranes into the biliary tree lead to obstructive jaundice, many times accompanied by severe acute cholangitis [44, 45]. The cysts can become infected as well and determine the formation of liver abscesses, which can sometimes lead to septicemia. In these cases, CT and MRC allow to achieve a better definition of the characteristics of the abscess and whether there is emptying of hydatid material into the

• The chronic inflammatory process of hepatic cysts located in segments of the liver dome determines firm adhesions to the right hemidiaphragm and even transits toward the pleuropulmonary space. As a result of this transphrenic transit, patients may present pleural empyemas or bronchopneumonia [46]. CT and MRC help to achieve a better definition and provide the most appropri-

• There are occasions in which large cysts may be more complicated due to the simultaneous rupture and emptying of hydatid material to the biliary tree and bronchi. A bilio-pleuro-bronchial fistula is established with the occurrence of the pathognomonic sign of bilioptysis, i.e., the expectoration of the bile. These patients present a fairly severe septic episode with hepatic and respiratory functional compromise [47]. The thorax-abdominal CT and MRC are useful to

These complications are better visualized with the use of CT.

intensity will appear. Among them, the following should be mentioned:

**4.3 Serology**

this parasitosis [41, 42].

pulmonary basal atelectasis.

**5. Complications**

biliary tree.

ate and safe management.

diagnose this complication [39].

**Figure 15.** *MR cholangiography. Residual cyst membranes (white star). Cholangiohydatidosis.*


*Article* **SARS-CoV-2 Seroprevalence Post-First Wave among Primary Care Physicians in Catania (Italy)**

**Caterina Ledda 1,\*, Flavia Carrasi 1, Maria Teresa Longombardo 2, Gianluca Paravizzini 2 and Venerando Rapisarda 1**

 **\*** Correspondence<EMAIL_ADDRESS>Tel./Fax: +39-095-3782049

**Abstract:** Family physicians or pediatricians and general practitioners (GPs) work in non-hospital settings. GPs usually visit many patients, frequently at their homes, with low potential, if any, to control the work setting. Particularly during the initial phases of the COVID-19 outbreak, they were not informed about the occurrence of SARS-CoV-2-infected patients, with inadequate information regarding the risk, a lack of suitable protective measures and, in some cases, deficient or poor accessibility to personal protective equipment (PPE). During the first wave of COVID-19, primary care physicians were on the front line and isolated the first cases of the disease. The present study aims to estimate the seroprevalence of SARS-CoV-2 in a cohort of 133 GPs working in Catania (Italy) after the first wave of COVID-19. Serological analysis revealed a low seroprevalence (3%) among GPs. The low seroprevalence highlighted in the results can be attributed to correct management of patients by GPs in the first wave. It is now hoped that mass vaccination, combined with appropriate behavior and use of PPE, can help further reduce the risk of COVID-19 disease.

**Keywords:** COVID-19; healthcare worker; healthcare personnel; SARS-CoV-2; general practitioners; medical doctor


local phenomenon, it has been shown to be of importance, especially in the transition regime. The physical mechanism of most of the widely used SGS models is the forward scatter of energy; in other words, the SGS models are absolutely dissipative, such as the eddy viscosity models (Smagorinsky and Vreman SGS models). In order to account for the backscatter of energy, several modifications to the eddy viscosity models have been proposed. An improvement is to calculate the eddy viscosity coefficient dynamically which is a function of space and

The dynamic Smagorinsky model for compressible turbulence is as follows, the model coeffi-

CD ¼ 〈LijMij〉

� 2Δ2 ρ

� � 2

The stretched-vortex model is based on an explicit structural modeling of small-scale dynamics. It can simulate the multiscale compressible turbulence and allows the anisotropy of the flow to be extended to the dissipation scale. In the stretched-vortex model, the flow within a computational grid cell is assumed to result from an ensemble of straight, nearly axisymmetric vortices aligned

S � � � � �

� � � � �

where Lij is the Leonard stress, 〈 〉 indicates the statistics averaging.

66 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

� � � � �

Lij ¼ b ρu~i u~j � 1 ρ^ b ρu~i b ρu~j � �

^ Sij � δij 3 ^ Skk � �

the overbars of "¯" and "^" denote the grid filter and test filter, respectively.

with the local resolved scale strain or vorticity. The resulting SGS flux terms are [16]

Δc 2

subgrid vortex axis, kc ¼ π=Δc is the cutoff wave number, E kð Þ¼ K0ε2=3k�5=3 exp �2k

τij ¼ ρ~k δij � e

~k½ δij � e

~k½ δij � e

represents the energy spectrum of subgrid vortices, and K0 is the Kolmogorov prefactor, ε is the

ν i e ν j � �

ν i e ν j � � ∂ ~cpT~ � �

ν i e ν j � � ∂Y~

∂xj

∂xj

〈MijMij〉 ð19Þ

3

ν is the unit vector aligned with the

2 ν=ð Þ 3j~aj � �

5 ð21Þ

b

Sij � δij 3 Skk ð20Þ

ð22Þ

ð23Þ

ð24Þ

time and can be negative in some regions [15].

Mij ¼ � 2^

2.2.4. Stretched-vortex model

where ~k ¼

ð∞ kc Δ 2 ρ^ ^ S � � � � �

> QT i ¼ �ρ

> > QY i ¼ � Δc 2

E kð Þdk is the subgrid turbulent kinetic energy, e

4

cient is

$$\begin{aligned} \frac{\partial \rho}{\partial t} &= 0\\ \frac{\partial \rho u\_i}{\partial t} &= \frac{\partial \sigma\_{ij}}{\partial x\_j} - \frac{\partial \tau\_{ij}}{\partial x\_j} \\ \frac{\partial \rho E}{\partial t} &= \frac{\partial (\sigma\_{ij} - \tau\_{ij})u\_i}{\partial x\_j} - \frac{\partial \left(q\_j^l + Q\_j^T\right)}{\partial x\_j} \\ \frac{\partial \chi^{(s)}}{\partial t} &= \frac{\partial}{\partial x\_j} \left(D \frac{\partial \chi^{(s)}}{\partial x\_j}\right) - \frac{\partial Q\_j^T}{\partial x\_j} \quad s = 1, 2, \dots, N - 1 \end{aligned} \tag{26}$$

For the inviscid flux, the three-dimensional problem can be simplified into one-dimensional (1D) problem in three directions of x, y, and z by the dimension splitting technique,

$$\begin{aligned} \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} &= 0 \\ \frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} &= 0 \\ \frac{\partial (\rho E)}{\partial t} + \frac{\partial (\rho u E + p u)}{\partial x} &= 0 \\ \frac{\partial \rho Y(s)}{\partial t} + \mu \frac{\partial Y(s)}{\partial x} &= 0 \qquad s = 1, 2, \dots, N - 1 \end{aligned} \tag{27}$$

And then the one-dimensional Eq. (27) in each direction is resolved by the two-step Lagrange-Remapping algorithm. Also, one time step is divided into four substeps: (i) the piecewise parabolic interpolating of physical quantities, (ii) solving the Riemann problems approximately, (iii) marching of the Lagrange equations, and (iv) remapping the physical quantities back to the stationary Euler meshes. The orders of accuracy of the spatial and temporal schemes are the third and second, respectively, for smooth flows.

#### 2.3.1. Lagrange step of finite volume method

The Lagrange equations in 1D for multi-fluid can be written as

$$\begin{cases} \frac{\partial \tau}{\partial t} - \frac{\partial (r^a u)}{\partial m} = 0 \\ \frac{\partial u}{\partial t} + r^a \frac{\partial p}{\partial m} = 0 \\ \frac{\partial E}{\partial t} + \frac{\partial (r^a u p)}{\partial m} = 0 \\ \frac{\partial Y^{(s)}}{\partial t} = 0 \qquad s = 1, 2, \dots, N - 1 \end{cases} \tag{28}$$

where τ is the specific volume, α is 0, 1, and 2 corresponding to the plane, axial symmetry, and spherical symmetry problems, r is the Lagrange spatial coordinates, m is the Lagrange mass coordinates m ¼ ðr r0 ρr αdr, u ¼ ∂r ∂t , mj�1/2 and mjþ1/2 are the mass coordinates of both sides of the jth grid, Δm ¼ mjþ1/2 � mj�1/2. The mass average of physical quantities τ, u, E, Yð Þs � �n j for the jth grid can be defined as

$$
\begin{pmatrix} \tau \\ u \\ E \\ Y(s) \end{pmatrix}\_j^n = \frac{1}{\Delta m\_j} \int\_{m\_{j-1/2}}^{m\_{j+1/2}} \begin{pmatrix} \tau(m\_\prime \ t\_n) \\ u(m\_\prime \ t\_n) \\ E(m\_\prime \ t\_n) \\ Y^{(s)}(m\_\prime \ t\_n) \end{pmatrix} dm \tag{29}
$$

Because of the calculations are referred to as scalars, here the averaged physical quantities in a cell are written as a uniformity Qn j at time t. The piecewise parabolic function Qn j ð Þ m in cells are constructed to compute the time average of physical quantity Q on both sides of the grid edge mjþ1/2,Q~ jþ1=2,L and Q~ jþ1=2,R. Then the Riemann problem at the grid edge mjþ1/2 is solved by using double shock wave approximation. After Lagrange marching step, we can obtain distribution of the physical quantities τnþ1 j n o, unþ1 j n o, Enþ1 j n o, and the position of new grids r nþ1 j�1=2 n o at time tnþ1, the pressure pnþ1 j n o can be computed by equation of state, the Yð Þs � �nþ1 j � � will not change. The marching formula is as follows

$$\begin{cases} r\_{j+1/2}^{n+1} = r\_{j+1/2}^{n} + \tilde{u}\_{j+1/2} \Delta t \\\\ \tau\_{j}^{n+1} = \frac{\left(r\_{j+1/2}^{n+1}\right)^{\alpha+1} - \left(r\_{j-1/2}^{n+1}\right)^{\alpha+1}}{(\alpha+1)\Delta m\_{j}} \\\\ u\_{j}^{n+1} = u\_{j}^{n} - \frac{\Lambda}{2} \frac{\Delta t}{\Delta m\_{j}} \left(\sigma\_{j+1/2} + \sigma\_{j-1/2}\right) \left(\tilde{p}\_{j+1/2} - \tilde{p}\_{j-1/2}\right) \\\\ E\_{j}^{n+1} = E\_{j}^{n} - \frac{\Lambda t}{\Delta m\_{j}} \left(\sigma\_{j+1/2}\tilde{u}\_{j+1/2}\tilde{p}\_{j+1/2} - \sigma\_{j-1/2}\tilde{u}\_{j-1/2}\tilde{p}\_{j-1/2}\right) \\\\ \sigma\_{j+1/2} = \frac{\left(r\_{j+1/2}^{n+1}\right)^{\alpha+1} - \left(r\_{j+1/2}^{n}\right)^{\alpha+1}}{(\alpha+1)\left(r\_{j+1/2}^{n+1} - r\_{j+1/2}^{n}\right)} \end{cases} \tag{31}$$

#### 2.3.2. Remap step

ð28Þ

j

j ð Þ m in cells are

dm ð29Þ

, and the position of new grids

can be computed by equation of state, the

for the

And then the one-dimensional Eq. (27) in each direction is resolved by the two-step Lagrange-Remapping algorithm. Also, one time step is divided into four substeps: (i) the piecewise parabolic interpolating of physical quantities, (ii) solving the Riemann problems approximately, (iii) marching of the Lagrange equations, and (iv) remapping the physical quantities back to the stationary Euler meshes. The orders of accuracy of the spatial and temporal

∂t ¼ 0 s ¼ 1, 2, …, N � 1

, mj�1/2 and mjþ1/2 are the mass coordinates of both sides of the

1

CCA

where τ is the specific volume, α is 0, 1, and 2 corresponding to the plane, axial symmetry, and spherical symmetry problems, r is the Lagrange spatial coordinates, m is the Lagrange mass

ð mjþ1=2

0

BB@

τð Þ m, tn u m, t ð Þn E m, t ð Þn YðsÞ

j at time t. The piecewise parabolic function Qn

, Enþ1 j n o

ð Þ m, tn

mj�1=2

Because of the calculations are referred to as scalars, here the averaged physical quantities in a

constructed to compute the time average of physical quantity Q on both sides of the grid edge mjþ1/2,Q~ jþ1=2,L and Q~ jþ1=2,R. Then the Riemann problem at the grid edge mjþ1/2 is solved by using double shock wave approximation. After Lagrange marching step, we can obtain distri-

> , unþ1 j n o

jth grid, Δm ¼ mjþ1/2 � mj�1/2. The mass average of physical quantities τ, u, E, Yð Þs � �n

schemes are the third and second, respectively, for smooth flows.

68 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

The Lagrange equations in 1D for multi-fluid can be written as

8

>>>>>>>>>>>><

>>>>>>>>>>>>:

αdr, u ¼ ∂r ∂t

> τ u E Yð Þs

1

n

¼ 1 Δmj

CCA

j

j n o

> j n o

will not change. The marching formula is as follows

0

BB@

at time tnþ1, the pressure pnþ1

∂τ ∂t � ∂ r a ð Þ u ∂m ¼ 0

∂u ∂t þ r a ∂p ∂m ¼ 0

∂E ∂t þ ∂ r a ð Þ up ∂m ¼ 0

∂Yð Þs

2.3.1. Lagrange step of finite volume method

coordinates m ¼

r nþ1 j�1=2 n o

Yð Þs � �nþ1 j � �

jth grid can be defined as

ðr r0 ρr

cell are written as a uniformity Qn

bution of the physical quantities τnþ1

In Lagrange step, the computational cells distort to follow the material motion, and in Remap step, the averaged physical quantities at the distorted Lagrange cells are remapped back to the stationary Euler meshes. The piecewise parabolic interpolation and integral methods in this step are the same as the ones in Lagrange step.

After Remap step, we define the physical quantities in Euler mesh as ρEuler � �nþ1 j , uð Þ Euler nþ1 j , ð Þ EEuler nþ1 j , Yð Þs � � Euler � �nþ1 j , the volume of fluid through the grid boundary j þ 1=2 as ΔV� jþ1=2, and the average of physical quantities as ρ� , u� , p� , Yð Þs � � � �� jþ1=2 . The density after Lagrange step is ρnþ1 j ¼ 1=τnþ1 j . In Remap step

$$
\Delta m\_{\rangle}^{n+1} = \rho\_{\rangle}^{n+1} V\_{\rangle}^{n+1} \tag{32}
$$

$$
\Delta m\_{\rangle}^{n} = \left(\rho\_{\text{Euler}}\right)\_{\rangle}^{n} V\_{\text{j}}^{n} \tag{33}
$$

$$
\Delta m\_{\dot{\jmath}+1/2}^{\*} = \rho\_{\dot{\jmath}+1/2}^{\*} \Delta V\_{\dot{\jmath}+1/2}^{\*}\tag{34}
$$

$$
\Delta E\_{j+1/2}^\* = \frac{1}{2} \left( u\_{j+1/2}^\* \right)^2 \Delta m\_{j+1/2}^\* + \Delta e\_{j+1/2}^\* \tag{35}
$$

$$\left(\rho\_{\text{Euler}}\right)\_{j}^{n+1} = \left[\Delta m\_{j}^{n+1} - \left(\rho\_{j+1/2}^{\*}\Delta V\_{j+1/2}^{\*} - \rho\_{j-1/2}^{\*}\Delta V\_{j-1/2}^{\*}\right)\right] \Big/ V\_{j}^{n} \tag{36}$$

$$\left(\left(\mu\_{\text{Euler}}\right)\_{j}^{n+1} = \left[\mathfrak{u}\_{j}^{n+1}\Delta m\_{j}^{n+1} - \left(\mathfrak{u}\_{j+1/2}^{\*}\Delta m\_{j+1/2}^{\*} - \mathfrak{u}\_{j-1/2}^{\*}\Delta m\_{j-1/2}^{\*}\right)\right] \Big/\Delta m\_{j}^{n} \tag{37}$$

$$\left( \left( Y^{\{\mathbf{s}\}} \right)\_{\text{Euler}} \right)\_{j}^{n+1} = \left[ \left( Y^{\{\mathbf{s}\}} \right)\_{j}^{n+1} V\_{j}^{n+1} - \left( \left( Y^{\{\mathbf{s}\}} \right)\_{j+1/2}^{\*} \Delta V\_{j+1/2}^{\*} - \left( Y^{\{\mathbf{s}\}} \right)\_{j-1/2}^{\*} \Delta V\_{j-1/2}^{\*} \right) \right] \Big/ V\_{j}^{n} \tag{38}$$

$$\left(\left(E\_{\text{Euler}}\right)\_{j}^{n+1} = \left[E\_{j}^{n+1}\Delta m\_{j}^{n+1} - \left(\Delta E\_{j+1/2}^{\*} - \Delta E\_{j-1/2}^{\*}\right)\right] / \Delta m\_{j}^{n} \tag{39}$$

where Δe� jþ1=2 is the advection of specific energy through the grid boundary j þ 1=2.

#### 2.3.3. Viscous flux, heat flux, and scalar flux

The viscous flux, heat flux, and scalar flux of Eq. (26) are calculated based on the computed inviscid flux by using second-order spatial center difference and two-step Runge-Kutta time marching. The first equation of Eq. (26) can be neglected, and which is written in conserved form

$$\frac{\partial \overline{\mathbf{U}}}{\partial t} + \frac{\partial \overline{\mathbf{F}}}{\partial x} + \frac{\partial \overline{\mathbf{G}}}{\partial y} + \frac{\partial \overline{\mathbf{H}}}{\partial z} = \mathbf{0} \tag{40}$$

and

$$\begin{aligned} \mathbf{U} &= \left(\rho u, \rho v, \rho w, \rho \mathbf{E}\right)^{T} \\ \mathbf{\overline{F}} &= \left(-\sigma\_{xx}, -\sigma\_{yx}, -\sigma\_{zx}, -u\sigma\_{xx} - v\sigma\_{yx} - w\sigma\_{zx} + q\_{x}\right)^{T} \\ \mathbf{\overline{G}} &= \left(-\sigma\_{xy}, -\sigma\_{yy}, -\sigma\_{zy}, -u\sigma\_{xy} - v\sigma\_{yy} - w\sigma\_{zy} + q\_{y}\right)^{T} \\ \mathbf{\overline{H}} &= \left(-\sigma\_{xz}, -\sigma\_{yz}, -\sigma\_{zz}, -u\sigma\_{xz} - v\sigma\_{yz} - w\sigma\_{zz} + q\_{z}\right)^{T} \end{aligned} \tag{41}$$

In the Cartesian coordinate system, the spatial derivative terms of Eq. (40) can be dispersed as follows:

$$\mathbf{L}\_{h}\left(\overline{\mathbf{U}}\_{i,j,k}^{n}\right) = \frac{\mathbf{F}\_{i-1/2,j,k}^{n} - \mathbf{F}\_{i+1/2,j,k}^{n}}{\Delta x} + \frac{\overline{\mathbf{G}}\_{i,j-1/2,k}^{n} - \overline{\mathbf{G}}\_{i,j+1/2,k}^{n}}{\Delta y} + \frac{\mathbf{H}\_{i,j,k-1/2}^{n} - \mathbf{H}\_{i,j,k+1/2}^{n}}{\Delta z} \tag{42}$$

Then, Eq. (40) is solved by two Runge-Kutta time marching,

$$\begin{cases} \overline{\mathbf{U}}\_{i,j,k}^{(l)} = \overline{\mathbf{U}}\_{i,j,k}^{\mathrm{E}} + \Delta t \mathbf{L}\_{h} \left( \mathbf{U}\_{i,j,k}^{n} \right) \\ \mathbf{U}\_{i,j,k}^{n+1} = \frac{1}{2} \left[ \overline{\mathbf{U}}\_{i,j,k}^{\mathrm{E}} + \overline{\mathbf{U}}\_{i,j,k}^{(l)} + \Delta t \mathbf{L}\_{h} \left( \overline{\mathbf{U}}\_{i,j,k}^{(l)} \right) \right] \end{cases} \tag{43}$$

The detailed descriptions of numerical algorithm are referred to Ref. [17].

#### 3. Verification and validation

In this section, the validity and reliability of our compressible multi-fluid dynamic programs are to be verified and validated by comparisons with analytical solutions and hydrodynamic interface instability experiments in a shock tube.

Figure 2. Comparison of postshock density (a), pressure (b), and velocity (c) in copper at t ¼ 1.0 μs when a copper pellet collides with a copper target at three different velocities.

## 3.1. Riemann problem of condensed matter

EEuler � �nþ1

2.3.3. Viscous flux, heat flux, and scalar flux

where Δe�

form

and

follows:

Lh Un i,j, k � � ¼ Fn

3. Verification and validation

interface instability experiments in a shock tube.

i�1=2,j, k � Fn

Δx

Then, Eq. (40) is solved by two Runge-Kutta time marching,

UðlÞ i,j, k ¼ UE

8 ><

>:

Unþ1 i,j, k ¼ 1

iþ1=2,j, k

þ Gn

2 UE

The detailed descriptions of numerical algorithm are referred to Ref. [17].

i,j, k þ ΔtLh Un

i,j, k þ UðlÞ

In this section, the validity and reliability of our compressible multi-fluid dynamic programs are to be verified and validated by comparisons with analytical solutions and hydrodynamic

j ¼ Enþ1

70 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

j Δmnþ1

∂U ∂t þ

U ¼ ρu,ρv,ρw,ρE � �T

∂F ∂x þ ∂G ∂y þ ∂H

F ¼ �σxx, � σyx, � σzx, � uσxx � vσyx � wσzx þ qx � �T

G ¼ �σxy, � σyy, � σzy, � uσxy � vσyy � wσzy þ qy � �T

H ¼ �σxz, � σyz, � σzz, � uσxz � vσyz � wσzz þ qz � �T

In the Cartesian coordinate system, the spatial derivative terms of Eq. (40) can be dispersed as

i,j�1=2, k � Gn

Δy

i,j, k � �

h i � �

i,j, k þ ΔtLh UðlÞ

i,jþ1=2, k

þ Hn

i,j, k

j � ΔE�

jþ1=2 is the advection of specific energy through the grid boundary j þ 1=2.

The viscous flux, heat flux, and scalar flux of Eq. (26) are calculated based on the computed inviscid flux by using second-order spatial center difference and two-step Runge-Kutta time marching. The first equation of Eq. (26) can be neglected, and which is written in conserved

h i � �

jþ1=2 � ΔE�

j�1=2

=Δmn

∂z ¼ 0 ð40Þ

i,j, k�1=2 � Hn

i,j, kþ1=2 Δz ð42Þ

j ð39Þ

ð41Þ

ð43Þ

A copper pellet collides with a copper target with three velocities of 2, 4, and 8 km/s. Using the reduction form of Gruneisen EOS for copper in simulations, the one-dimensional numerical results of postshock density (a), pressure (b), and velocity (c) compare with theoretical solutions, as shown in Figure 2, the solid lines corresponding to numerical results and the dot lines corresponding to the theoretical solutions.

## 3.2. Riemann problem of one-dimensional gas/liquid

At initial time, the region [0, 1.0 cm] is filled with gas with high pressure 1.0 � 108 Pa and density 1.29 g/cm3 . The region [1.0 cm, 5.0 cm] is filled with water with pressure 1.0 � 105 Pa and density 1.0 g/cm3 . The gas and water are all described by Stiffen gas EOS. The left and right boundaries are flow. Figure 3 shows the distributions of the density (a), pressure (b), and velocity (c) at 20 μs. There are a forward shock wave and a backward rarefaction wave after interaction. The pressure and velocity around the interface are well continuous and have no nonphysical oscillation.

## 3.3. Single-mode Richtmyer-Meshkov instability

The two- and three-dimensional single-mode Richtmyer-Meshkov instabilities are numerical simulated by MVPPM program, which also compare with the theoretical model [18]. The initial small perturbation is a sinusoidal one with wavelength 60 mm (global wavelength for a three-dimensional case) and amplitude 1.0 mm. The incident air shock wave with Mach number 1.2 impacts the air/SF6 single-mode interface. Figure 4 shows the comparisons of amplitude of single-mode perturbed interface with the linear impulsive and nonlinear models, the left one corresponding to the two-dimensional (2D) results, and the right one corresponding to the three-dimensional results.

#### 3.4. Simulations of shock tube experiments

The compressible multi-fluid dynamic programs are used to simulated several shock tube experiments of interface instability for validation. These experiments include planar interface and gas cylinder shock tube experiment and planar and cylindrical jelly experiments.

Figure 5 shows the comparison of 2D calculated width of turbulent mixing zone (TMZ) [19] of Leinov's shock tube experiment with reshock with experiment [20] in which the Mach number

Figure 3. Distributions of the density (a), pressure (b), and velocity (c) at 20 μs for one-dimensional Riemann problem of gas/liquid interface.

Figure 4. Comparisons of amplitude of single-mode perturbed interface with the linear impulsive model and nonlinear models, left: two-dimension and right: three-dimension.

Figure 5. TMZ width versus time (t ¼ 0 denotes the reshock arrival to the interface).

of incident air shock wave is 1.2, which impacts the air/SF6 interface, the grid size is 50 μm. Figure 6 shows the 3D calculated TMZ width [21] of Poggi's multi-mode shock tube experiment [22] with reshock in which the Mach number of incident SF6 shock wave is 1.453, which impacts the SF6/air interface. Figures 7 [10] and 8 [23] show the calculated and experimental interface images of AWE's SF6 half-height and double-bump shock tube experiment [24, 25], the Mach number of incident air shock wave is 1.26, respectively. Figure 9 [26] shows the SF6 gas cylinder evolution at different times under the initial air shock wave with the Mach number 1.2, Figure 10 shows the width and height of gas cylinder at different times for experiment and numerical simulations.

Figure 6. Time history of TMZ width.

Figure 3. Distributions of the density (a), pressure (b), and velocity (c) at 20 μs for one-dimensional Riemann problem of

72 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 4. Comparisons of amplitude of single-mode perturbed interface with the linear impulsive model and nonlinear

gas/liquid interface.

models, left: two-dimension and right: three-dimension.

Figure 5. TMZ width versus time (t ¼ 0 denotes the reshock arrival to the interface).

Figure 7. Two dimension calculated from MVFT and experimental half-height interface at different times.

Figure 8. Three dimension calculated from MVFT and experimental double-bump interface at different times, left column corresponding to experimental images, right column corresponding to 3D calculated images, and middle column corresponding to span-wise average images of 3D calculated results.

Figure 9. Evolutions of SF6 gas cylinder, upper row corresponding to experiment, lower row corresponding to numerical simulation by MVFT program, the time sequences are 0, 50, 190, 330, 470, 610, and 750 μs severally.

For the problems of interface instability with a high density ratio, our hydrodynamic programs are also applicable. For the LLNL's jelly ring experiments [27], the jelly ring only has a sinusoidal periodic initial perturbation at the outer interface and is driven by expansion of the explosion products of a gaseous mixture of C2H2 and O2. The jelly is mainly made of water. The thickness of ring is 15 mm, the mode number of perturbation is 6 and 36, respectively, and the initial amplitude is 1.0 mm. Figure 11 shows the evolutions of the jelly ring at different times from simulated from LLNL' CALE program (a) and our MVPPM program (b), and the experimental image (c) at 600 μs [28]. Figure 12 shows the time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of the jelly ring simulated from CALE and

Figure 10. Width and height of gas cylinder at different times for experimental and numerical simulations.

Figure 11. Evolutions of jelly ring at different simulated from CALE (a: left group) and MVPPM (b: middle group) programs, and the experimental image at 600 μs (c).

Figure 12. Time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of jelly ring simulated from CALE and MVPPM programs.

For the problems of interface instability with a high density ratio, our hydrodynamic programs are also applicable. For the LLNL's jelly ring experiments [27], the jelly ring only has a sinusoidal periodic initial perturbation at the outer interface and is driven by expansion of the explosion products of a gaseous mixture of C2H2 and O2. The jelly is mainly made of water. The thickness of ring is 15 mm, the mode number of perturbation is 6 and 36, respectively, and the initial amplitude is 1.0 mm. Figure 11 shows the evolutions of the jelly ring at different times from simulated from LLNL' CALE program (a) and our MVPPM program (b), and the experimental image (c) at 600 μs [28]. Figure 12 shows the time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of the jelly ring simulated from CALE and

Figure 9. Evolutions of SF6 gas cylinder, upper row corresponding to experiment, lower row corresponding to numerical

simulation by MVFT program, the time sequences are 0, 50, 190, 330, 470, 610, and 750 μs severally.

Figure 8. Three dimension calculated from MVFT and experimental double-bump interface at different times, left column corresponding to experimental images, right column corresponding to 3D calculated images, and middle column

corresponding to span-wise average images of 3D calculated results.

74 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

MVPPM programs [28]. Figure 13 shows the images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE (b), and MVPPM (c) programs [28].

From the above comparisons between our numerical simulations and theoretical model, experiments for the problems of hydrodynamic interface instability from the density ratio low

Figure 13. Images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE program (b), and MVPPM program (c).

to the high density ratio, the numerical results agree well with theory and experiments, so the validity and reliability of our compressible multi-fluid dynamic programs have been verified and validated.

## 4. Interface instability and turbulent mixing

As we know, the physical mechanism for the occurrence of Richtmyer-Meshkov instability is the baroclinic vorticity deposition at the interface resulting from the misalignment of the pressure gradient across the shock front and the local density gradient across the interface. The evolution equation of vorticity is as follows:

$$\frac{d\boldsymbol{\omega}}{dt} = \frac{\nabla\rho \times \nabla p}{\rho^2} + \boldsymbol{\omega} \cdot \nabla \mathbf{u} - \boldsymbol{\omega} \nabla \cdot \mathbf{u} \tag{44}$$

where ω ¼ ∇ � u is the vorticity and viscous terms are neglected. The first term on the right side of Eq. (44) is the baroclinic vorticity production term. The second term on the right side of Eq. (44) is the vortex-stretching term, which is zero in the two-dimensional case, as the vorticity and velocity fields are orthogonal. This term enhances dissipation, resulting in more diffuse and smaller scale structures in the turbulent mixing zone. The third term on the right side is the compression term and does not contribute to the vorticity evolution significantly. The baroclinic vorticity production is much larger when the shock wave impacts the interface and pass through it and constitutes the principal mechanism of the Richtmyer-Meshkov instability.

In this section, we will introduce the growth laws of the Richtmyer-Meshkov instability and the induced turbulent mixing and its dynamic behavior by numerical simulations.

#### 4.1. Growth laws of Richtmyer-Meshkov instability and turbulent mixing

Figure 14 shows the amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock for initial perturbation amplitude 1 and 12 mm

Figure 14. Amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock.

Figure 15. Schematic of shock tube and computational model.

to the high density ratio, the numerical results agree well with theory and experiments, so the validity and reliability of our compressible multi-fluid dynamic programs have been verified

Figure 13. Images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE program (b), and MVPPM

As we know, the physical mechanism for the occurrence of Richtmyer-Meshkov instability is the baroclinic vorticity deposition at the interface resulting from the misalignment of the pressure gradient across the shock front and the local density gradient across the interface.

where ω ¼ ∇ � u is the vorticity and viscous terms are neglected. The first term on the right side of Eq. (44) is the baroclinic vorticity production term. The second term on the right side of Eq. (44) is the vortex-stretching term, which is zero in the two-dimensional case, as the vorticity and velocity fields are orthogonal. This term enhances dissipation, resulting in more diffuse and smaller scale structures in the turbulent mixing zone. The third term on the right side is the compression term and does not contribute to the vorticity evolution significantly. The baroclinic vorticity production is much larger when the shock wave impacts the interface and pass through it and constitutes the principal mechanism of the Richtmyer-Meshkov

In this section, we will introduce the growth laws of the Richtmyer-Meshkov instability and

Figure 14 shows the amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock for initial perturbation amplitude 1 and 12 mm

the induced turbulent mixing and its dynamic behavior by numerical simulations.

4.1. Growth laws of Richtmyer-Meshkov instability and turbulent mixing

ρ2 þ ω � ∇u � ω∇ � u ð44Þ

4. Interface instability and turbulent mixing

dω

76 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

dt ¼ ∇ρ � ∇p

The evolution equation of vorticity is as follows:

and validated.

program (c).

instability.

and wavelength 60 mm. After initial shock, the perturbation enters the nonlinear stage quickly. The growth rate increases fast and reaches the highest peak, then will reduce owing to the reduction of the effect of compressibility and the dominant role of flow nonlinearity. At the late times, the amplitude is increasing linearly, and the growth rate remains a constant. For the larger initial amplitude, the perturbation enters the later linear growth stage earlier.

For the multi-mode Richtmyer-Meshkov instability and the induced turbulent mixing with reshock, the initial air shock wave with Mach number 1.2 impacts the air/SF6 interface. The multi-mode initial perturbation is composed of eight dominant mode wavelengths of 0.8, 1.0, 1.25, 1.6, 2.0, 2.5, 3.2, and 4.0 mm superimposed with a very small random disturbance. The shock tube experiment can be referred to Ref. [29]. Figure 15 shows the schematic of shock tube and computational model. The transmitted wave rebounds between the interface and the end-wall of shock tube and produces multiple shock-interface interactions. At about 1.7 ms,

Figure 16. Time history of the TMZ width.

the transmitted shock wave reflected back off the end-wall impacts the interface. Figure 16 shows the time history of TMZ width, the black line denotes the numerically calculated width and the red line is obtained from fitting of the numerical results. As can be seen, after the initial shock, the TMZ width starts to grow as a power law ~t θ with the value θ to be determined as 0.352. After the reshock, more energy is deposited onto the interface to promote the development of TMZ, and the TMZ width evolves in time as a negative exponential law ~�e�t=t � where the value of t \* is 0.519. Then after the following interaction of the reflected rarefaction wave with the interface, the TMZ also evolves as a negative exponential law but with a different factor t \* ¼ 0.875. Under the subsequent impingements with lower and lower intensity, the TMZ width, after a slight reduction caused by the reflected compression wave, evolves in an approximate linear fashion with a velocity of 2.05 m/s. Figure 17 shows the instantaneous images of the TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times, the TMZ exhibits a very complex spatial structure.

#### 4.2. Evaluation of different subgrid-scale stress models

In large-eddy simulation, the effect of small scales on large-scale motions is represented by the SGS stress model. Most of the commonly used SGS models assume that the main function of subgrid scales is to remove energy from the large scales and dissipate it through the action of the viscous forces. But, as we know, in fact the energy is also transferred from the small scales to the large scales (backscatter) in a small and local range. The SGS turbulent dissipation, which is the work of SGS stress, represents the energy transfer between resolved and subgrid scales,

$$
\varepsilon\_{\rm scs} = \tau\_{\rm j\overline{}} \overline{\mathbf{S}}\_{\rm ij} \tag{45}
$$

If it is negative, the subgrid scales remove energy from the resolved scales (forward scatter); if it is positive, they release energy to the resolved scales (backscatter). It is easy to see that the eddy viscosity models such as Smagorinsky model, Vreman model, etc. are absolutely dissipative.

Figure 17. Instantaneous images of TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times.

the transmitted shock wave reflected back off the end-wall impacts the interface. Figure 16 shows the time history of TMZ width, the black line denotes the numerically calculated width and the red line is obtained from fitting of the numerical results. As can be seen, after the initial

0.352. After the reshock, more energy is deposited onto the interface to promote the development of TMZ, and the TMZ width evolves in time as a negative exponential law ~�e�t=t

wave with the interface, the TMZ also evolves as a negative exponential law but with a

the TMZ width, after a slight reduction caused by the reflected compression wave, evolves in an approximate linear fashion with a velocity of 2.05 m/s. Figure 17 shows the instantaneous images of the TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times, the TMZ exhibits a very complex spatial structure.

In large-eddy simulation, the effect of small scales on large-scale motions is represented by the SGS stress model. Most of the commonly used SGS models assume that the main function of subgrid scales is to remove energy from the large scales and dissipate it through the action of the viscous forces. But, as we know, in fact the energy is also transferred from the small scales to the large scales (backscatter) in a small and local range. The SGS turbulent dissipation, which is the work of SGS stress, represents the energy transfer between resolved and subgrid scales,

If it is negative, the subgrid scales remove energy from the resolved scales (forward scatter); if it is positive, they release energy to the resolved scales (backscatter). It is easy to see that the eddy viscosity models such as Smagorinsky model, Vreman model, etc. are absolutely

\* is 0.519. Then after the following interaction of the reflected rarefaction

\* ¼ 0.875. Under the subsequent impingements with lower and lower intensity,

θ with the value θ to be determined as

εSGS ¼ τijSij ð45Þ

�

shock, the TMZ width starts to grow as a power law ~t

78 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

4.2. Evaluation of different subgrid-scale stress models

where the value of t

Figure 16. Time history of the TMZ width.

different factor t

dissipative.

Figure 18. Distribution of calculated SGS turbulent dissipation in a streamwise direction for Smagorinsky and Vreman models and for AWE's SF6 half-height shock tube experiment [24] at two times. Red line corresponding to the Vreman model, green line corresponding to the Smagorinsky model.

Figure 18 shows the distribution of calculated SGS turbulent dissipation in streamwise direction for Smagorinsky and Vreman models and for AWE's SF6 half-height shock tube experiment [24] at two times [11]. The SGS dissipation of the Smagorinsky model is much greater than the Vreman model over 1.5 times; therefore, the dissipation is too great for the Smagorinsky SGS model. The SGS turbulent dissipations of the Vreman model (VM), the dynamic Smagorinsky model (DSM), and the stretched-vortex model (SVM) based on a planar Richtmyer-Meshkov instability with incident Mach number 1.2 are shown in Figure 19 [30]. Before the interfacial flow has developed to be turbulent completely, the dynamic and stretched-vortex models have all predicted the energy backscatter, but the energy backscatter predicted by the dynamic model is

Figure 19. Distribution of the SGS turbulent dissipation in a streamwise direction for the Vreman, dynamic Smagorinsky and stretched-vortex models at four times.

larger and its range is wider. After reshock, the turbulent fluctuations are stronger extremely the turbulent dissipation also increases. The dynamic model's dissipation is the highest, then followed by the Vreman model, and the stretched-vortex model's dissipation is the lowest. At the late time, the SGS dissipation of the dynamic and Vreman models is the same basically, and the dynamic model is also to be absolutely dissipative, yet the stretched-vortex model is still able to predict the local phenomenon of energy backscatter in a small range. So, the dynamic model is poor in representing the energy backscatter. The Vreman and stretched-vortex models are all robust, but the former is absolutely dissipative.

#### 4.3. Effects of the initial conditions on the growth of RMI and the turbulent mixing

#### 4.3.1. Effects of the initial conditions on the growth of single-mode RMI

First, we consider the effects of initial conditions of perturbation amplitude and wavelength on the growth of single-mode Richtmyer-Meshkov instability without reshock [18]. The incident shock wave with Mach number 1.2 impacts the air/SF6 interface. The initial perturbation amplitude and wavelength are listed in Table 1. The calculations are carried out in two dimensions, and the RMI does not develop into turbulent mixing completely. Figure 20 shows the effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ0 ¼ 60 mm. The perturbation amplitude (a) and growth rate (b) increase gradually with the increasing of initial amplitude. And when the initial amplitude is much larger, the growth of amplitude enters a linear stage earlier at the late times, and the growth rate remains a constant. The growth rate increases fast and reaches the highest peak at the early times. After the peak, the effect of compressibility is reduced and the flow nonlinearity starts to play a dominant role and causes the growth rate to decay with time. Figure 21 shows the effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial

Table 1. Initial conditions of perturbation for single-mode RMI.

larger and its range is wider. After reshock, the turbulent fluctuations are stronger extremely the turbulent dissipation also increases. The dynamic model's dissipation is the highest, then followed by the Vreman model, and the stretched-vortex model's dissipation is the lowest. At the late time, the SGS dissipation of the dynamic and Vreman models is the same basically, and the dynamic model is also to be absolutely dissipative, yet the stretched-vortex model is still able to predict the local phenomenon of energy backscatter in a small range. So, the dynamic model is poor in representing the energy backscatter. The Vreman and stretched-vortex models are all

Figure 19. Distribution of the SGS turbulent dissipation in a streamwise direction for the Vreman, dynamic Smagorinsky

4.3. Effects of the initial conditions on the growth of RMI and the turbulent mixing

First, we consider the effects of initial conditions of perturbation amplitude and wavelength on the growth of single-mode Richtmyer-Meshkov instability without reshock [18]. The incident shock wave with Mach number 1.2 impacts the air/SF6 interface. The initial perturbation amplitude and wavelength are listed in Table 1. The calculations are carried out in two dimensions, and the RMI does not develop into turbulent mixing completely. Figure 20 shows the effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ0 ¼ 60 mm. The perturbation amplitude (a) and growth rate (b) increase gradually with the increasing of initial amplitude. And when the initial amplitude is much larger, the growth of amplitude enters a linear stage earlier at the late times, and the growth rate remains a constant. The growth rate increases fast and reaches the highest peak at the early times. After the peak, the effect of compressibility is reduced and the flow nonlinearity starts to play a dominant role and causes the growth rate to decay with time. Figure 21 shows the effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial

4.3.1. Effects of the initial conditions on the growth of single-mode RMI

80 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

robust, but the former is absolutely dissipative.

and stretched-vortex models at four times.

Figure 20. Effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ0 ¼ 60 mm.

Figure 21. Effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial perturbation amplitude a0 ¼ 3 mm.

perturbation amplitude a0 ¼ 3 mm. The perturbation amplitude (a) and growth rate (b) reduce gradually at the early times and increase gradually at the late times with the increasing of initial wavelength. And when the initial perturbation strength (ration of initial amplitude to wavelength) is much smaller, the perturbation growth is mainly dependent on the initial perturbation amplitude and slightly dependent on the initial perturbation wavelength at the late times of RMI. The same conclusion can be obtained from three-dimensional calculations.

## 4.3.2. Effects of the initial conditions on the growth of multi-mode RMI and the induced turbulent mixing

For the effects of three-dimensional initial multi-mode conditions, the case of Richtmyer-Meshkov instability is same as the above Section 4.1. Table 2 lists the initial condition of perturbations, the perturbation strength (PS) is defined as the ratio of initial amplitude to wavelength. Turbulent kinetic energy K, dissipation rate ε, and enstrophy Ω are defined as follows [31]:

$$K = \frac{\langle \rho u''\_i u''\_i \rangle}{2 \langle \rho \rangle} + \frac{\langle \tau\_{ii} \rangle}{2 \langle \rho \rangle} \tag{46}$$

$$\varepsilon = \frac{\langle \sigma\_{\vec{\eta}} \eth \mathsf{u} \nu'' , \langle \mathsf{Ox}\_{\vec{\eta}} \rangle }{\langle \rho \rangle} - \frac{\langle \pi\_{\vec{\eta}} \eth \mathsf{u} \nu'' , \langle \mathsf{Ox}\_{\vec{\eta}} \rangle}{\langle \rho \rangle} \tag{47}$$

$$
\Omega = \frac{1}{2} \|\boldsymbol{\omega}\|^2 \tag{48}
$$

where u00 i is the velocity fluctuation and 〈 〉 denotes the transverse plane-average. For largeeddy simulations, the turbulent kinetic energy and dissipation rate both include two parts named as the resolved-scales (the first term on the right side of Eqs. (46) and (47)) and the subgrid-scales (the second term on the right side of Eqs. (46) and (47)).

Figure 22 shows the growth history of TMZ width. Figures 23 and 24 show the time histories of the perk values of turbulent kinetic energy and enstrophy, respectively. For the larger initial perturbation strength, the TMZ grows faster, the turbulent kinetic energy is also larger or the turbulence strength is also stronger, the deposited vorticity is larger too. The development of turbulent mixing has a strong dependence on the initial conditions between the initial shock and the impingement of the first reflected rarefaction wave, after that the evolution of the turbulent mixing has lost the memory of the initial conditions.

Table 2. Model parameters.

Figure 22. Growth history of the TMZ width for different models.

perturbation amplitude a0 ¼ 3 mm. The perturbation amplitude (a) and growth rate (b) reduce gradually at the early times and increase gradually at the late times with the increasing of initial wavelength. And when the initial perturbation strength (ration of initial amplitude to wavelength) is much smaller, the perturbation growth is mainly dependent on the initial perturbation amplitude and slightly dependent on the initial perturbation wavelength at the late times of

For the effects of three-dimensional initial multi-mode conditions, the case of Richtmyer-Meshkov instability is same as the above Section 4.1. Table 2 lists the initial condition of perturbations, the perturbation strength (PS) is defined as the ratio of initial amplitude to wavelength. Turbulent kinetic energy K, dissipation rate ε, and enstrophy Ω are defined as

> iu00 i〉 2〈ρ〉 þ 〈τii〉

> > i=∂xj〉

i=∂xj〉 〈ρ〉 � 〈τij∂u00

Ω ¼ 1 2 j ω j

eddy simulations, the turbulent kinetic energy and dissipation rate both include two parts named as the resolved-scales (the first term on the right side of Eqs. (46) and (47)) and the

Figure 22 shows the growth history of TMZ width. Figures 23 and 24 show the time histories of the perk values of turbulent kinetic energy and enstrophy, respectively. For the larger initial perturbation strength, the TMZ grows faster, the turbulent kinetic energy is also larger or the turbulence strength is also stronger, the deposited vorticity is larger too. The development of turbulent mixing has a strong dependence on the initial conditions between the initial shock and the impingement of the first reflected rarefaction wave, after that the evolution of the

η0 (mm) 0.07 0.14 0.28 0.56 1.12 PS 0.035 .0.7 0.14 0.28 0.56

i is the velocity fluctuation and 〈 〉 denotes the transverse plane-average. For large-

CASE1 CASE2 CASE3 CASE4 CASE5

2〈ρ〉 ð46Þ

〈ρ〉 ð47Þ

2 ð48Þ

RMI. The same conclusion can be obtained from three-dimensional calculations.

K ¼ 〈ρu00

ε ¼ 〈σij∂u00

subgrid-scales (the second term on the right side of Eqs. (46) and (47)).

turbulent mixing has lost the memory of the initial conditions.

4.3.2. Effects of the initial conditions on the growth of multi-mode RMI and the

82 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

induced turbulent mixing

follows [31]:

where u00

Table 2. Model parameters.

Figure 23. History of turbulent kinetic energy for different models.

Figure 24. History of enstrophy for different models.

#### 4.4. Dynamic characters of multi-mode RM instability and induced turbulent mixing

Figures 25–27 show the spatial profiles of the turbulent kinetic energy, dissipation rate and enstrophy along the motion direction of shock wave at different times before and after reshock individually. They all have a spatial profile similar to Gaussian distribution. The strongest turbulent intensity is located in the center of TMZ. Figures 28–30 show the temporal evolutions of the peak values of the turbulent kinetic energy, dissipation rate, and enstrophy in the spatial profiles, along with their fitted results, respectively. The turbulent kinetic energy, dissipation rate, and enstrophy decay gradually because of dissipation and diffusion. After the initial shock and before the reshock, the turbulent kinetic energy and enstrophy decay with time as a power law ~t θ except the dissipation rate which decays with time as an exponential law ~et=t . One reason is that the TMZ is not fully developed turbulence before reshock and the other reason is that the flow in TMZ is stronger anisotropic between the transverse direction and the axis direction. After reshock and the first reflected rarefaction wave, they all decay with time as the negative exponential law with the closed decay factors at the same

Figure 25. Spatial profiles of the turbulent kinetic energy at different times. (a) Before reshock and (b) after reshock.

Figure 26. Spatial profiles of the turbulent kinetic energy dissipation rate at different times. (a) Before reshock and (b) after reshock.

Figure 27. Spatial profiles of the enstrophy at different times. (a) Before reshock and (b) after reshock.

Figure 28. Time history of the peak value of turbulent kinetic energy.

4.4. Dynamic characters of multi-mode RM instability and induced turbulent mixing

84 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

time as a power law ~t

law ~et=t

after reshock.

Figures 25–27 show the spatial profiles of the turbulent kinetic energy, dissipation rate and enstrophy along the motion direction of shock wave at different times before and after reshock individually. They all have a spatial profile similar to Gaussian distribution. The strongest turbulent intensity is located in the center of TMZ. Figures 28–30 show the temporal evolutions of the peak values of the turbulent kinetic energy, dissipation rate, and enstrophy in the spatial profiles, along with their fitted results, respectively. The turbulent kinetic energy, dissipation rate, and enstrophy decay gradually because of dissipation and diffusion. After the initial shock and before the reshock, the turbulent kinetic energy and enstrophy decay with

θ except the dissipation rate which decays with time as an exponential

. One reason is that the TMZ is not fully developed turbulence before reshock and

the other reason is that the flow in TMZ is stronger anisotropic between the transverse direction and the axis direction. After reshock and the first reflected rarefaction wave, they all decay with time as the negative exponential law with the closed decay factors at the same

Figure 25. Spatial profiles of the turbulent kinetic energy at different times. (a) Before reshock and (b) after reshock.

Figure 26. Spatial profiles of the turbulent kinetic energy dissipation rate at different times. (a) Before reshock and (b)

stage, and be similar to the growth of TMZ width. And then, they all decay asymptotically due to no remarkable energy deposition. Therefore, the turbulent mixing behaves in a statistical self-similar pattern. Figure 31 shows the one-dimensional global energy spectra of three velocity components on a log-log scale at three times. The energy spectra of two transverse components of velocity are too close, and there is a difference between transverse and axis components. The turbulent mixing flow is continuous anisotropic yet the anisotropy weakens gradually. That is to say, the development of the turbulent mixing presents a trend of isotropy.

#### 4.5. Numerical study of the elliptic gas cylinder instability

Our group first performed the experimental and numerical investigations of the elliptic gas cylinder instability. As we know that the initial density distribution of the gas cylinder is hard

Figure 29. Time history of the peak value of turbulent kinetic energy dissipation rate.

Figure 30. Time history of the peak value of enstrophy.

to determine in the experiment, which can only give the one-dimensional radial concentration distribution for circular gas cylinder as an approximate Gaussian function. Drawing to the experience on the one-dimensional diffusive interfacial transition layer with finite thickness for circular gas cylinder, we constructed a two-dimensional diffusive interfacial transition layer with finite thickness for elliptic gas cylinder through numerical simulation [32],

$$\rho(\alpha, \beta) = \chi\_0 \rho\_{SF\_6} e^{-\left[\left(a - a\_{\rm min}\right)^2 + \left(\beta - \beta\_{\rm min}\right)^2\right]/\delta^2} \tag{49}$$

$$\frac{(\mathbf{x} - \mathbf{x}\_0)}{\alpha^2} + \frac{(\mathbf{y} - \mathbf{y}\_0)}{\beta^2} = 1\tag{50}$$

Figure 31. One-dimensional global energy spectra of three components of velocity on a log-log scale.

where x0 ¼ y0 ¼ 0, α ∈ [αmin, αmax], β ∈ [βmin, βmax], αmin ¼ βmin ¼ 1.0 � 10�5 mm, αmax ¼ 6.30 mm, βmax ¼ 2.30 mm, and δ ¼ 6.16 mm. χ0 ¼ 0.71 is the concentration of the elliptic center. The density distribution is plotted in Figure 32. The Mach number of air shock wave is 1.25. The simulations (see Figure 33) reproduce the elliptic gas cylinder instability experiment very well, they achieve to a good agreement qualitatively and quantitatively, some salient features of the vortex pairs are obtained clearly.

Figure 34 shows the vorticity at the center of the core and the distance between the two vortex cores of these two simulations. When the incident shock accelerates the elliptic gas cylinder along the major axis, the absolute vorticity in the vortex core |ωcore| is larger; but the distance between two vortex cores is larger when the incident shock accelerates the cylinder along the minor axis. So, the effect of convergence is stronger when the incident shock accelerates the elliptic gas cylinder along the major axis, for which the gas jet appears at the symmetry center.

#### 4.6. Foundation of new RM instability and mechanism

to determine in the experiment, which can only give the one-dimensional radial concentration distribution for circular gas cylinder as an approximate Gaussian function. Drawing to the experience on the one-dimensional diffusive interfacial transition layer with finite thickness for circular gas cylinder, we constructed a two-dimensional diffusive interfacial transition layer

> e �½ðα�αminÞ 2 þðβ�βminÞ 2 �=δ2

α2 þ ðy � y0Þ

ð49Þ

β2 ¼ 1 ð50Þ

with finite thickness for elliptic gas cylinder through numerical simulation [32],

ðx � x0Þ

ρðα, βÞ ¼ χ0ρSF6

Figure 30. Time history of the peak value of enstrophy.

Figure 29. Time history of the peak value of turbulent kinetic energy dissipation rate.

86 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

At present, the investigation for the Richtmyer-Meshkov instability is performed in a uniform flow field. We first study the Richtmyer-Meshkov instability and turbulent mixing in a

Figure 32. Density images and distributions of the (a) circular and (b) elliptic SF6 gas cylinder initially constructed.

Figure 33. Experimental evolution images and numerical simulation results by MVFT at t ¼ 200, 300, 400, 500, and 600 µs, the experiments corresponding (a), (c), and (e), and the simulations corresponding (b), (d), and (f).

nonuniform flow field by experiment and numerical simulations and find some new phenomena. Figure 35 shows the test section schematic of shock tube. The air shock wave with Mach number 1.27 accelerates the two-mode sinusoidal air/SF6 interface (amplitude A01 ¼ 5 m, A02 ¼ 7.5 mm). The experimental Schlieren images (gray images in Figure 38) show the transmitted shock wave and the interface all incline which is different from the familiar RMI. We think this

Figure 34. Vorticity at the center of the cores and the distance between the two vortex cores from simulations.

Figure 35. Initial structure diagram in the shock tube.

nonuniform flow field by experiment and numerical simulations and find some new phenomena. Figure 35 shows the test section schematic of shock tube. The air shock wave with Mach number 1.27 accelerates the two-mode sinusoidal air/SF6 interface (amplitude A01 ¼ 5 m, A02 ¼ 7.5 mm). The experimental Schlieren images (gray images in Figure 38) show the transmitted shock wave and the interface all incline which is different from the familiar RMI. We think this

Figure 33. Experimental evolution images and numerical simulation results by MVFT at t ¼ 200, 300, 400, 500, and 600 µs,

the experiments corresponding (a), (c), and (e), and the simulations corresponding (b), (d), and (f).

Figure 32. Density images and distributions of the (a) circular and (b) elliptic SF6 gas cylinder initially constructed.

88 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 36. Density profiles of nonuniform flows with Gaussian function and uniform flows in a vertical direction.

may be owing to the nonuniformity of flow field. Also, drawing to the experience on the diffusive interfacial transition layer of circular gas cylinder, we construct a nonuniform flow field with a Gaussian distribution of density along the direction perpendicular to the shock motion direction (as shown in Figure 36). The numerical result (shown in Figure 37) confirms this idea [12, 33, 34].

Figure 37. Simulated density distribution in a nonuniform flow field by using Gaussian function.

Figure 38. Experimental schlieren photography images and numerical simulation results by MVFT2D at a certain time series (the sizes of the pictures are ones of the test window [0.038, 0.25 m] [0.0, 0.2 m]).

Figure 38 shows the evolution of the interface and the propagation of the transmitted shock wave, the calculated results agree well with the experiment. Due to the nonuniform flow field of SF6 gas, the propagating velocity of transmitted shock wave in the upper part of shock tube is faster than in the bottom of shock tube, and it forms an oblique shock wave front and the interface. Figure 39 shows the shock front locations between experiment and numerical simulations, Figure 40 shows the perturbation amplitude of experiment, numerical simulation, and theories, they are in good agreement. Figure 41 shows the calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields. As we can see that, at late times, the growth of small perturbation in a low-density zone catches up and exceeds the large perturbation in a high-density zone, which is opposite to the case of uniform flow field. Figure 42 shows the perturbation amplitudes of four different initial amplitudes 2.5, 5.0, 7.5, and 10 mm in low- and high-density zones of nonuniform flow fields. It shows that the effect of initial amplitude on the growth of RM instability in the nonuniform flow field is different from the case of uniform one, which can be explained by the baroclinic vorticity shown in Figures 43 and 44, the baroclinic vorticity produced in the low-density zone is larger than that in the high-density zone.

Figure 39. Shock front locations of the experiment and calculated results on the three test lines.

Figure 40. Perturbation amplitudes of the experiment, simulation, and comparison with theories.

Figure 38 shows the evolution of the interface and the propagation of the transmitted shock wave, the calculated results agree well with the experiment. Due to the nonuniform flow field of SF6 gas, the propagating velocity of transmitted shock wave in the upper part of shock tube is faster than in the bottom of shock tube, and it forms an oblique shock wave front and the interface. Figure 39 shows the shock front locations between experiment and numerical simulations, Figure 40 shows the perturbation amplitude of experiment, numerical simulation, and theories, they are in good agreement. Figure 41 shows the calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields. As we can see that, at late times, the growth of small perturbation in a low-density zone catches up and exceeds the large perturbation in a high-density zone, which is opposite to the case of uniform flow field. Figure 42 shows the perturbation amplitudes of four different initial amplitudes 2.5, 5.0, 7.5, and 10 mm in low- and high-density zones of nonuniform flow fields. It shows that the effect of initial amplitude on the growth of RM instability in the nonuniform flow field is different from the case of uniform one, which can be explained by the baroclinic vorticity shown in Figures 43 and 44, the baroclinic vorticity produced in the low-density zone is larger than

Figure 38. Experimental schlieren photography images and numerical simulation results by MVFT2D at a certain time

Figure 37. Simulated density distribution in a nonuniform flow field by using Gaussian function.

90 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

series (the sizes of the pictures are ones of the test window [0.038, 0.25 m] [0.0, 0.2 m]).

that in the high-density zone.

Figure 41. Calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields.

Figure 42. Perturbation amplitudes of four different initial amplitude 2.5, 5.0, 7.5 and 10 mm in the low and high density zones of nonuniform flow field.

Figure 43. Average vorticity when initial amplitude is 5.0 mm which is in the low and high density zones of nonuniform flow field at 1 ms.

Figure 44. Baroclinic vorticity in the uniform and nonuniform flow fields with the initial amplitude group (A01 ¼ 2.5 mm, A02 ¼ 7.5 mm) at 1.0 ms.

Figure 42. Perturbation amplitudes of four different initial amplitude 2.5, 5.0, 7.5 and 10 mm in the low and high density

92 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 43. Average vorticity when initial amplitude is 5.0 mm which is in the low and high density zones of nonuniform

Figure 44. Baroclinic vorticity in the uniform and nonuniform flow fields with the initial amplitude group (A01 ¼ 2.5 mm,

zones of nonuniform flow field.

A02 ¼ 7.5 mm) at 1.0 ms.

flow field at 1 ms.

Figure 45. Density contour of the numerical simulated by MVFT at a time series. The left column with a uniform initial condition, the middle column with a δ1 nonuniform Gaussian function, and the right column with a δ2 nonuniform Gaussian function. The small arrow denotes the direction of propagation of the shock wave fronts before reshock the interface.

Figure 46. TMZ width of RM instability in the initial uniform and nonuniform flow fields.

For the above investigations of RM instability in the nonuniform flow field, they are all without reshock. In the following section, we will study the effect of nonuniformity of flow field and the reshock on the RM instability. Here, we consider two kinds of nonuniformity: the nonuniform coefficient δ1 ¼ 0.6162 m and δ2 ¼ 0.4961 m. Figure 45 shows the density contour of the numerical simulated by MVFT at a time series, the left column with a uniform initial conditions, the middle column with a δ1 nonuniform Gaussian function, and the right column with a δ2 nonuniform Gaussian function. There is a significant difference between the uniform and nonuniform flows before reshock, but the difference decreases in evidence after reshock. Figure 46 shows the TMZ width of RM instability in the initial uniform and nonuniform flow fields. It points out that the growth of the TMZ width for the initial nonuniform flow field is greater than that for the uniform flow field, and the lesser the nonuniform coefficient, the higher the growth rate of TMZ width, but the difference for the three different flow configurations diminishes after reshock.

## 5. Prospect

Now, we investigated the Richtmyer-Meshkov instability and the induced turbulent mixing in fluid flow by using large-eddy simulation, but the turbulent mixing is a complicated threedimensional problem with multiple time-space scales, and the more engineering applications involve the materials mixing with strength. So, the future work will be carried out in the following aspects: (a) the direct numerical simulations with high precision and high resolution on the platform of supercomputer, (b) the interface instability and turbulent materials mixing with strength, this may involve the material properties such as the deformation, fracture, melt, phase transition, material microstructures, etc., and it may make this problem to be more difficult.

## Acknowledgements

Professor Jingsong Bai and Associate professor Tao Wang from the Institute of Fluid Physics, China Academy of Engineering Physics, People's Republic of China and some other from the same institute should be thanked for the contributions and helps to these works, such as Prof. Ping Li, Dr. Liyong Zou, Kun Liu, Mr. Jinhong Liu, Min Zhong and Bing Wang, and Jiaxin Xiao.

The authors would like to thank the support from "Science Challenge Project" (no. TZ2016001), the National Natural Science Foundation of China (nos. 11372294, 11202195, 11072228, and 11532012), the Science Foundation of China Academy of Engineering Physics (nos. 2011B020 2005, 2011A0201002, and 2008B0202011), and the Foundation of National Key Laboratory of Shock Wave and Detonation Physics (no. 9140C670301150C67290).

## Author details

For the above investigations of RM instability in the nonuniform flow field, they are all without reshock. In the following section, we will study the effect of nonuniformity of flow field and the reshock on the RM instability. Here, we consider two kinds of nonuniformity: the nonuniform coefficient δ1 ¼ 0.6162 m and δ2 ¼ 0.4961 m. Figure 45 shows the density contour of the numerical simulated by MVFT at a time series, the left column with a uniform initial conditions, the middle column with a δ1 nonuniform Gaussian function, and the right column with a δ2 nonuniform Gaussian function. There is a significant difference between the uniform and nonuniform flows before reshock, but the difference decreases in evidence after reshock. Figure 46 shows the TMZ width of RM instability in the initial uniform and nonuniform flow fields. It points out that the growth of the TMZ width for the initial nonuniform flow field is greater than that for the uniform flow field, and the lesser the nonuniform coefficient, the higher the growth rate of TMZ width, but the difference for the three different flow configura-

94 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Now, we investigated the Richtmyer-Meshkov instability and the induced turbulent mixing in fluid flow by using large-eddy simulation, but the turbulent mixing is a complicated threedimensional problem with multiple time-space scales, and the more engineering applications involve the materials mixing with strength. So, the future work will be carried out in the following aspects: (a) the direct numerical simulations with high precision and high resolution on the platform of supercomputer, (b) the interface instability and turbulent materials mixing with strength, this may involve the material properties such as the deformation, fracture, melt, phase transition, material microstructures, etc., and it may make this problem to be more

Professor Jingsong Bai and Associate professor Tao Wang from the Institute of Fluid Physics, China Academy of Engineering Physics, People's Republic of China and some other from the same institute should be thanked for the contributions and helps to these works, such as Prof. Ping Li, Dr. Liyong Zou, Kun Liu, Mr. Jinhong Liu, Min Zhong and Bing Wang, and Jiaxin Xiao. The authors would like to thank the support from "Science Challenge Project" (no. TZ2016001), the National Natural Science Foundation of China (nos. 11372294, 11202195, 11072228, and 11532012), the Science Foundation of China Academy of Engineering Physics (nos. 2011B020 2005, 2011A0201002, and 2008B0202011), and the Foundation of National Key Laboratory of

Shock Wave and Detonation Physics (no. 9140C670301150C67290).

tions diminishes after reshock.

5. Prospect

difficult.

Acknowledgements

Jingsong Bai1,2 and Tao Wang1 \*

\*Address all correspondence to<EMAIL_ADDRESS>
1 Institute of Fluids Physics, China Academy of Engineering Physics, Mianyang, Sichuan, People's Republic of China

2 National Key Laboratory of Shock Wave and Detonation Physics (LSD), Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang, Sichuan, People's Republic of China

## References

[24] Bates KR, Nikiforakis N, Holder D. Richtmyer-Meshkov instability induced by the interaction of a shock wave with a rectangular block of SF6. Physics of Fluids. 2007;19:036101

[10] Bai JS, Wang T, Li P, Zou LY, Liu CL. Numerical simulation of the hydrodynamic instability experiments and flow mixing. Science in China Series G: Physics, Mechanics &

[11] Wang T, Bai JS, Li P, Liu K. Large-eddy simulations of the Richtmyer-Meshkov instability of rectangular interface accelerated by shock waves. Science in China Series G: Physics,

[12] Bai JS, Liu JH, Wang T, Zou LY, Li P, Tan DW. Investigation of the Richtmyer-Meshkov instability with double perturbation interface in nonuniform flows. Physical Review E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics. 2010;81(2):

[13] Smagorinsky J. General circulation experiments with the primitive equations. Monthly

[14] Vreman AW. An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic

[15] Moin P, Squires K, Cabot W, Lee S. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A (Fluid Dynamics). 1991;3(11):

[16] Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multiscale modelling of a Richtmyer-Meshkov instability with reshock. Journal of Fluid Mechanics. 2006;557:29–61

[17] Bai JS. High resolution numerical methods and adaptive mesh refinement algorithm for compressible multi-fluid dynamics [PhD thesis]. China Academy of Engineering Physics;

[18] Wang T, Bai JS, Li P, et al. Two and three dimensional numerical investigations of the single-mode Richtmyer-Meshkov instability. Chinese Journal of High Pressure Physics.

[19] Wang T, Bai JS, Li P, Zhong M. The numerical study of shock-induced hydrodynamic

[20] Malamud G, Elabz Y, Leinov E, et al. Two- and three-dimensional numerical simulations of the interaction of a Richtmyer-Meshkov instability induced turbulent mixing zone

[21] Bai JS, Wang T, Liu K, et al. Large-Eddy simulation of the Three-dimensional experiment on Richtmyer-Meshkov instability induced turbulence. International Journal of Astron-

[22] Poggi F, Thorembey MH, Rodriguez G. Velocity measurements in turbulent gaseous mixtures induced by Richtmyer–Meshkov instability. Physics of Fluids. 1998;10:2698–2712 [23] Wang T, Li P, Bai JS, Tao G, Wang B, Zou LY. Large-eddy simulation of the Richtmyer-

Meshkov instability. Canadian Journal of Physics. 2015;93(10):1124–1130

with several re-shocks. In: 10th IWPCTM; 17–21 July 2006, Paris, France. 2006.

theory and applications. Physics of Fluids. 2004;16(10):3670–3681

instability and mixing. Chinese Physics B. 2009;18(3):1127–1135

Astronomy. 2009;52(12):2027–2040

Weather Review. 1963;91:99–164

056302

2746–2757

2003. (in Chinese)

2013;27(2):277–286

omy and Astrophysics. 2012;2:28–36

Mechanics & Astronomy. 2010;53(5):905–914

96 Turbulence Modelling Approaches - Current State, Development Prospects, Applications