Lucas Tallois
email: lucas.tallois@cea.fr

Simon Peluchon
email: simon.peluchon@cea.fr

Philippe Villedieu
email: philippe.villedieu@onera.fr

A second-order extension of a robust implicit-explicit acoustic-transport splitting scheme for two-phase flows

Keywords: Two-phase flows, MUSCL, Unstructured meshes, Implicit, Compressive limiter

Diffuse interface methods have proven their ability to simulate complex two-phase flows. A number of robust numerical schemes have been developed to simulate such flows involving large density and pressure ratios. Diffusion induced by these methods, however, makes it difficult to localize the interface between the two fluids. To overcome this issue, while retaining the advantages of diffuse interface methods, a second-order extension using a Monotonic Upstream-centered Scheme for Conservation Laws-type (MUSCL-type) method of the implicit-explicit acoustic-transport splitting scheme introduced in [40] is presented. A specific compressive limiter is used for the volume fraction in order to limit the diffusion of the interface between the two fluids.

Numerical simulations are presented to illustrate the capability of the proposed new method to simulate highly complex compressible two-phase flows.

Introduction

This work takes place in the context of liquid ablation. When an object enters the atmosphere at very high speed (see Fig. 1), it is subjected to significant heat flux which may degrade its structure. Depending on its composition, the object may sublimate or liquefy. Thus, when a liquid phase appears, it is necessary to deal with the flow of both the liquid and gas phases by an appropriate model. The modelling of two-phase flows has been the subject of numerous studies. Several methods have been devel-oped, each with advantages and disadvantages. They can be grouped into two distinct families: the sharp interface methods and the diffuse interface methods.

Among these methods, the first one is the so-called Lagrangian method [START_REF] Vonneumann | A method for the numerical calculation of hydrodynamic shocks[END_REF]. Each phase can be modeled with different equations, and the interface is followed explicitly by moving the mesh at the material speed. In the so-called ALE methods (Arbitrary Lagrangian-Eulerian methods), the mesh follows the displacement of the interface. The mesh can, however, be subject to strong distortions that can impact the robustness of the calculations.

Moreover, among the Eulerian methods with no or almost no diffusion, there are the Front capturing methods as Volume Of Fluid (VOF) [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and Moment Of Fluid (MOF) [START_REF] Ahn | Multi-material interface reconstruction on generalized polyhedral meshes[END_REF] methods, the Front tracking methods [START_REF] Caramana | Elimination of Artificial Grid Distortion and Hourglass-Type Motions by Means of Lagrangian Subzonal Masses and Pressures[END_REF] and the Level Set method [START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. Although widely used for their satisfactory results, these methods are complex to implement and some may be not conservative, which can be a serious drawback depending on the targeted applications.

Diffuse interface methods, on the other hand, are based on two-phase models, i.e. models containing conservation equations for each of the two phases that are artificially assumed to be present at any point in space. The equations are solved on a fixed mesh in a Eulerian manner. The same equations are solved over the whole domain. These methods allow the diffusion of the interface between each fluid. Thus, mixing zones appear, corresponding to a numerical spreading of the interface over a few cells. These areas require special processing to preserve the thermodynamic consistency of the model. Many studies have been carried out with the aim of developing two-phase models, following the work of Baer & Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] and their seven-equation model. Indeed, there is a whole hierarchy of two-phase models deduced from the seven-equation model [START_REF] Flåtten | Relaxation two-phase flow models and the subcharacteristic condition[END_REF][START_REF] Lund | A hierarchy of relaxation models for twophase flow[END_REF]. A local velocity equilibrium between the two fluids [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] allows us to reduce to a six-equation model. If the equilibrium of the pressures is assumed, the five-equation model of Kapila et al. [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Kreeft | A new formulation of Kapila's five-equation model for compressible two-fluid flow, and its numerical treatment[END_REF] or the reduced five-equation model of [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF]] can be obtained.

Another approach to simulate two-phase flow consists in writing the same set of equations for each fluid, and to close the model by defining the mixing quantities as the quantities of each fluid weighted by its volume fraction [START_REF] Chanteperdrix | A Compressible Model for Separated Two-Phase Flows Computations Volume 1[END_REF][START_REF] Grenier | An accurate low-mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF][START_REF] Blanchard | Modélisation et simulation multi-échelles de l'atomisation d'une nappe liquide cisaillée[END_REF]]. An isobaric closure assumption and the resolution of a 2 × 2 system allows the volume fraction in the mixing region to be determined.

Another approach existing in the class of diffuse interface methods is the phase field method [START_REF] Anderson | Diffuse-Interface Methods in Fluid Mechanics[END_REF][START_REF] Heider | A review on phase-field modeling of hydraulic fracturing[END_REF]. This method introduces an auxiliary field (the phase field or color function) that acts as an order parameter taking two distinct values in each of the phases, with a smooth change between the two values in the area around the interface. In two-phase flow models, the notion of color function or phase field parameter is not explicitly introduced but its role is implicitly played by the volume fraction of one of the two fluids.

In this work, a diffuse interface method based on a two-phase flow model is considered. The advantages of such method are numerous: the same equations are solved over the whole domain, interfaces do not require any specific processing and appearance of new interfaces as well as the changes of topology are achieved naturally. In the case where one wishes to model the two-phase flow between a gaseous phase and a liquid phase, where the second phase is not present initially and appears from solid fusion, these methods seem the most relevant and the easiest to implement. The counterpart of these methods is the diffusion of the interface between the two fluids, which can deteriorate the solution when the fluids are immiscible. Each of the fluids is assumed to be present at any point in space, with a volume fraction very close to 0 or 1 depending on whether we are in fluid 1 or fluid 2. It is only in the vicinity of the interface between the fluids that the volume fraction of each of the fluids pass through the value 1/2. Since we do not consider that fluids slide relative to each other at the mesh scale, a two-phase model with one velocity and one pressure is used. This avoids having to consider a model with two velocities and two pressures [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] which poses greater theoretical and numerical difficulties and is not necessary when the mixing zone is purely numerical and as fine as possible, and does not model a zone with inclusions such as drops or bubbles of a size smaller than the scale resolved by the mesh. As the two phases are not initially mixed, the non-conservative term in the Kapila et al. model [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Kreeft | A new formulation of Kapila's five-equation model for compressible two-fluid flow, and its numerical treatment[END_REF] is null. The reduced five-equation model of [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF] will therefore be used in this work. Several methods have been developed in order to limit this numerical diffusion of the interface. Kokh and Lagoutière used the anti-diffusive scheme [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF] developed in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF], which consists in using a flux as downwind as possible, while guaranteeing the stability and consistency of the scheme. This method limits the diffusion of the interface to 2 mixture cells. Jung et al. [START_REF] Helluy | OpenCL numerical simulations of two-fluid compressible flows with a 2D random choice method[END_REF] combined the Glimm's random projection method near the material interface with the upwind scheme to simulate two-phase flow without any diffusion cells. They extended this method to the second order in space and time. Extension to higher dimension is achieved with directional splitting. Shyue and Xiao [START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach[END_REF] adapted the THINC method of [START_REF] Xiao | A simple algebraic interface capturing scheme using hyperbolic tangent function[END_REF] as a sharpening technique to simulate two-phase flows. The main idea is to reconstruct the solution with the hyperbolic tangent function, controlling the compressive character of the reconstruction with a βparameter. This is widely combined with the Boundary Variation Diminishing algorithm [START_REF] Sun | Boundary variation diminishing (BVD) reconstruction: A new approach to improve Godunov schemes[END_REF][START_REF] Deng | High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces[END_REF][START_REF] Cheng | Lowdissipation BVD schemes for single and multi-phase compressible flows on unstructured grids[END_REF] to simulate two-phase flows on structured and unstructured grids. In the context of MUSCL-type schemes [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], where flux are computed with piecewise linear reconstruction with slope limitation, slope limiters can be designed to have compressive properties and give a sharp interface when applied to the volume fraction transport equation. This approach has been widely used in the literature as for instance in Qian et al. [START_REF] Qian | A freesurface capturing method for two fluid flows with moving bodies[END_REF], Blanchard [START_REF] Blanchard | Modélisation et simulation multi-échelles de l'atomisation d'une nappe liquide cisaillée[END_REF], Chiapolino et al. [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF] and De Vuyst et al. [START_REF] De Vuyst | A geometrically accurate low-diffusive conservative interface capturing method suitable for multimaterial flows[END_REF].

This paper is an extension of the work presented in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF] where the splitting strategy of [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF] is used to solve the fiveequation system of [START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF][START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF]. The strategy consists in solving the acoustic part with an implicit scheme while the transport step is solved explicitly. The effects related to heat dissipation and viscosity can be easily added to the first step [START_REF] Peluchon | Development of numerical methods to simulate the melting of a thermal protection system[END_REF], which justifies its implicit treatment. We will focus here only on the hyperbolic part of the equations. The aim is to increase the accuracy of the first-order method described in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF] and to improve the resolution of the material interface. A MUSCL-type scheme on an unstructured grid is used, in order to improve the accuracy of the numerical scheme. A strategy to improve the resolution of the two-phase interface on an unstructured grid is employed. Indeed, the method must be able to handle a body-fitted mesh composed of quadrilaterals and deformed by the ablation. In the context of the splitting strategy, a new and more accurate implicit-explicit time-scheme is presented.

The paper is organized as follows. First, the governing equations used to simulate two-phase flow are introduced and the first order numerical scheme is briefly recalled in Section 2. Then, a second-order in space extension and improvement of the time scheme are derived in Section 3. Some numerical results are presented in Section 4.

Model and numerical scheme

The five-equation system

We denote by ρ k , ε k and p k the density, internal energy and pressure of fluid k = 1, 2. Each fluid is equipped with an Equation Of State (EOS) of the form p k = p k (ρ k , ε k ). The volume fraction z k , such that z 1 + z 2 = 1, allows the position of the interface to be located. In the sequel, we denote by z = z 1 the volume fraction of the first fluid. The mixture density and internal energy are given by

ρ = ρ 1 z + ρ 2 (1 -z), ρε = ρ 1 ε 1 z + ρ 2 ε 2 (1 -z).
(2.1)

Both fluids share the same velocity u and the same pressure p. The five-equation system of [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF] for non-miscible fluids with isobaric closure reads

                   ∂ t (ρ 1 z) + ∇ • (ρ 1 zu) = 0, ∂ t (ρ 2 (1 -z)) + ∇ • (ρ 2 (1 -z)u) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + ∇ • (ρEu + pu) = 0, ∂ t z + u • ∇z = 0, (2.2)
where

E = ε + ||u|| 2 2
is the total energy of the two-fluid media.

Each fluid is considered to be governed by a stiffened gas EOS, since we consider flows with a liquid phase and a gas phase. The EOS of phase k thus reads

p k = ρ k ε k (γ k -1) -γ k π k , (2.3)
where γ k > 1 is the adiabatic exponent and π k ≥ 0 is a reference pressure. Note that perfect gas EOS is obtained with π k = 0. The mixture pressure p is the solution of the system

p 1 (ρ 1 , ε 1 ) = p 2 (ρ 2 , ε 2 ), ρε = ρ 1 ε 1 z + ρ 2 ε 2 (1 -z).
(2.4)

In the case of two stiffened gases, (2.4) can be solved explicitly and we get the following expression for the mixture pressure p = ρε(γ -1) -γπ,

where the adiabatic exponent mixture γ and the reference pressure π are given by

γ = 1 + 1 2 k=1 z k γ k -1 and π = γ -1 γ 2 k=1 z k γ k π k γ k -1 . (2.6)
The speed of sound c k of each phase is defined by

c 2 k = ∂p k ∂ρ k | s k
, where s k is the entropy. The mixture sound speed for two stiffened gas reads

c 2 = γ p + π ρ . (2.7)
Note before going further that the system (2.2) can be written, introducing the mass fraction of the first fluid y = ρ1z ρ , in the following form

                   ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρy) + ∇ • (ρyu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + ∇ • (ρEu + pu) = 0, ∂ t z + u • ∇z = 0.
(2.8)

The evolution equation of ρ is obtained by summing the evolution equation of ρ 1 z 1 and ρ 2 z 2 of system (2.2). The evolution equation of ρy is nothing but the evolution equation of ρ 1 z 1 by definition. This form allows us to write in the Lagrangian form the acoustic step, described below.

Remark 1

If the volume fraction only takes the values 0 and 1, then (2.8) is equivalent to the 4-equation model i.e (2.8) without the equation on ρy.

Splitting strategy

We will briefly recall the method of solving the fiveequation system, where a detailed approach is described in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. The expansion of system (2.8) reads

                   ∂ t ρ + ρ∇ • u + u • ∇ρ = 0, ∂ t (ρy) + ρy∇ • u + u • ∇(ρy) = 0, ∂ t (ρu) + ρu∇ • u + ∇p + u • ∇(ρu) = 0, ∂ t (ρE) + ρE∇ • u + ∇ • (pu) + u • ∇(ρE) = 0, ∂ t z + u • ∇z = 0.
(2.9) The five-equation system is then split into two subsystems [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. The first system contains only the acoustic waves, while the second system takes into account the propagation of material waves through the fluid.

The system corresponding to the acoustic step is hyperbolic with eigenvalues 0 and ±c. It is rewritten and then solved in its Lagrangian form, which reads

∂ t V + ϑ∇ • G(V ) = 0, (2.10)
where V = (ϑ, y, u, E, z) are the Lagrangian variables, G(V ) = (u, 0, p, pu, 0) is the flux and ϑ = 1/ρ is the specific volume. Note that during this step, volume fraction and mass fraction are not modified. The transport step system, where the propagation speed is only u, reads

∂ t U + u∇ • U = 0, (2.11)
with U = (ρ, ρy, ρu, ρE, z) the conservative variables.

Numerical scheme

Let us first consider a domain Ω ∈ R 2 discretized in N cells Ω i such that v(i) is the set of cells neighboring the cell Ω i by the edges, |Ω i | is the area of the cell Ω i , |Γ ij | is the length of the edge common to the cells Ω i and Ω j and n ij is the unit vector normal to Ω i (see Fig. 2). The mass centers x i of cell Ω i is defined as

x i = 1 |Ω i | Ωi xdx.
(2.12) The time is discretized by t n = n∆t for n ∈ N, where ∆t > 0 is the time step. In the finite volume approach that is used, the notation φ n i is an approximation of

Ω i Ω j • x i • x j • x f Γ ij n ij
1 |Ω i | Ωi φ(x, t n )dx,
for any quantity φ(x, t).

In the sequel, the first order in space and time scheme is described.

Acoustic step

An approximate four-state Riemann solver [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF], or the Lagrangian scheme EUCCLHYD [START_REF] Maire | A cellcentered lagrangian scheme for two-dimensional compressible flow problems[END_REF] can be used to solve this system. The numerical scheme for the acoustic step reads

V † = V n - ∆t |Ω i | ϑ n i j∈v(i) |Γ ij | G # ij , (2.13)
where

G # ij = (-ū ij , 0, pij n ij , pij ūij , 0)
# is the numerical flux and † corresponds to the intermediate state between the acoustic step and the transport step. The scheme can be explicit by taking # = n or implicit with # = †. For the approximate four-state Riemann solver of Gallice [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF], the flux reads

ūij = C- ij u i + C+ ij u j C- ij + C+ ij • n ij - p j -p i C- ij + C+ ij , pij = C+ ij p i + C- ij p j C- ij + C+ ij -C- ij C+ ij u j -u i C- ij + C+ ij • n ij (2.14)
where Cij and C+ ij are the Riemann solver slopes. The choice of slopes is crucial to guarantee the positivity of the intermediate states of the Riemann solver and thus of the associated Godunov type scheme [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF] (see also [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF]). Their computation is detailed in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. In the explicit case, the Courant-Friedrichs-Lewy (CFL) condition ensures that no wave cross within the same cell. It is given by

∆t max 1≤i≤N ϑ i |Ω i | max j∈v(i) |Γ ij | Cij < 1 2 . (2.15)
When the two phases have large pressure or density ratios, the CFL (2.15) related to the acoustic system can be very constraining. This is the case for liquid/gas interactions. Thus, the acoustic step is solved with an implicit time-scheme [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. This approach is extensively detailed in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. We will recall here the main ideas. For the sake of simplicity, a one-dimensional problem is considered. The extension in dimension 2 follows the same lines. The speed of sound and the Riemann solver slopes are frozen at time n. Since the numerical flux G ij depends only on the velocity and the pressure, we solve the following subsystem

   ∂ t p + a 2 ϑ∂ x u = 0, ∂ t u + ϑ∂ x p = 0, (2.16)
where a = ρc is the Lagrangian sound speed. We now look for the solution of the system

X † -X n + ∆t ∆x ϑ n [[H]] = 0, (2.17)
where X = (p i , u i ) and

[[H]] = (a 2 i (ū † i+1/2 -ū † i-1/2 ), p † i+1/2 -p † i-1/2 ), (2.18)
are the jumps of the numerical flux, computed with (2.14).

Since the system (2.16) is linear in (p, u), it can be rewritten as

X † -X n + ∆t ∆x M X † = 0, (2.19)
where

M X = ϑ n [[H]]
. The matrix M depends on the choice made to evaluate the flux. It is given in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF] for the four-state Riemann solver of [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF]. Basically, it depends on the Riemann solver slope and the sound speed. In practice, the system obtained is solved under the delta-form

I d + ∆t ∆x M X † -X n = - ∆t ∆x M X n . (2.20)
The resolution of this linear system then gives a new velocity and a prediction of the pressure. These pressure predictions are used to evaluate the numerical flux G † ij and update the acoustic variables ϑ † and E † explicitly with formulae (2.13). The real pressure is then determined by the EOS (2.5). It does not seem possible to prove that the pressure obtained at the end of the acoustic step is close to the pressure given by the EOS.

Transport step

The transport system is resolved under the equivalent form

∂ t ψ + ∇ • (uψ) -ψ∇ • u = 0, (2.21)
where ψ corresponds to the variables (ρ, ρy, ρu, ρE, z) of the system (2.11). The transport equation (2.21) consists of a first conservative term ∇ • (uψ) followed by a second non-conservative term ψ∇ • u. The first-order numerical scheme of the transport step is as follows

ψ n+1 i = ψ † i - ∆t |Ω i | j∈v(i) |Γ ij | ψ † ij u ij + ψ † i ∆t |Ω i | j∈v(i) |Γ ij | u ij .
(2.22) It is solved with an explicit time scheme. The most natural choice is to use the upwind scheme to solve the conservative term. A judicious choice on the material velocity u ij gives a globally conservative scheme. This choice consists in taking the opposite of the first component of the acoustic flux [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. In the case of the Gallice's Riemann solver [START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF], we have

u ij = ū# ij . The transport step is stable under the following CFL condition ∆t max 1≤i≤N   1 |Ω i | j∈v(i) |Γ ij | |u ij |   < 1.
(2.23)

Global scheme

The overall algorithm for a time step between t n and t n+1 reads

Step 1: From a state (ρ, ρy, ρu, ρε, z) n , compute (ρ, ρy, ρu, ρε, z) † , the approximation of the acoustic system (2.13).

Step 2: Find the fluid state (ρ, ρy, ρu, ρε, z) n+1 by solving the transport system (2.11) with the initial state (ρ, ρy, ρu, ρε, z) † .

The global scheme reads

                                           ρ n+1 i = ρ n i - ∆t |Ω i | j∈v(i) |Γ ij | ρ † ij ū# ij , (ρy) n+1 i
= (ρy) n i - ∆t |Ω i | j∈v(i) |Γ ij | (ρy) † ij ū# ij , (ρu) n+1 i
= (ρu) n i - ∆t |Ω i | j∈v(i) |Γ ij | (ρu) † ij ū# ij + p# ij , (ρE) n+1 i = (ρE) n i - ∆t |Ω i | j∈v(i) |Γ ij | (ρE) † ij ū# ij + p# ij ū# ij , z n+1 i = z n i - ∆t |Ω i | j∈v(i) |Γ ij | z n ij ū# ij + z n i ∆t |Ω i | j∈v(i) |Γ ij | ū# ij . (2.24)
This scheme is clearly conservative for the densities, momentum, and total energy equations. It preserves the contact discontinuities, namely constant u and p states.

Remark 2

Numerical flux depends on ∆t. Indeed, the variables of the flux with the exponent † are updated by the acoustic step (2.13) which depends on the time step. Thus, stationary solutions depend on the time step. We can see the two steps as an implicit method to calculate the numerical flux associated with (2.24).

Second-order extension

Second-order in space extension

When using diffuse interface methods, the volume fraction will naturally be diffused at the interface between the two fluids considered. The simple first-order in space scheme used in both the acoustic and the transport step will produce a numerical diffusion of the solution that is far too large to accurately follow the interface. To overcome this problem, we use a second-order MUSCL-type method [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF] called U-MUSCL [START_REF] Burg | Higher Order Variable Extrapolation for Unstructured Finite Volume RANS Flow Solvers[END_REF]. It consists in using in the numerical flux the polynomial reconstructions αij and αji of the solution on Γ ij such that

αij = α i + κ 2 (α j -α i ) + (1 -κ)∇α i .(x f -x i ), αji = α j + κ 2 (α i -α j ) + (1 -κ)∇α j .(x f -x j ), (3.1)
with κ ∈ [-1, 1] a parameter. The use of the mass center x i given by (2.12), instead of any point in cell Ω i , ensures that the reconstructions αn i of variable α n i are conservative, i.e.

1 |Ω i | Ωi αn i (x)dx = α n i . (3.2)
By writing (3.1) it in the form

αij = κ (α j + α i ) 2 + (1 -κ) (α i + ∇α i .(x f -x i )) , (3.3)
it can be seen as a combination of the linear interpolation (first term) and the linear extrapolation (second term) of αij . It is an extension of van Leer's κ-scheme on unstructured grids. Note that κ = 0 returns the conventional MUSCL reconstruction. As shown in [START_REF] Nishikawa | On the loss and recovery of second-order accuracy with U-MUSCL[END_REF], the U-MUSCL scheme with κ = 1/3 gives more accurate results than κ = 0. However, this method is second-order only if the center of the face x f lies exactly between the mass centers of cells x i and x j (see Fig. 2). Therefore, a modification from [START_REF] Nishikawa | On the loss and recovery of second-order accuracy with U-MUSCL[END_REF] must be employed to preserve the second-order accuracy. This consists in modifying the reconstruction (3.1) such that

αij = α i + κ 2 (α p -α i ) + (1 -κ)∇α i .(x f -x i ), (3.4)
where

α p = α j + ∇α j .(2x f -x i -x j ). (3.5)
The U-MUSCL scheme (3.4) with the second-order preserving modification will be used as the reconstruction method due to its simplicity and accuracy when using κ = 1/3. As shown in [START_REF] Padway | Resolving Confusion Over Third-Order Accuracy of U-MUSCL[END_REF], U-MUSCL with κ = 1/3 can be third-order accurate in cell-averaged solutions on regular grids for one-dimensional problem, and in the pointvalued solution on regular grids for linear equations in all dimensions. However, only second-order can be obtained in our configuration, namely the two-dimensional case on unstructured grid.

Gradient computation on unstructured meshes

The MUSCL-type scheme presented previously requires the computation of the gradient of the unknowns on unstructured meshes. There are two methods main for gradient computation: the Green-Gauss method (GG) and the Least SQuares method (LSQ). The first one is based on the application of the Green-Ostrogradski theorem. Although easy to implement and inexpensive, it is not very accurate on deformed meshes [START_REF] Syrakos | A critical analysis of some popular methods for the discretisation of the gradient operator in finite volume methods[END_REF]. Thus, the LSQ method is preferred here for its robustness and accuracy. It is based on the Taylor development of the solution within a cell i on a stencil v(i) composed of its neighbors. It can be composed of neighbors per edge, neighbors per vertex, or from an even more extended stencil. These first two configurations are represented in Fig. 3. An enlarged stencil for gradient computation enables greater accuracy, especially when the mesh is highly distorted.

The LSQ method consists in solving an over-determined system AX = b where X = ∇α i (see Appendix A). By construction, the method is exact for any linear function, regardless of the mesh. This is essential to obtain secondorder accuracy with MUSCL-type schemes. The accuracy of the gradient computation can be improved by using a weighted matrix W , usually diagonal, with coefficients ω kk such as

ω kk = ||x j k -x i || -q , (3.6)
where q is a positive parameter. When q = 0, we find the unweighted method. With q = 1, the weight corresponds to the inverse of the distance between the mass center of cell i and cell j. In this work, q is set to 3/2 [START_REF] Syrakos | A critical analysis of some popular methods for the discretisation of the gradient operator in finite volume methods[END_REF]. We thus need to resolve W AX = W b. This over-determined system is solved by going through the normal equations, i.e., by multiplying the system by (W A) T , and the solution reads

• x i • x j1 • x j2 • x j3 • x j4 • x j5 • x j6 • x j7 • x j8
X = (A T W T W A) -1 A T W T W b. (3.7)
In the case where the matrix is ill-conditioned however, the gradient calculation could be imprecise. Another approach to avoid such a problem is to use the QR decomposition. The matrix W A can be decomposed as W A = QR where Q is an orthogonal matrix and R is an upper triangular matrix. The solution is then written

X = R -1 Q T W b. (3.8)
The second method is thus preferred. The QR decomposition is performed with the modified Gram-Schmidt algorithm [START_REF] Gander | Algorithms for the QR-Decomposition[END_REF]. Note that the matrix R -1 Q T W consists only of geometrical parameters. When the mesh is not moving, it is then sufficient to calculate it once when initializing the calculations and to store it. This method can give secondorder accurate gradient reconstruction for smooth curvilinear mesh or on meshes of identical parallelograms [START_REF] Syrakos | A critical analysis of some popular methods for the discretisation of the gradient operator in finite volume methods[END_REF].

Gradient limiter

In order to eliminate oscillations inherent to high-order schemes and capture discontinuities, these reconstructions are limited by a function of the local gradient, called "limiter". Many limiters exist in the literature, designed to have different properties. In the case of the five-equation system, the volume fraction is a discontinuous variable, which follows a transport equation. Following the lines of [START_REF] Qian | A freesurface capturing method for two fluid flows with moving bodies[END_REF][START_REF] Blanchard | Modélisation et simulation multi-échelles de l'atomisation d'une nappe liquide cisaillée[END_REF][START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF], the idea here is to apply a compressive limiter to this variable, which allows a better representation of the interface between the two fluids. For the other variables, namely densities, velocity and pressure, the limiter should be as robust as possible, given the implicit treatment of the acoustic step.

The reconstruction αn ij is rewritten in the form

αij = α i + φ i ∆α i,j , (3.9)
where φ i is the gradient limiter, and

∆α i,j = κ 2 (α p -α i ) + (1 -κ)∇α i .(x f -x i ). (3.10)

Proposition 1

Let φ i the limiter of the reconstruction of α i in cell i, such that

φ i =      min β, α i -α min i |∆α i | max , α max i -α i |∆α i | max if |∆α i | max = 0, 0 otherwise. (3.11) with |∆α i | max = max j∈v(i) |∆α i,j |, α min i = min j∈v(i) (α i , ᾱj ) and α max i = max j∈v(i) (α i , ᾱj ).
Depending on the parameters β and ᾱj , this limiter displays different properties :

• If β = 1 and ᾱj = α i ω i + α j ω j ω i + ω j with ω i = ||x j k - x i || -1 then the limiter is monotonicity preserving, i.e. ∀j ∈ v(i), α i ≤ α j =⇒ αij ≤ αji . (3.12)
• If β > 1 and ᾱj = α j , the limiter is compressive.

Proof. The design of the previous limiter is provided in Appendix B.

This approach is used in [START_REF] Qian | A freesurface capturing method for two fluid flows with moving bodies[END_REF] in the framework of onedimensional limiter for two-phase flows, where the compressive limiter is used on the density whereas the minmod limiter is used on the velocity and the pressure. In [START_REF] Blanchard | Modélisation et simulation multi-échelles de l'atomisation d'une nappe liquide cisaillée[END_REF], a similar β-limiter is designed without the monotonicity preserving condition. The compressive limiter (3.11) with β = 2 is used on the volume fraction of a two-phase model while the less compressive limiter with β = 1 is used on the others variables. In the work presented in [START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF], a family of β-limiters is designed from the approach of [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] where as previously, the compressive limiter is used on the volume fraction only.

Monotonicity preserving. When dealing with large pressure or density ratio, as in the study of liquid-gas interaction, the monotonicity preserving property can be essential, in order to improve robustness. The well-known minmod limiter verifies this condition, which explains why it is very popular in industrial codes.

Compressive property. The compressive nature of a method or limiter reflects the ability to accurately reproduce a discontinuity. For example, the superbee limiter [START_REF] Toro | Riemann solvers and Numerical Methods for Fluid Dynamics[END_REF] is wellknown to sharpen sinusoidal profiles and to behave well on discontinuity transport. This property is particularly valuable in the context of immiscible fluid flows, as it improves the representation of the interface between these two fluids [START_REF] Qian | A freesurface capturing method for two fluid flows with moving bodies[END_REF][START_REF] Blanchard | Modélisation et simulation multi-échelles de l'atomisation d'une nappe liquide cisaillée[END_REF][START_REF] Chiapolino | Sharpening diffuse interfaces with compressible fluids on unstructured meshes[END_REF]. When the limiter (3.11) is allowed to take values greater than 1, it is then compressive.

Remark 3

When a second-order in space scheme is used, the maximum time step given by CFL conditions (2.15) and (2.23) must be divided by two (see [START_REF] Berthon | Stability of the MUSCL Schemes for the Euler Equations[END_REF] for example).

3.1.3. Second-order in space scheme in the splitting strategy The limiter (3.11) is used in the context of the fiveequation system, while reconstructions are performed with the U-MUSCL scheme with κ = 1/3. The objective is to sharpen the fluid's discontinuities by reducing the number of mixture cells. The reconstructed variables are the primitive variables, in order to limit velocity and pressure oscillations. They are reconstructed as follows Acoustic step: (u, p) n are reconstructed with β = 1.

Transport step: (ρ 1 , ρ 2 , u, p) # are reconstructed with β = 1 and z # with β = 2. The conservative variables are deduced from the primitive ones with

                     ρ = zρ 1 + (1 -z)ρ 2 , ρy = zρ 1 , ρu = (zρ 1 + (1 -z)ρ 2 ) u, ρε = k z k ρ k ε k = k z k p + γ k π k γ k -1 , ρE = ρε + 1 2 ρu • u. (3.13)

Remark 4

The limiter needs the computation of α min i and α max i , the local extrema. The 5-point stencil, i.e., with the edges neighbors, is used for the acoustic step, while the 9-point stencil i.e., with the vertex neighbors is utilized for the transport step (see Fig. 3). Better convergence properties have been observed numerically on the acoustic step with only 5 points.

Time-scheme 3.2.1. The Heun method

The time scheme is important to increase the accuracy of the overall scheme, especially when high-order in space schemes are used. This prevents the appearance of numerical artifacts that will degrade the solution, as shown in [START_REF] De Vuyst | A geometrically accurate low-diffusive conservative interface capturing method suitable for multimaterial flows[END_REF] for instance. To this end, we could use the Heun method also called Runge-Kutta TVD, which we apply on the global scheme (2.24). It reads

U * = U n + ∆tF(U n , ∆t), U * * = U * + ∆tF(U * , ∆t), U n+1 = 1 2 (U n + U * * ) , (3.14)
with U being the conservative variables (ρ, ρy, ρu, ρε, z) and F is a numerical flux discretized in (2.24). It is not clear that this scheme is second-order in the case of the splitting strategy because of the time step dependence of the numerical flux (see Rem. 2). Indeed, the numerical flux is a first-order in time approximation of the spatial operator. However, it does improve the numerical solution, which will be shown in the results section.

An alternative implicit time-scheme

The scheme described above is simple to implement and robust. When an implicit approach is used in the acoustic step however, it must be repeated twice per time step, which induces a significant cost in computation time. For this, a new implicit scheme based on the Crank-Nicolson scheme has been developed. The modified scheme proposed here reads

V † -V n ∆t = - 1 2 ϑ n ∇ • G(V † ) + ∇ • G(V n ) . (3.15)
This scheme is not the Crank-Nicolson scheme because it is not second-order since we kept ϑ n for the explicit and implicit terms. It can be seen as the average of the flux given by the explicit (2.13) and the implicit (2.13) versions of the acoustic step. The truncation error of this modified Crank-Nicolson scheme is smaller than that of the first implicit scheme (see Appendix C). The numerical scheme for the acoustic step reads

V † = V n - ∆t |Ω i | ϑ n i j∈v(i) |Γ ij | 1 2 G n ij + G † ij . (3.16)
The resolution of (3.16) follows the same line as the firstorder implicit time scheme. In the one-dimensional case, it is easy to show that the new scheme leads to solve under the delta-form the system

I d + ∆t 2∆x M X † -X n = - ∆t ∆x M X n . (3.17)
From an implementation point of view, this represents a small modification of the velocity-pressure system (2.16) resolution algorithm. Once again, the other equations are solved explicitly. In order to have a globally conservative scheme, the material velocity used in the transport step must read

u ij = ūn ij + ū † ij 2 . (3.18)
This comes from the same arguments used in [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF] in order to have a globally conservative scheme.

Extension with the second-order in space scheme

When using the second-order scheme in space, nonlinearities introduced by the limiter make it impossible to write the flux jump as a matrix-vector product of X. A conventional way to solve non-linear equations is to use a quasi-Newton method. Using this method, the solution is obtained by iterating over the following equation

I d + ∆t ∆x M X k+1 -X k = -X k + X n - ∆t ∆x ϑ n [[H k ]], (3.19)
with k = 0, ..., ∞. The matrix M is built with the firstorder flux while the right-hand side of (3.19) is computed with the second order flux. Note that a single iteration over (3.19) gives exactly (2.20), the system solved when first-order in space scheme is used. Through this method, the implicit time scheme converges to the second-order in space solution. When using the second implicit scheme, equation (3.19) becomes

I d + ∆t 2∆x M X k+1 -X k = -X k + X n - ∆t 2∆x ϑ n [[H k ]] + [[H n ]] . (3.20)

Remark 5

In practice, we use a convergence threshold set to 10 -5 with a maximum of 10 iterations. These values have proven to be sufficient in the simulations performed.

Numerical results

In this section, a selection of numerical results for the one and two-dimensional cases are presented. The first cases will illustrate the accuracy of the compressive limiter (β = 2) on the classical pure transport equation. The RK2-TVD time scheme presented in the previous section is always used for pure transport cases.

The next cases will show the ability of the second-order method to simulate complex two-phase flows. We denote by EXEX the explicit scheme for both the acoustic and transport steps, while IMEX is the implicit-explicit strategy of the splitting scheme. When adding RK2, the RK2-TVD scheme is thus used on the overall algorithm. The CNEX scheme denotes the new implicit-explicit scheme built from a modification of the Crank-Nicolson scheme. For two-dimensional liquid-gas or gas-gas cases, the density gradient will be computed in log scale during postprocessing, as log(|∇ρ|), in order to reproduce numerical Schlieren.

Zalesak's disk

We consider here the classical Zalesak's disk test case [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF], namely a disk with a slot rotating in the velocity field u(x, y) = (0.5 -y, x -0.5) m.s -1 on the domain The present limiter gives very good results. The interface is very sharp and the shape of the disk is wellpreserved, almost symmetrical. The use of κ = 1/3 (see Fig. 4.a) gives a good resolution of the slot and the corners, compared to the conventional MUSCL method with κ = 0 (see Fig. 4.b). This might be because the method with κ = 1/3 is third-order accurate on the transport equation on Cartesian mesh.

Ω = [0, 1] m × [0, 1] m. The disk is defined at the initial time by z(x, y) =      1 if (||x -c||< r) and

Kothe-Rider advection

We now consider the Kothe-Rider forward-backward advection test case [START_REF] Rider | Reconstructing Volume Tracking[END_REF] on the domain Ω = [0, 1] m × [0, 1] m. This test case consists in transport a circle in a backward-forward velocity field such that u(x, y, t) = -sin(πy) cos(πy) sin(πx) 2 cos(πt/t f ) sin(πx) cos(πx) sin(πy) 2 cos(πt/t f ) .

The circle of radius R = 0.15 m is set at (0.5, 0.75) m. The CFL is set to 0.2 and 256 × 256 cells are used with constants boundary conditions. The results in Fig. 5 are obtained at the intermediate time t = 3s and at the final time t f = 6s. At first, until t = 3s, the velocity field deforms the disk into a filament. The filament obtained (see Fig. 5.a) is wellcaptured in regions where it is large enough regarding the cell's size. In a second time, the velocity field rewinds the filament to reform a disk. At the end, the shape of the disk is quite well-preserved (see Fig. 5.b). Some diffusion from the forward part deteriorates the results. A mesh with more cells gives a better disk resolution.

One-dimensional liquid-gas shock tube

We consider here the classical one-dimensional liquidgas shock tube test on Ω = [0, 1] m, studied in [START_REF] Vilar | Positivity-preserving cell-centered lagrangian schemes for multi-material compressible flows: From first-order to high-order. Part I: The onedimensional case[END_REF] for instance. The stiffened gas EOS is used for the liquid, while the gas is modelled by the perfect gas EOS. The initial conditions read

(ρ, u, p, γ, π) = (10 3 , 0, 10 9 , 4.4, 6 × 10 8 ) if x ≤ 0.7 m, (5, 0, 10 5 , 1.4, 0) else.
The CFL number for the transport step is set to 0.25, and we use 400 cells with constant boundary conditions. The final time is set to t f = 240µs. This case serves to test the accuracy of the different time schemes to represent the three waves of the problem, that are rarefaction wave, contact discontinuity and shock. This is a challenging test case, since large pressure and density ratios are involved. The results are presented in Fig. 6. The results are in very good agreement with the exact solution. The three waves are well-represented, especially the contact discontinuity, thanks to the compressive limiter used on the volume fraction. One can see an undershoot on contact discontinuity of the density (see Fig. 6.a). This observation is quite common when sharp interface techniques are used, as the anti-diffusive scheme [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF][START_REF] Friess | A second order anti-diffusive Lagrange-remap scheme for two-component flows[END_REF]. Regarding velocity and pressure, the CNEX time scheme gives more accurate results than the IMEX scheme or the IMEX+RK2 strategy, especially near the rarefaction wave. It also reduces the overshoot on the density and better localizes the shock position.

Two-dimensional Air-R22 interaction

We shall now consider the experiment of Haas & Sturtevant [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF]. In this test case, a shock wave at Mach 1.22 propagating through the air hits an R22 gas cylinder of radius R = 25 mm located at (130 mm, 44.5 mm). Many simulations of this case can be found in the literature, such as [START_REF] Shyue | A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions[END_REF][START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF][START_REF] Helluy | OpenCL numerical simulations of two-fluid compressible flows with a 2D random choice method[END_REF]. The computational domain Ω = [0, 200] mm × [0, 89] mm is described in Fig. 7.

The initial states of both fluids in the pre-and postshock regions are given in Tab. 1. Both fluids are modeled by perfect gas EOS with γ = 1.4 for the air and γ = 1.249 for the R22. The pre-shock density, velocity, and pressure are computed with the shock wave relation (see [START_REF] Toro | Riemann solvers and Numerical Methods for Fluid Dynamics[END_REF] for instance) such that the Mach shock is 1.22. Each phase is modeled by perfect gas EOS.

Location

ρ (kg.m -3 ) u (m.s -1 ) p (10 Because of the symmetry of the problem, we use a domain of Ω = [0, 200] mm × [44.5, 89] mm with boundary symmetric condition at the lower side, a wall boundary condition at the upper side, an outflow boundary conditions on the left side and a supersonic inlet boundary condition on the right side. The mesh is composed of 1600 × 356 cells. In this case, the implicit time scheme is not useful since the ratio of densities and pressures is approximately 1. The explicit version of the time scheme is therefore used. The CFL number for both steps is set to 0.25. The final time is t f = 1020µs.

The numerical Schlieren of this case is presented in Fig. 8. For the lower half, the results are obtained with the EXEX time-scheme, while the EXEX+RK2 strategy is used for the upper half. For both simulations, the second-order in space scheme with compressive limiter is used. When the shock impacts the R22 gas, the cylinder is deformed. This velocity difference, combined with the interaction of the interface with the transmitted shock waves produce hydrodynamic instabilities, and curling vortices are generated at the interface between the two gases. The scheme is able to reproduce a large scale and finer flow structure, while maintaining a very sharp interface throughout the simulation. The results are in good agreement with the experiments of [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF] and with numerical results from the literature [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF][START_REF] Kreeft | A new formulation of Kapila's five-equation model for compressible two-fluid flow, and its numerical treatment[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluide[END_REF][START_REF] Cheng | Lowdissipation BVD schemes for single and multi-phase compressible flows on unstructured grids[END_REF]. The use of the EXEX-RK2 scheme appears to attenuate the spurious waves that appear with the EXEX scheme.

Two-dimensional liquid-gas interaction

We shall now consider a two-dimensional liquid-gas interaction, widely studied in the literature [START_REF] Saurel | A simple method for compressible multifluid flows[END_REF][START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluide[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. The rectangular computational domain Ω = [0, 2] m × [0, 1] m is described in Fig. 9. A gas cylinder of radius R = 0.4 m is initially at rest in a liquid at (0.5, 0.5) m while a shock at x = 0.04 m propagates in the liquid from left to right.

The initial conditions of each phase are provided in Tab. 2. The air is modeled by perfect gas EOS with γ = 1.4 while the liquid is modeled by stiffened gas EOS with γ = 4. Table 2: Two-dimensional liquid-gas interaction : initial data condition at the lower side, a wall boundary condition at the upper side and inflow/outflow boundary conditions on the other sides. The mesh is composed of 1600 × 400 cells.

The CNEX time scheme is used with the second-order in space strategy. Since the acoustic step is implicit, the time step is computed with the transport condition (2.23) with a CFL number set to 0.25. Volume fraction and numerical Schlieren are presented on Fig. 10. The results are in good agreement with the results presented in [START_REF] Saurel | A simple method for compressible multifluid flows[END_REF][START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of five-equation model[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluide[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF] shock wave hits the gas cylinder and compresses it while propagating inside. When the shock reaches at the left of the cylinder, the gas is split into two different cylinders.

Remark 6

Many schemes are unable to cope with this test case because of the high ratio of pressure and density, even with explicit time schemes. The robustness of the method come from the exact conditions on the Riemann solver slopes to ensure positivity of the intermediate states [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for twophase flows[END_REF]. However, these conditions are derived from the explicit time-scheme of the acoustic step. In our case, when an implicit timescheme is used, no conditions can be derived to guarantee that the solution remains in the convex set of admissible solutions, which ensures that the sound speed given by (2.7) is real, and the mixture density is positive. When using the compressive limiter however, it happens in this test case that p + π becomes negative in mixture cells where the volume fraction of liquid is approximately 10 -4 and the mixture density is about 10 -3 . To overcome this issue, when the square of the velocity of the mixing sound speed in cell i is negative, it is calculated with

c 2 i = max(c 2 i,liq , c 2 i,gas ) (4.1)
where c i,liq and c i,gas are the sound speed of the liquid phase and gas phase in cell i. This fix can be seen as adding diffusion to the implicit resolution of the acoustic step, where the slopes of the Riemann solver are locally increased. It is the only part where the sound speed is needed.

4.6. Two-dimensional liquid-gas interaction on deformed mesh All the cases presented above were performed on Cartesian grids. Thus, the last liquid-gas test case is reproduced on a deformed mesh, in order to show that the method works on unstructured grids. The deformation is performed from a Cartesian mesh, where the nodes are randomly moved. A small portion of the mesh is shown in Fig. 11. The same parameters, scheme and number of cells are used. This time, the calculation is performed on the complete domain, without using any symmetry conditions.

The numerical Schlieren at several times are presented in Fig. 12. The results are in good agreement with the previous simulation on a Cartesian grid. One can observe that the shape of the gas cylinders are not perfectly symmetrical with the horizontal axis. This comes from the non-symmetry of the randomly deformed grid. The different waves are, however, well-captured, and the two-phase interface remains sharp.

Conclusion

In this work, several methods for reducing numerical diffusion when studying compressible two-phase flows with immiscible fluids were presented and combined. Since nonmiscible fluids are considered, the five-equation system of [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF] was used. A robust implicit-explicit acoustictransport splitting scheme for two-phase flows is used to solve the governing equations. The main contribution of this work is to improve the accuracy of numerical methods used on meshes potentially distorted by ablation. A new and more accurate implicit-explicit time scheme was developed. The second-order extension of the method was performed with a MUSCL-type scheme, that serves to improve the resolution of the numerical scheme. The U-MUSCL scheme [START_REF] Burg | Higher Order Variable Extrapolation for Unstructured Finite Volume RANS Flow Solvers[END_REF] with κ = 1/3 is used, with the correction of [START_REF] Nishikawa | On the loss and recovery of second-order accuracy with U-MUSCL[END_REF] for deformed grids. A multidimensional limiter is designed to capture shock and contact discontinuity without oscillations. A β parameter serves to control the compressive property of the limiter, used only on the volume fraction. It sharpens the interface between phases, reduce the numerical diffusion and reproduces finer flow structures. The continuation of this work consists in enriching the model by taking into account capillary effects. These physics will be included in the acoustic step, still with an implicit time-scheme.

Appendix A. Least-squares method

The least-squares method used to compute the local gradient is described in this appendix. We consider the solution α i in x i , the center of the cell i. We then define the following linear model αi (x) = α i + ∇α i .(x -x i ).

(A.1)

The method consists in finding the minimum of the sum

S i = j∈v(i) ( αi (x j ) -α j ) 2 , (A.2)
where v(i) is a set of neighbors of cell i. The minimization problem (A.2) can be written in the form

S i = ||AX -b|| 2 (A.3) with X =    ∂α i ∂x ∂α i ∂y   , A =      x j1 -x i y j1 -y i x j2 -x i y j2 -y i . . . . . . x j N -x i y j N -y i      and b =      α j1 -α i α j2 -α i . . . α j N -α i      .
where N is the number of neighboring cells. Thus, finding the minimum of (A.2) consists in finding the minimum of the function f (X) = ||AX -b|| 2 . Its gradient is zero when AX = b.

Appendix B. Design of the limiter

In this appendix, the main lines used to obtain the limiter from the inequalities on the local monotonicity preserving principle are detailed.

First, the monotonicity principle is derived. To satisfy the monotonicity principle, the reconstruction of α i need to hold between min j∈v(i) (α i , ᾱj ) and max j∈v(i) (α i , ᾱj ) where ᾱj is the average value between α i and α j on their common face. It is such that

ᾱj -α i ||x j k -x i || = α j -ᾱj ||x j k -x j || . (B.1)
Thus, ᾱj = α i ω i + α j ω j ω i + ω j (B.2)

with ω i = ||x j k -x i || -1 .

Lemma 1

Let αij be the extrapolation of α i on Γ ij such that αij = α i + φ i ∆α i,j , (B.3)

where φ i is the gradient limiter, and Proof. From B.5, we have ∀j ∈ v(i)

∆α i,j = κ 2 (α p -α i ) + (1 -κ)∇α i .(x f -x i ). (B.
   α min i ≤ α i + φ i ∆α i,j ≤ α max i , α min i ≤ α i -φ i ∆α i,j ≤ α max i , which gives    α min i -α i ≤ φ i ∆α i,j ≤ α max i -α i , α i -α max i ≤ φ i ∆α i,j ≤ α i -α min i .
Then, two cases should be distinguished

• if ∆α i,j > 0 then

       α min i -α i ∆α i,j ≤ φ i ≤ α max i -α i ∆α i,j , α i -α max i ∆α i,j ≤ φ i ≤ α i -α min i ∆α i,j ,
• if ∆α i,j < 0 then

       α i -α max i |∆α i,j | ≤ φ i ≤ α i -α min i |∆α i,j | , α min i -α i |∆α i,j | ≤ φ i ≤ α max i -α i |∆α i,j | .
Because φ i has to be positive, the left inequalities are always true. The right inequalities have to be true ∀j ∈ v(i) so the minimum value is chosen. The limiter thus reads The constraint φ i ≤ 1 must be added in order to ensure the reconstruction to be exact for linear solution on admissible mesh. The final expression of limiter is thus

φ i = min j∈v(i)              min α max i -α i ∆α i,
φ i =      min 1, α i -α min i |∆α i | max , α max i -α i |∆α i | max
if ∆α i,j = 0, 0 otherwise, (B.8) Secondly, if the limiter must be compressive, the only condition that need to be satisfy is the maximum principle ∀j, min j∈v(i) (α i , α j ) ≤ αij ≤ max j∈v(i) (α i , α j ).

(B.9)

The limiter is thus

φ i =      min β, α i -α min i |∆α i | max , α max i -α i |∆α i | max
if ∆α i,j = 0, 0 otherwise, (B.10) with α min i = min j∈v(i) (α i , α j ) and α max i = max j∈v(i) (α i , α j ). The β parameter is added to control the compressive property of the limiter. For the volume fraction which is theoretically piecewise constant, preservation of the maximum principle is the only desired property. It is not necessary to preserve the linearity of the solution so β > 1 can be employed.

Appendix C. Truncation error of implicit acoustic time-schemes

In this appendix, the calculation of the truncation error of both implicit time-scheme used for the acoustic step (2.10) is detailed.

Lemma 2

The truncation error of the implicit time-scheme (2.13) that reads

V n+1 -V n ∆t = -ϑ n ∇ • G(V n+1 ), (C.1)
is given by

τ n = ∆t ∂ 2 t V (t) 2 + ϑ(t)∇ • ∂ t G(V (t)) + O(∆t 2 ), (C.2)
while the truncation error of the modified Crank-Nicolson time-scheme (3.15) that reads From the truncation errors calculated earlier, one can see that the new implicit scheme, called the modified Crank-Nicolson scheme, produces a smaller error than the first implicit scheme. We then expect to achieve more accurate numerical results with the second scheme.

V n+1 -V n ∆t = - 1 2 ϑ n ∇ • G(V n ) + ∇ • G(V n+1

Figure 1 :

 1 Figure 1: Schematic representation of an object entering the atmosphere.

Figure 2 :

 2 Figure 2: Notations associated with the unstructured mesh.

Figure 3 :

 3 Figure 3: Schematic representation of the neighborhood by edge (gray) and the neighborhood by vertex (white and gray) for the calculation of the gradient on cell i by the LSQ method.

(

  |x -cx |≥ e or y -c y > l) , 0 otherwise, with x = (x, y), c = (0.5, 0.75) m, r = 0.15 m, e = 0.025 m and l = 0.1 m. The CFL is set to 0.4 and 256×256 cells are used with constant boundary conditions. The results in Fig.4are obtained after one revolution at the final time t f = 2πs.

Figure 4 :

 4 Figure 4: Zalesak's disk test case. Initial condition (a), solution after one rotation with κ = 1/3 (b) and with κ = 0 (c).

Figure 5 :

 5 Figure 5: Kothe-Rider advection test case. Initial condition (a), solution obtained at the intermediate time t = 3s (b) and solution at the final time t = 6s (c).

  4 and π = 6.8 × 10 8 . Because of the symmetry of the problem, a domain of Ω = [0, 2] m×[0.5, 1] m is used with a boundary symmetric Location ρ (kg.m -3 ) u (m.s -1 ) p (Pa) Liquid (shock) 1030.9

Figure 6 :

 6 Figure 6: One-dimensional liquid-gas shock tube test case. The density (upper left), velocity (upper right), pressure (lower left) and volume fraction (lower right) are compared at the final time t f = 240µs with the exact solution. Except for the volume fraction, only one cell in 10 is represented by points.

Figure 7 :

 7 Figure 7: Two-dimensional Air-R22 interaction : geometric description of the test case.

  (a) t = 0µs. (b) t = 204µs. (c) t = 408µs.

  (d) t = 612µs. (e) t = 816µs. (f) t = 1020µs.

Figure 8 :

 8 Figure 8: Numerical Schlieren of the air-R22 interaction test case. Results of the EXEX scheme (lower half) and results with the EXEX+RK2 scheme (upper half) at several times.

Figure 9 :

 9 Figure 9: Two-dimensional liquid-gas interaction : geometric description of the test case.

  (a) t = 100µs. (b) t = 200µs. (c) t = 300µs. (d) t = 400µs. (e) t = 500µs. (f) t = 600µs.

Figure 10 :

 10 Figure 10: Two-dimensional liquid-gas interaction test case. Numerical Schlieren (upper half) and volume fraction (lower half) at several times.

Figure 11 :

 11 Figure 11: Mesh used for the Two-dimensional liquid-gas interaction on deformed mesh.

  (a) t = 100µs. (b) t = 200µs. (c) t = 300µs. (d) t = 400µs. (e) t = 500µs. (f) t = 600µs.

Figure 12 :

 12 Figure 12: Two-dimensional liquid-gas interaction test case on random mesh. Numerical Schlieren at several times.

4 ). (B. 5 )

 45 The reconstruction (B.3) and its symmetrical preserve the monotonicity of the solution i.e.∀j ∈ v(i), α min i ≤ α i ± φ i ∆α i,j ≤ α max i where α min i = min j∈v(i) (α i , ᾱj ), α max i = max j∈v(i) (α i , ᾱj ), the limiter φ i reads |∆α i | max if ∆α i,j = 0, 0 otherwise, (B.6) with |∆α i | max = max j∈v(i) |∆α i,j |.

Table 1 :

 1 Two-dimensional Air-R22 interaction : initial data.

	5 Pa)

  . The compressive limiter gives a sharp interface of the volume fraction. The propagation waves are shown on numerical Schlieren. The

		1000								solution IMEX			500	solution IMEX
									IMEX RK2				IMEX RK2
		800								CNEX			400	CNEX
	ρ (kg.m -3 )	400 600	20 30 40 50									u (m.s -1 )	200 300			488 490 492 494
		200	10										100			486
																	0.35	0.526	0.702	0.878
			0
		0	0.8		0.826	0.852	0.878					0
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						x (m)										x (m)
						(a)										(b)
	p (Pa)	0 2e+08 4e+08 6e+08 8e+08 1e+09				-1e+06 0 1e+06 2e+06 0.28 3e+06	0.433	0.586 solution 0.739 IMEX IMEX RK2 CNEX	z	0.2 0.4 0.6 1 0.8	0 0.2 0.4 0.6 0.8 1				solution IMEX IMEX RK2 CNEX
													0	0.814	0.817	0.82		0.823
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						x (m)										x (m)
						(c)										(d)

  (t) + ϑ(t)∇ • G(V (t) +∆t∇ • (∂ t G(V (t))) + O(∆t 2 ),

					) ,
					(C.3)
	is given by
	τ n = ∆t	∂ 2 t V (t) 2	+	ϑ(t)∇ • ∂ t G(V (t)) 2	+O(∆t 2 ). (C.4)
	Proof. The truncation error of the first implicit time-scheme
	scheme is
	τ n = = ∂ t V (t) + V (t + ∆t) -V (t) ∆t ∆t 2 ∂ 2 t V = ∆t ∂ 2 t V (t) 2 + ϑ(t)∇ • ∂ (C.5) + ϑ(t)∇ • G(V (t + ∆t)),
	The truncation error of the modified Crank-Nicolson
	scheme is
	V (t + ∆t) -V (t) ∆t + 1 ϑ(t) (∇ • G(V (t + ∆t)) + ∇ • G(V (t))) , 2 = ∂ t V (t) + τ n = ∆t ∂ 2 t V (t) + ϑ(t)∇ • G(V (t) 2 + ∆t ∇ • (∂ t G(V (t))) + O(∆t 2 ), 2 = ∆t 2 ∂ 2

t G(V (t)) + O(∆t 2 ). t V (t) + ϑ(t)∇ • ∂ t G(V (t)) + O(∆t 2 ). (C.6)


Mohammed H Geesi

Yassine Riadi
email: y.riadi@psau.edu.sa

Abdellah Kaiba

El Elmutasim O Ibnouf

Hassane Anouar

Oussama Dehbi

Saïd Lazar

Philippe Guionneau

Elhassane Anouar

of 4-[2-(1-methyl-1H-imidazol-2-ylsulfanyl)-acetylamino]benzenesulfonic acid

Keywords: Synthesis, antimicrobial evaluation, crystal structure, Hirschfeld surface analysis and docking studies Synthesis, Antimicrobial Evaluation, Crystal Structure, Hirschfeld Surface Analysis synthesis, antimicrobial activity, docking, crystal structure, Hirshfeld surface analysis 1

may come

Introduction

Heterocyclic chemistry is a vast field, considering that the number of heterocyclic compounds continues to increase. Among the different classes of heterocyclic compounds, nitrogenous structures are the most common and are present in many natural compounds of plant, animal and synthetic origin. Heterocyclic structures are sometimes associated with each other, but are, in most cases, linked to very diverse structural scaffolds. Some hybrid compounds, comprising mainly heterocycles containing nitrogen, sulphur and/or oxygen atoms, have shown remarkable pharmacological activity [1][2][3][4].

For example, the imidazole nucleus is responsible for the synthesis of a large number of active agents used in various fields. In the medical and therapeutic fields, it is crucial for anticancer treatments [5,6], antihistamines [7], anti-inflammatories [8], bactericides [9], fungicides [10,14] anti-hypertensives [START_REF] Faber | Bio-transformations in Organic Chemistry[END_REF][START_REF] Yan | [END_REF][13][14], anti-tuberculosis drugs [15], anti-parasitic agents [16,17], etc. It is also used in the agrochemical sector [START_REF] Katritzky | Comprehensive Hetrocyclic Chemistry II[END_REF][START_REF] Reis | [END_REF] and in the industry sector (e.g., catalysts, flame retardants, steel corrosion inhibitors, etc.) [START_REF] Seddon | [END_REF][21][22][23].

This research was performed because of our experience with the synthesis and characterisation of novel active heterocyclic compounds [START_REF] Riadi | Ecofriendly synthesis of 3cyanopyridine derivatives by multi-component reaction catalyzed by animal bone meal[END_REF][START_REF] Dehbi | Water-mediated synthesis of disubstituted 5-aminopyrimidines from vinyl azides under microwave irradiation[END_REF][START_REF] Ishak | An efficient synthesis and antimicrobial and antifungal activities of disubstituted 3, 4-dihydro-2H-1, 3-thiazin-4-ones using lemon juice: A natural approach[END_REF][START_REF] Riadi | An efficient and reusable heterogeneous catalyst animal bone meal for facile synthesis of benzimidazoles, benzoxazoles, and benzothiazoles[END_REF][START_REF] Dehbi | Synthesis and optimization of an original V-shaped collection of 4-7disubstituted Pyrido [3, 2-d] pyrimidines as CDK5 and DYRK1A inhibitors[END_REF][START_REF] Riadi | Green, rapid and efficient synthesis of new antibacterial pyridopyrimidinone mediated by eutectic mixture of Urea/CuCl2[END_REF][START_REF] Ouerghi | Synthesis, antibacterial evaluation, crystal structure and molecular interactions analysis of new 6-Bromo-2-chloro-3-butylquinazolin-4(3H)one[END_REF][START_REF] Geesi | Synthesis, antibacterial evaluation, Raman, Crystal Structure and Hirshfeld Surface analysis of a new 3-(4-fluorophenyl)-6-methyl-2-(propylthio) quinazolin-4 (3H)-one[END_REF]. The main objective was the preparation, biological evaluation and crystal structural identification of a new, imidazolecontaining molecule: 4-[2-(1-methyl-1H-imidazol-2-ylsulfanyl)-acetylamino]-benzenesulfonic acid (hereafter, "5").

Methodology and Results

Synthesis

First, an intermediate compound (3) was prepared from the reaction of (1) 2-mercapto-1methylimidazole (1 mmol) and (2) ethyl 2-bromoacetate (1 mmol) in the presence of potassium carbonate (0.5 mmol) under refluxed acetone for 10 h. The intermediate compound (3) was obtained after filtration and recrystallisation in ethanol. Then, 1 mmol of the intermediate compound (3) and 1.1 mmol of benzenesulfonic acid were refluxed in the presence of 20 mg of trimethylamine and 10 mL of ethanol for 3 h. After the reaction's completion (monitored using TLC), the mixture was filtered and the resulting solid was recrystallised in a mixture of ethanol and water to achieve a pure compound (5) with an 83% yield (Scheme 1).

Antimicrobial Performance

Herein, the minimum inhibitory concentration (MIC) is referred to as the concentration of 5, for which bacterial growth was completely inhibited. MIC was determined by two-fold microdilution in sterile flat-bottom 96-well polystyrene plates, as detailed in [START_REF] Andrews | Determination of minimum inhibitory concentrations[END_REF].

Dimethyl sulfoxide (DMSO) was used as the negative control. The MIC was determined as follows: first, concentrations of 5 in DMSO were adjusted to 100, 75, 50, 38.5, 25, 18.75 and 12.5 μg/mL. Second, 200 μL of the former serially diluted test samples were dispensed to 96well microtrays. Next, test bacteria were inoculated into the wells by adjusting the final volume to 400 μL, followed by the incubation of the plates at 37 °C for 24 h. The concentration, for which a clear reduction in turbidity was observed, was referred to as the MIC. All experiments were carried out in triplicate.

The prepared 5 was tested for antimicrobial efficiency against different microbial pathogens, such as methicillin-resistant Staphylococcus aureus (MRSA), Klebsiella pneumoniae (Kp), Aspergillus niger (An) and Candida albicans (Ca). As shown in Table 1, the synthesised derivative was found to be effective against all tested microbes.

X-ray Diffraction Analysis

Single-crystal data for C 12 H 13 N 3 O 4 S 2 were collected on a Bruker Apex II diffractometer equipped with a molybdenum anticathode and graphite monochromator. The colourless crystal compound had a plate form about 0.45 × 0.12 × 0.02 mm 3 in size. Absorption corrections were empirical. The ϕ and ω scan modes were used and the crystal-detector distance was 35 mm. The DENZO-SMN programme was used to reduce the data. The SIR97 programme was used to determine a structural hypothesis [START_REF] Altomare | SIR97: a new tool for crystal structure determination and refinement[END_REF]. The structure solution was solved by the SHELX97 programme, and refinement was performed on F 2 via the full-matrix least squares technique using the SHELXL programme [START_REF] Sheldrick | SHELXL-97[END_REF]. Hydrogen atoms were theoretically located. The MERCURY programme for Windows was used for structural visualisation [START_REF] Macrae | Mercury CSD 2.0-new features for the visualization and investigation of crystal structures[END_REF]. All programmes mentioned above were included in the WINGX package [START_REF] Farrugia | WinGX suite for small-molecule single-crystal crystallography[END_REF].

At the end of refinement, the reliability factors and goodness of fit were determined to have final values of R 1 = 0.0570 and wR 2 = 0.1514 and 1.001, respectively. The residual electronic density ranged from -0.474-0.489 e. Å -3 . These three values prove the good quality of this refinement. The Cambridge Crystallographic Data Centre (CCDC) now contains crystallographic data for 5 (CCDC2151896), which can be obtained via the following link:

www.ccdc.cam.ac.uk/data_request/cif.

In addition, the structure refinement showed that 5 crystallises with a monoclinic system A complex network of intermolecular hydrogen bonds, π-stacking and π-π stacking interactions ensured cohesion between molecules (see Figure 2). Indeed, a strong hydrogen-bond interaction exists between the acetamide moiety via C=O⋅⋅⋅H-N with a distance of 2.513 Å, as well as between the sulfonic acid via S=O⋅⋅⋅OH=S with a distance of 2.882 Å. Also, there exists a hydrogen-bond interaction in the methyl group of imidazole (CH 3 ) with a distance of 1.739 Å.

The distances of π-π stacking interactions around the aromatic ring are in the range of 3.387-3.392 Å, and the distance of the π-stacking interaction in the imidazole ring is 2.882 Å. The angle between the plane of the hydroxamic acid and the molecule itself is ~86.80°). The hydrogen bonds, pi-stacking and C-H….pi interactions between neighbouring molecules make crystal packing possible. Indeed, each molecule is connected to nine other neighbouring molecules with acetamide, sulfonic acid, imidazole and a methyl group (CH 3 ).

Via 5's Hirshfeld surface, which was generated using CrystalExplorer [START_REF] Turner | CrystalExplorer. Version[END_REF], we observed and analysed the behaviour of the intermolecular interactions in the crystal. The analysis of these intermolecular contacts around the asymmetric unit revealed different interactions computed by the Hirschfeld surface d norm mapping (see Figure 5). Indeed, Figure 5a

Docking Studies

The binding interactions between 5 and the active residues of sterol 14-demethylase (CYP51) were emphasised using a molecular docking study. Sterol 14-demethylase was chosen as the target enzyme for C. albicans inhibition. The docking of 5 into the active sites of the sterol 14-demethylase was performed using the PDB code 6AYB [START_REF] Debnath | CYP51 is an essential drug target for the treatment of primary amoebic meningoencephalitis (PAM)[END_REF]. Calculations were performed using the AutoDock package [START_REF] Morris | AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility[END_REF]. Validation of molecular docking was done by re-docking the original ligand ketoconazole into the active site of sterol 14-demethylase, which was reproduced with an RMSD value of 1.203 Å. Stepwise molecular docking studies were reported in our previous study [START_REF] Riadi | Synthesis, characterization, biological evaluation and molecular docking of a new quinazolinone-based derivative as a potent dual inhibitor for VEGFR-2 and EGFR tyrosine kinases[END_REF].

As seen in Table 3, 5 showed relatively potent antifungal activity against C. albicans. In an attempt to determine the binding interactions between 5 and the original docked ligand ketoconazole into the sterol 14-demethylase binding site, the binding energies and number of conventional intermolecular hydrogen bonds established in the stable complexes 5-sterol 14demethylase and ketoconazole-sterol 14-demethylase were determined. The different binding interactions of 5 and ketoconazole into the active binding site of sterol 14-demethylase are represented in Figure 6. The complexes 5-sterol 14-demethylase and ketoconazole-sterol 14-demethylase display negative binding energies, which may explain their potent antifungal activity. Furthermore, by observing the molecular docking outputs (see Table 3 and Figure 6), the original docked ligand ketoconazole demonstrated a higher (almost double) enzymatic inhibition than 5. This is mainly due to the high number of interactions established between ketoconazole and the active residues of sterol 14-demethylase. Indeed, ketoconazole interacts with 13 amino acids and the co-factor porphyrin (HEM A501), while 5 interacts with only five amino acids. Interestingly, 5 established strong π-cation and π-σ interactions with the co-factor porphyrin (see Figure 6): the π-cation interaction was established between porphyrin's ferrous ion (Fe 2+ ) and the π electrons of 5's imidazole.

Conclusion

We synthesised a new imidazole derivative 4-[2-(1-methyl-1H-imidazol-2-ylsulfanyl)acetylamino]-benzenesulfonic acid-via a simple and efficient two-step process. The synthesised derivative was evaluated for its antimicrobial efficiency against four microbial pathogens. In addition, the compound was characterised using a single crystal technique to confirm its geometric structure. Finally, the inhibitory mechanism was determined using docking studies.

Scheme 1 .

 1 Scheme 1. Synthesis of 4-[2-(1-methyl-1H-imidazol-2-ylsulfanyl)-acetylamino]-benzenesulfonic(5)

  in a P21/c space group and has the following unit-cell parameters: a = 17.069 (2) Å, b = 5.033 (3) Å, c = 17.298 (2) Å, β = 106.60° (2) and Z = 4. The unit cell contained four molecules of C 12 H 13 N 3 O 4 S 2 .

Figure 1 .

 1 Figure 1. View of asymmetric unit of the titled compound C 12 H 13 N 3 O 4 S 2

Figure 2 .

 2 Figure 2. Visualisation of the hydrogen bonds, pi -stacking and C-H….pi and interactions

Figure 5 .

 5 Figure 5. d norm mapped on the Hirschfeld surface for visualizing intermolecular interactions of C 12 H 13 N 3 O 4 S 2 .

Figure 6 .

 6 Figure 6. 3D (right) and 2D (left) closest interactions established 5-Sterol 14demethylase (up) and 5-Sterol 14-demethylase (bottom) complexes.

Table 1 .

 1 Antibacterial activity of synthesized derivative.

				Microbial Strains
			MRSA	Klebsiella pneumoniae	Aspergillus niger	Candida albicans
		MIC (μg/ml)	6.25	25	6.25	6.25
	Compound 5	MBC (μg/ml) Inhibitio	9.5	18.75	9.5	9.5
		n zone in mm (100	25	17	24	26
		μg/ml)
		IC 50 (µM)	19.11	76.45	19.11	19.11

a MIC: Minimum inhibitory concentration in μg/mL. MBC: Minimum bactericidal concentration b : No inhibition zone, no activity.

Table 2 .

 2 Docking binding energies, conventional hydrogen bonds and the number of residues closest to the docked compounds (5 and ketoconazole) into the active sites of sterol 14demethylase.

		Free binding energy (kcal/mol)	Conventional H-Bonds (HBs)	Number of closest residues to the docked ligand in the active site
	5	-7.57	0	MET A360, MET A362, LEU A358, TYR A120, ALA A39
				MET A360, PHE A216, MET A362, LEU
	ketoconazole	-13.14	0	A467, LEU A358, THR A297, CYS A430, ALA A293, PHE A292, PHE A114, VAL
				A119, ALA A289, TYR A120