A Structure-Preserving Finite Volume Scheme for a Hyperbolic Reformulation of the Navier–Stokes–Korteweg Equations

: In this paper, we present a new explicit second-order accurate structure-preserving ﬁnite volume scheme for the ﬁrst-order hyperbolic reformulation of the Navier–Stokes–Korteweg equations. The model combines the uniﬁed Godunov-Peshkov-Romenski model of continuum mechanics with a recently proposed hyperbolic reformulation of the Euler–Korteweg system. The considered PDE system includes an evolution equation for a gradient ﬁeld that is by construction endowed with a curl-free constraint. The new numerical scheme presented here relies on the use of vertex-based staggered grids and is proven to preserve the curl constraint exactly at the discrete level, up to machine precision. Besides a theoretical proof, we also show evidence of this property via a set of numerical tests, including a stationary droplet, non-condensing bubbles as well as non-stationary Ostwald ripening test cases with several bubbles. We present quantitative and qualitative comparisons of the numerical solution, both, when the new structure-preserving discretization is applied and when it is not. In particular for under-resolved simulations on coarse grids we show that some numerical solutions tend to blow up when the curl-free constraint is not respected.


Introduction
Many models of mechanics and physics in general are described by a set of evolution equations, in which some state variables are also subject to stationary differential constraints (involutions), restricting the set of admissible solutions.Some of the most notable examples are constraints on the divergence or the curl of a given field, which are encountered for example when dealing with Maxwell's equations of electrodynamics or with magnetohydrodynamics, which require the magnetic field to be divergence-free for all times.In many cases, the constraints are supplied as an extra equation to the system, i.e., a given field obeys a time-evolution equation, but is also bound independently by an additional property.This is the case, for example, in the incompressible Navier-Stokes equations.The time evolution of the velocity field is governed by the momentum conservation equation, but has to also satisfy a supplementary divergence-free constraint due to the incompressibility assumption.We call the constraint supplementary because it is not a consequence of the evolution equations, but an additional property.In the incompressible Navier-Stokes equations the divergence-free condition is actually the pressure equation, see, e.g., [1][2][3]: even if the velocity is divergence-free initially, there is no guarantee that it remains so for later times by means of its evolution equation alone.
Yet, in some other examples, it occurs that the constraint is intrinsic to the mathematical structure of the evolution equation, that is exact solutions verify it by nature.For instance, a potential flow has an intrinsic curl-free constraint on the velocity, since the latter is expressed as the gradient of a scalar field.Thus, at the theoretical level, the solutions verify the curl-free constraint by nature.In Maxwell's equations of electromagnetism, if the magnetic field is initially divergence-free, its time evolution equation ensures it remains so for later times.In such cases, we are talking about involution constraints, meaning they are direct consequences of the structure of the system and they do not need to be supplied independently.Several systems of different scopes and applications fall into this category, like the general symmetric hyperbolic and thermodynamically compatible (SHTC) class of PDEs forwarded by Godunov and Romenski and studied in [4][5][6][7][8][9], the multi-phase flow model of Romenski et al. [10] or first-order reductions of the Einstein field equations like those studied in [11][12][13].
While it might seem somewhat comforting that involution constraints are automatically satisfied by exact solutions of the PDE, provided that compatible initial conditions are given, this is a completely different story at the numerical level.Particular attention has to be directed to such constraints, even if they are an intrinsic property of the system at the continuous level.Indeed, the accumulation of discretization errors over the course of a numerical simulation can lead to the violation of curl-free or divergence-free constraints.This leads to numerical solutions that are neither compatible with the physics of the studied system, nor with its mathematical structure, see for example [14][15][16].In some cases, the buildup of such errors can lead to unstable numerical simulations and to finite-time blowups [17,18].The aforementioned issues have stimulated a lot of research in numerical analysis for hyperbolic PDE in recent years in order to develop better techniques and specific numerical methods that are able to address these notorious problems and provide numerical results that are compatible with the mathematical structure of the equations and with the underlying physics.Such techniques include for instance the so-called constrained transport methods [19][20][21][22][23][24][25][26], the Generalized-Lagrangian multiplier divergence-cleaning approach of Munz et al. [27][28][29][30], which was later also extended to hyperbolic curl cleaning in [13,18,31,32], as well as mimetic discretizations and reconstructions [16,[33][34][35][36][37].
The hyperbolic reformulation of the NSK system incorporates a hyperbolic relaxation model of the EK Equations [71][72][73] into the unified Godunov-Peshkov-Romenski (GPR) model of continuum mechanics [4][5][6][7].It admits a curl-involution constraint, which will be the main obstacle to address.The current paper is a natural sequel to the previous work of the authors [32], where a hyperbolic Generalized Lagrangian Multiplier (GLM) curl cleaning technique was employed [13,18,31].In this paper, we design a new second-order accurate finite volume scheme, based on a carefully chosen vertex-based staggered grid, allowing to define a discrete gradient operator for which the discrete curl is identically null at the machine level.The scheme inherits most of its foundations from the works [17,74].Similar techniques were developed previously for hyperbolic systems with divergence constraints, see for example [16,23,33,[75][76][77].Advantages of this type of method include for instance the fact that the PDE system in itself remains unchanged (unlike in GLM cleaning techniques for example) and that the discrete curl-errors are zero up to machine accuracy at the numerical level.This comes at the price of additional efforts in the implementation, as the method generally requires a well-chosen spatial discretization and a careful (staggered) positioning of the relevant variables.
The rest of this paper is organized as follows.In Section 2, we briefly recall the hyperbolic reformulation of the NSK model; we detail its structure and mathematical properties.In Section 3, we lay out the details of our new structure-preserving finite volume scheme.In particular, we describe the structure of the staggered computational grid, the splitting of the system, how each subsystem is evolved in time and how the compatible discrete differential operators are defined.Finally, in Section 4, we show how the presented numerical method performs on a set of different test cases in one and two space dimensions.In particular, we test the scheme on a standard 1D Ostwald ripening test case.We provide a convergence table for a stationary 2D bubble solution.Then, an extensive study of the curl errors is performed for a non-stationary non-condensing bubble test case, where we compare curl errors with and without staggering, even for alternative discrete curl operators.The last test-case is the standard 2D Ostwald ripening for which we show that our simulations perform well when strong topological changes in the solutions are present.The paper is rounded-off with some conclusions and an outlook to future work in Section 5.

Hyperbolic Reformulation of the Navier-Stokes-Korteweg Equations
The NSK equations can be seen as a dispersive modification of the Navier-Stokes equations, based on the Van der Waals-Korteweg gradient theory [78].It serves as a general model for compressible viscous and dispersive fluid flows.In our case, we will be considering an isothermal barotropic version of the equations, similar to [50].In its general form and under these considerations, the NSK system reads as Here and in the rest of this paper, repeated index summation is implied.x ∈ R d is the vector of Cartesian coordinates in d space dimensions and t is the time.In terms of notations, any system variable that is introduced in the following manuscript is assumed to depend on space and time.The notation (x, t) will be omitted to ease notation.In system (1), ρ is the density field and u = (u 1 , • • • , u d ) is the velocity field.P = P(ρ) is the pressure field, given as a function of the Helmholtz free energy by P = ρW (ρ) − W(ρ).δ ik is the Kronecker symbol and also denotes the identity matrix.The tensors σ = [σ ik ] and κ = [κ ik ] are the viscous Navier-Stokes stress tensor under Stokes' hypothesis and the dispersive Korteweg stress tensor, respectively.Their expressions are given by and In these expressions, µ is the dynamic viscosity and γ is the capillarity coefficient, both assumed constant.The dissipationless part of the momentum Equation (1b) can be obtained as an Euler-Lagrange equation, derived from a variational principle applied to the following Lagrangian under the differential constraint of total mass conservation

The Van der Waals Equation of State
The NSK system is commonly used to describe viscous compressible multiphase flows via a diffuse interface approach.It is required in this context to have an equation of state allowing multiple admissible thermodynamic equilibria, such as the double-well potential introduced by Van der Waals [79].In this case, the Helmholtz free energy admits two distinct local minima, each corresponding to a different state of the fluid, typically liquid and vapor phases.The free energy and the corresponding pressure are given in their dimensionless form by where R is the ideal gas constant, T is the temperature, assumed constant, c v is the heat capacity at constant volume and a and b are constant parameters.In the rest of this paper, the following values will be taken for which a representative graph of the pressure P(ρ) is given in Figure 1.V and ρ M L are the Maxwell states for the vapor and liquid phases, respectively.ρ * V and ρ * L are local extrema of the pressure and delimit the non-convexity region.

First-Order Hyperbolic Reformulation of the NSK Equations
The system considered here for numerical simulation is the first-order hyperbolic relaxation model approximating the NSK equations presented in [32] and whose structure and properties are briefly recalled in this section.The model can be seen as the combination of the unified hyperbolic and thermodynamically compatible Godunov-Peshkov-Romenski (GPR) model of continuum mechanics [4][5][6][7], with an augmented Lagrangian approach [71,72], which allows to cast dispersive systems of the Euler-Korteweg type into first-order hyperbolic PDEs with stiff relaxation terms.The full system of equations, is given by Additional degrees of freedom are introduced here with respect to the original system (1) and whose meanings are recalled hereafter.A = [A ik ] is the so-called distortion field.It describes the deformation of material elements in the GPR model framework.In the absence of source terms in Equation (5c), A would be equivalent to the inverse deformation gradient.η, w and p are auxiliary variables introduced in the setting of the augmented Lagrangian method.η can be seen as the proxy for ρ as the order parameter whose gradient is used in the Korteweg theory.The gradient and material derivative of η are promoted to new independent variables, respectively, p and w and whose evolution equations are established from their definitions [71,72].Naturally, their evolution equations need also to be supplemented with the appropriate initial conditions w(x, t = 0) = η(x, t = 0) and p(x, t = 0) = ∇η(x, t = 0).Being the gradient of a scalar field, p inherently admits a curl involution constraint ∇ × p = 0.One can trivially verify it by applying the curl operator on Equation (5e) and recall the corresponding initial condition to obtain leading to ∇ × p = 0 for all times.One could remark that, for the same reasons, the non-conservative part of Equation (5e) is theoretically zero as well.However it must be kept for structural purposes and in particular to have Galilean invariance.The total momentum conservation Equation (5b) has a seemingly similar structure as its counterpart in the original NSK model (1b).However, both equations are fundamentally different since the viscous and dispersive stress tensors in (5b), denoted by Σ = [Σ ik ] and K = [K ik ], respectively, do not involve any high-order derivative of the thermodynamic degrees of freedom.Indeed, their expressions are given by where the tensor G = [G ik ] is the so-called metric tensor obtained from the distortion field as G mm δ ik is its deviatoric (trace-free) part.Furthermore, c s is a relaxation speed for the GPR part of the model and can be seen as the characteristic speed for shear waves, as evidenced by the eigenstructure of the system.P is an additional pressure coming from the relaxation of the dispersive terms.The expression for the source term S A = [S A ik ] that appears in Equation (5c) is given by where τ is the characteristic time of the strain relaxation.Finally, one should note that the dissipationless part of system (5), i.e for S = 0, can be derived from Hamilton's principle of stationary action, applied to the Lagrangian under the constraints (5a) and (5c).( 5b) and (5f) are then the Euler-Lagrange equations, (5e) and (5d) are closure relations that are immediate consequences of the definitions.The total energy density E = ρE of this system, with E the specific total energy, is obtained immediately from the Lagrangian and is given by One can derive an additional conservation law for it as a consequence of system ( 5) by summing its equations, each multiplied by the corresponding entropic variable to obtain

Hyperbolicity
The hyperbolicity of system ( 5) has been addressed in [32], albeit in a slightly different case, where an additional evolution equation for the so-called curl-cleaning field is also present.In some similar models in the literature [18,31], GLM curl-cleaning is not only a means to control discrete curl errors, but it is also necessary in order for the equations to be strongly hyperbolic in multiple dimensions of space.We show hereafter that this is not the case for system (5).Since the latter is invariant by rotations, it is enough to consider the onedimensional case, i.e all the variables are only dependent on the coordinate x and the time t and all vector components are kept.This allows to simplify the algebra.We consider the vector of primitive variables V = (ρ, u 1 , u 2 , u 3 , η, w, p 1 , p 2 , p 3 , A 11 , A 12 , A 13 , A 21 , A 22 , A 23 , A 31 , A 32 , A 33 ) T and we linearize the system around the reference state A = I and p = (p 1 , 0, 0).This allows to reduce the system in its quasilinear form which reads where S are the sources and M is the quasilinear matrix given by Direct computations show that the eigenvalues are given by where Z 1 and Z 2 are auxiliary quantities introduced to lighten the expressions, in the same fashion as in [32] and whose definitions are recalled here for completeness: Trivial algebra shows that the eigenvalues λ 1−18 are real for a judicious, but not very restrictive choice of the relaxation parameters.In fact, in the region of non-convexity of the Helmholtz free energy, even though a 2 0 is negative, one can choose for example α sufficiently small so that Z 1 remains positive.The reader is referred to [32] for further details.In this case, the eigenvectors associated with λ 1−18 form a basis of R 18 , even in the regions where the Helmholtz free energy W(ρ) is not convex.Therefore, under reasonable assumptions, the system of equations given by ( 5) is strongly hyperbolic.

Setting and Notations
In order to numerically solve the system of Equations ( 5), we propose here a discretization based on a staggered grid, that is compatible with the intrinsic curl-free constraints present in the system.Such an exactly curl-free discretization was recently introduced in similar contexts but for different PDE systems in [17,74] and we recall here the necessary details of the numerical scheme.Before proceeding any further, it seems necessary to briefly define the setting as well as the corresponding notations that will be used throughout this section.In an attempt to keep the text clear and concise, we shall provide a thorough description in two-dimensions of space, from which the extrapolation to a third-dimension is rather straightforward.In all what follows, the coordinates shall be denoted by x and y instead of x 1 and x 2 for less cumbersome notations.The computational domain will be denoted by Ω c and spans over [−L x /2, L x /2] × [−L y /2, L y /2], where L x and L y are prescribed in advance.It is discretized along the x-axis and the y-axis, forming a uniform Cartesian grid of N x × N y cells, whose barycenter coordinates are denoted by (x p , y q ) and whose expressions are given in a classic form by where the mesh spacings are given by ∆x = L x /N x and ∆y = L y /N y , respectively.An arbitrary grid cell thus spans the subdomain ].In order to avoid confusion between tensor indices and discretization indices, throughout this paper we will use the subscripts i, j, k, l, m for tensor indices, the indexes p, q for the discretization indices in space and the index n for the time discretization.
In the setting of the exactly curl-free method developed here, staggered grid points need to be defined, which are required to compute discrete differential operators that are compatible with the system's structure.The staggered grid points are defined in the vertices of the grid cells.For an arbitrary cell of center (x p , y q ) these are the points of coordinates (x p± 1 2 , y q± 1 2 ).A representation of the grid points is given in Figure 2. Last, the time is discretized and the discrete time t n = n∆t at time step n will be denoted as usually by t n .For every variable φ, the value at coordinates (x p , y q ) at time t n shall be denoted as φ p,q,n .The time superscript will be omitted in most parts, and will be explicitly written only when necessary.
Schematic of the computational domain featuring the grid points and the staggered dual grid points.Red squares are barycenters and blue circles are the vertexes of the computational cells.

Flux Splitting
In order to solve system of Equations ( 5) numerically, we cast it into the following split form where Q is the vector of conserved variables.F b (Q) and S b (Q) are the fluxes and the sources, respectively, which will be discretized on the barycenter-based grid points.The terms F v (Q), ∇G v (Q) and B v (Q) • ∇Q are, respectively, the flux terms, the gradient terms and the non-conservative products that will be discretized on the cell vertices.S v (Q) are the source terms coming for the GPR model and which will be also evolved on the vertices.If we denote by F k the k th column of a flux tensor F, then all the previously described terms are given by , and The fluxes F k b contain the convective part of the system together with the pressure terms, both coming from the NSK model itself, as well as from the augmented Lagrangian relaxation procedure.The fluxes F k v only contain the stress tensors, from which the pressure P is subtracted.G v contains the terms whose gradient will be computed by using a compatible discrete gradient operator, to be defined in the next section.Note that in our case, there is only one curl-involution constraining the vector p.In general, if the relaxation source term S v is zero, meaning that dispersive solid mechanics are considered, a curl-free constraint would also affect each of the lines of the tensor A, which are evolved in time.Therefore, regardless of the considered applications, one loses nothing by discretizing the equation on A accordingly, so that a future extension to such applications would be straightforward.Under the previous notations, the system of Equations ( 9) will be discretized explicitly, only for the stiff relaxation source term, S v which will be integrated using an implicit Euler scheme.The fluxes F b and F v as well the sources S b will be discretized using an unlimited second-order accurate MUSCL-Hancock-type finite volume scheme [80], based on a Rusanov-type approximate Riemann solver.The terms ∇G v , B v • ∇Q and S v will be discretized separately, in a structure-preserving way, hence ensuring that the curl-free constraint on p is respected exactly at the discrete level.Details of this dicretization are given in the next section.

Compatible Discrete Operators and Curl-Free Discretization
The main foundation of the structure-preserving feature of the scheme presented in this paper relies on the proper definitions of the discrete differential operators, making use of the carefully chosen grid points.First, it is of utmost importance to recall that the discrete variables, written here in primitive form, are stored either on the main grid points (x p , y q ) or on the dual grid points (x p± 1 2 , y q± 1 2 ) according to the following distribution (see also Figure 2): Main grid : ρ p,q , u p,q i , η p,q , w p,q , Dual grid : Some update formulas, as will be seen further down in this section, will require either a value of the dual variables in the cell centers or vice-verse.For that, and for a generic variable denoted by ψ, we shall make use of the following averaging formulas from the corners to the cell barycenters and from the barycenters to the corners.The formulas are, respectively, given by the simple arithmetic averages as follows: Now, it is necessary to introduce the definitions of the discrete gradient and discrete curl operators, which are at the core of our structure-preserving numerical scheme.For a given scalar field stored in the centers of the control volumes and denoted by φ p,q , one can compute a natural discrete gradient ∇ h φ in the corners and whose components in the x and y directions are computed from the surrounding cell centers as follows where Based on this corner gradient, one can now define a compatible discrete curl operator ∇ h × ∇ h φ computed in the cell centers, from the surrounding values of the discrete gradient and whose component in the z direction is given as follows It is easy to see that under this definition, substituting the expressions of the discrete gradient ( 11) into ( 12) results into the discrete identity for all cells in the computational domain, see also [74].Clearly, this means that any gradient field defined by means of the discrete operator ( 11) is exactly curl-free for the discrete curl operator (12).It now remains to set up a compatible time-evolution that maintains this exactly curl-free property.Thus, for each component k of the gradient field p, the time update equation is given by The last term in particular is a compatible discretization of the term u × (∇ × p) and straightforward computations [74] show that if ∇ h × p p+ 1 2 ,q+ 1   2 ,n = 0, then also ∇ h × p p+ 1 2 ,q+ 1   2 ,n+1 = 0 holds.A similar update formula will be considered for the distortion field A. It reads where the last term on right hand side corresponds to the potentially stiff relaxation source term, which is discretized implicitly in time.In general, as long as this relaxation source term is present in the evolution equation, there is no curl involution affecting the rows of A.
It is however the case when τ → ∞.Although this will never be the case in the scope of this paper, the time evolution for A was also discretized in the same compatible way as the time evolution of p, in an attempt to keep the scheme as general as possible and extensible to different regimes, different from the ones considered here.For both update formulas presented above, the scheme relies solely on a central discretization.This means that a compatible artificial viscosity must be introduced, in the same fashion as in [17,74].At the continuous level, we recall that the vector Laplacian of p reads as While we have already defined a discrete analogue to the gradient and curl operators, we have not defined a discrete divergence operator yet.In order to have second-order accurate numerical dissipation terms, one should consider a piece-wise linear reconstruction of p in the cell centers as the stencil for the discrete divergence operator, computed in the corners.First, for every vertex of coordinates (x p , y q ), one needs to compute vertexextrapolated values of p, denoted by p ne , p nw , p se and p sw , where the letters in subscript refer unambiguously to the four cardinal directions.The discrete divergence operator at the vertex (x p+ 1 2 , y q+ 1 2 ) is then given by 2∆y .
Based on this definition, a discrete analogue of the vector Laplacian ( 16), at the vertex Next, one multiplies this Laplacian by the mesh size h = max(∆x, ∆y) and by an appropriate representative velocity c * , then inserts it to into the update formula (14), to obtain where ε ijk denotes the classical Levi-Civita tensor.In a similar manner, the same procedure is adapted to include numerical viscosity into the evolution of the A at the numerical level, which now reads as In practice, we take the value c * = c s .One last remark concerns the algebraic constraint det(A) = ρ/ρ 0 , which is a consequence of the time evolution equation of the distortion field A, see [6,7].It is possible to enforce this constraint simply by uniformly rescaling the object A at every time step according to the update formula with Ã the rescaled distortion field and A the distortion field before rescaling.This completes the description of the discretization of the terms ∇G v (Q), B v (Q) • ∇Q and S v (Q).

Second Order MUSCL-Hancock-Type Scheme for the Remaining Terms
The subsystem which contains the remaining terms F b , F v and S b , i.e., is discretized explicitly in time according to a conservative but unlimited MUSCL-Hancocktype scheme, see [80].The time update formula reads as usual, i.e., where the columns of the numerical flux (F 1 ) and (F 2 ) are defined in the corresponding cell edges as follows: (F 1 ) .
In these expressions, the terms s x max and s y max are the maximum signal speeds in their respective direction, computed from the local maximum eigenvalues λ max = max k∈ [1,...,18] λ k of the quasilinear matrix M. Q ± are the boundary extrapolated values, obtained from a piecewise linear reconstruction from the neighboring cell centers according to the formulas and where ∂ x Q p,q,n and ∂ y Q p,q,n are the unlimited slopes in the x and y directions, respectively, given by and the approximation of the time derivative ∂ t Q p,q,n is computed from the PDE as follows This completes the description of the unlimited second-order MUSCL-Hancock-type finite volume scheme which is used for the discretization of subsystem (19).

Numerical Results
In all what follows, whenever a configuration is described for the HNSK (Hyperbolic reformulation of the NSK) system, we shall provide the initial conditions for the density ρ and the velocity field u, as would have been the case for the original NSK system.For all the remaining variables, the corresponding initial values are computed accordingly, as a function of ρ and u as follows: A(x, 0) = I, η(x, 0) = ρ(x, 0), p(x, 0) = ∇ρ(x, 0), w(x, 0) = −ρdiv(u(x, 0)).
In particular, at the numerical level, we recall that it is necessary to initialize p using the compatible gradient ∇ h p given in (11).Concerning boundary conditions, and unless specified otherwise, periodic boundaries are considered.We recall that in the scope of the NSK model, only one fluid is under consideration and bulk phases are only identified by the values of the order parameter (ρ for NSK and η for HNSK).Although no interface per se exists, we might refer to the continuous transition between the liquid and vapor states densities as a diffuse interface, or just interface for shortness.In the scope of the GPR model, we remind that unlike the original NSK equations, the model does not depend explicitly on the dynamic viscosity µ.However, in the stiff relaxation limit, it can be shown via formal asymptotic analysis that at leading order, the equivalent viscosity is given by µ = τρ 0 c 2 s /6, see [6] for details.Therefore, in what follows and for each test case, we shall report the equivalent µ and the shear velocity c s .The reference density ρ 0 shall be always taken as ρ 0 = 1.The time-step is constrained by a usual Courant-Friedrichs-Lewy condition of the form In what follows, we take ∆t = CFL ∆x/λ M where λ M = max Ω λ max and CFL < 1 is a constant coefficient whose value is reported in every test case.

One-Dimensional Ostwald Ripening
We show a first test case in order to validate the numerical scheme in one dimension of space.Ostwald ripening occurs when a set of different-sized bubbles indirectly merge, without actually moving one towards the other.To summarize the process, smaller bubbles disappear and free their masses into the medium which is absorbed by bigger bubbles.The process continues until one large bubble remains.We consider an initial configuration, adapted from [32,50,[53][54][55] described by the following generic initial data where ρ l represents the initial liquid phase density and ρ v is the initial vapor phase density.n b is the number of bubbles initially present.For each bubble of index k, x k is the position of the bubble center and r k is the approximate radius, so that each bubble spans approximately over the region . is a parameter controlling the initial thickness of the 'interface', that is the steepness of the jump from ρ l to ρ v and vice versa.Note that this affects only the transient states as the equilibrium interface width only depends on the capillarity coefficient γ.For this simulation, we take the same parameters that were used in [32] and which are recalled here for completeness: The computational domain and the corresponding discretization parameters are The numerical results are plotted in Figure 3, along the reference numerical solution obtained in [32].The comparison shows an excellent agreement between both solutions.

2D Stationary Droplet
It is possible to compute numerically a quasi-exact stationary solution for the HNSK system in two space dimensions, either by solving an initial value problem or a boundary value problem.We shall use the latter in this case and consider the reference solution computed in [32] using the COLNEW BVP solver, see [81][82][83].We recall that this solution is obtained from solving the following BVP and is provided for reference in Figure 4.The structure of the solution involves a fast transition separating the liquid and vapor phases and which constitutes the diffuse interface.The peak observed in the gradient field corresponds naturally to the inflection point in the density profile.In order to simulate numerically this solution we take the following parameters over different-sized grids.The final time is set to t = 0.025 and we take CFL = 0.9.The numerical result is then compared with the initially provided exact solution.A convergence table for this solution is provided in Table 1 for a set of variables.The obtained computational results show that our new structure-preserving scheme achieves second-order of accuracy, as expected.Table 1.Numerical convergence table for the 2D stationary bubble reference solution obtained from solving (21).N x and N y denote the number of elements in the x and y directions.The L 2 norm of the errors for a set of representative variables with respect to the exact solution, are reported in the upper part and the corresponding convergence order O 2 is reported in the lower one.Parameters used in the simulations are γ = 10 −3 , µ = 10 −2 , c s = 1, α = 10 −2 and β = 10 −3 .Final time is t = 0.025.

Non-Condensing Bubble
In order to study quantitatively the numerical curl errors, we consider here a noncondensing 2D bubble.The initial configuration is taken in the form This describes a diffuse bubble centered at the origin with an approximate radius r 0 .The domain is taken Ω c = [−0.25,0.25] × [−0.25, 0.25] discretized over 1024 × 1024 cells.The remaining physical parameters are taken as follows: The CFL number is set to CFL = 0.9 and all boundary conditions are periodic.The initial condition is not at equilibrium and the dynamics involve radial fluctuations around the diffuse interface region until a steady state is reached.The time evolution of curl errors is plotted in Figure 5 in two different configurations: one with our new structure-preserving and exactly curl-free discretization, and one without, that is all the variables are simply evolved in the barycenters of the control volumes, using the MUSCL-Hancock-type scheme.
One can clearly see that the discrete curl errors are kept constantly around machine precision when the new compatible and exactly curl-free scheme is used, whereas the curl errors are about eight orders of magnitude higher and growing in time, otherwise.This is a natural consequence of the choice of the grid points used for the computations that ensure such a result for the definition of the proposed discrete curl operator.However, a discrete curl does not have a unique definition.Hence, the following question arises: does the staggered scheme still perform well regarding curl errors when an alternative arbitrary discrete curl definition is employed for a posteriori measuring the curl errors?Thus, let us consider a discrete curl computed by centered finite differences, using barycenter-based values b p,q , interpolated a posteriori from the vertices accordingly to (12).This alternative curl in a cell of coordinates x p , y q , denoted by ∇ h alt × p p,q , reads as We plot hereafter the time evolution of the L 1 norm of the discrete curl errors according to the alternative definition (24) and the compatible discrete curl operator (12), for the same parameters reported in (23) and for different mesh sizes.Note that the alternative curl (24) is only used for the postprocessing of the data and for measuring the curl errors.It is not used inside the numerical scheme itself.The obtained computational results are reported in Figure 6.The results clearly show that for the new structure-preserving finite volume scheme proposed in this paper which is based on the discrete curl ∇ h × p, for which the method is proven to be exactly curl-free, the curl errors are stable and remain only slightly above machine accuracy, due to round-off errors and which obviously increase with mesh refinement.However, when a simply barycenter-based MUSCL-Hancock method is used, which is not exactly curl-free, the numerical solution becomes unstable and blows up in finite time.The time evolution is thus reported up to the time when this blowup occurs.Similar blowups have been reported previously for hyperbolic PDE with curl involutions in [17,18].One can observe a significant increase in the curl errors before that happens.On the other hand, for the postprocessing with ∇ h alt , it is interesting to observe that even when one changes the definition of the discrete curl operator to an arbitrary one which the method is not meant to preserve exactly, the structure-preserving method still performs much better concerning the curl errors compared to the standard MUSCL-Hancock scheme.The curl errors of the structure-preserving scheme are always lower and essentially do not display any increasing trend and remain bounded over time.In order to investigate the long-time behavior further, which is very important for structure-preserving schemes, we provide in Figure 7 the discrete-curl evolution for the same simulation, but over a longer time for two different meshes.For this simulation, all the model parameters are kept the same as in Figure 6.The figure shows that the discrete curl-errors remain around the machine zero range, even for long times and millions of time steps.The slight increase is of course due to the buildup of round-off errors, which do not necessarily cancel, but which accumulate in time and which are thus naturally greater for more refined meshes.At this point, it seems also interesting to investigate whether the use of an exactly curl-free scheme has any notable repercussions on the overall behavior of the numerical solution.Thus, we plot in Figure 8 a comparison of the density profile and of the gradient field component p 1 for the same non-condensing bubble test, with and without a structurepreserving discretization, at an intermediate time given by t = 2.The results show that a deformation of the bubble can be observed around the interface region and a glitch can be observed in the center of the bubble when no curl-preservation is in effect.These deformations grow over time and are seemingly what causes the failure of the numerical simulation at later times.To illustrate this further, we also plot the norm of the discrete curl field |∇ h × p| for both simulations in Figure 9.A quick comparison with Figure 8 shows that the curl errors are mainly concentrated around the same locations where the deformations are observed and especially around the middle glitch.

Two-Dimensional Ostwald Ripening
A standard test case for the NSK system is the 2D Ostwald ripening, see, e.g., [32,50,[53][54][55].The principle and the setting is similar to what is proposed in the one-dimensional counterpart and we consider an initial condition expressed in Cartesian coordinates by the general expression where the center of the bubble k is located at (x k , y k ), c k is the corresponding local radial coordinate c k = (x − x k ) 2 + (y − y k ) 2 and r k is the approximate radius of the bubble k.
We take here a configuration with n b = 10 bubbles, whose initial positions and approximate radii are summarized in Table 2. Since the dynamics of the solution involve the vanishing of the smaller bubbles and the expansion of the bigger ones, the time evolution involves important topological changes, implying significant dynamics in the order parameters and in the gradient field p.We consider a computational domain  Figure 11 shows that the gradient field evolved independently, is also qualitatively in agreement with the density field.In order to investigate the results more quantitatively, in particular to check the compatibility of ρ with the order parameter η and if the gradient field p is in line with the gradient of η, we plot the variables of interest at t = 5 on the lines defined by x = −0.05 and y = 0.1, both of which cross the final bubble around its center.We compare both ρ and η and p with ∇η computed numerically with finite differences.The comparisons are shown in Figure 12.The results show that the order parameter remains an excellent approximation to the density field.Both components of the independent gradient field p are also in agreement with the true gradient ∇η in both directions, even after significant deformations occur in the solution.comparisons are shown in Figure 12.The results show that the order parameter remains an excellent approximation to the density field.Both components of the independent gradient field p are also in agreement with the true gradient ∇η in both directions, even after significant deformations occur in the solution.

Conclusions
In this paper, we have presented a new second-order accurate exactly curl-free scheme on staggered grids to solve a first-order hyperbolic reformulation of the Navier-Stokes-Korteweg system presented in [32] while also preserving its natural curl involution constraints exactly at the discrete level.The tests in two dimensions of space show that the discrete curl errors remain around machine precision and are stable in time.We also showed that in the absence of an exactly curl-free discretization, erroneous behaviors can be observed in the solution, leading to finite-time blowups, especially for under-resolved meshes.
Future research will concern the extension to three space dimensions, which is out of scope of the present paper.Another future direction of research would be to improve the current scheme by a semi-implicit time discretization, similar to the one forwarded in [17,74,77].This requires first to find an appropriate splitting of the hyperbolic system, allowing to separate the fast relaxation characteristics so that a semi-implicit discretization would improve the overall performance and reduce the time-stepping constraints.The development of arbitrary high-order structure preserving methods is also considered on the longer term, following the ideas presented in [84,85].Future work will also concern the extension of the presented exactly curl-free discretization to a thermodynamically compatible one, similar to the ideas on HTC schemes recently presented in [86][87][88].The framework of HTC schemes allows the construction of provably nonlinearly stable schemes in the energy norm and to the best of our knowledge, there are no exactly curl-free HTC schemes available yet.To date, in HTC schemes the involution constraints have been taken into account via a hyperbolic and thermodynamically compatible GLM cleaning approach, see, e.g., [89].

Figure 1 .
Figure 1.The Van der Waals pressure as a function of ρ for the values given in (4).ρ M V and ρ M L are the Maxwell states for the vapor and liquid phases, respectively.ρ * V and ρ * L are local extrema of the pressure and delimit the non-convexity region.

Figure 4 .
Figure 4. Boundary value problem reference solution for the HNSK system.The left figure shows the radial profiles of ρ and η while the right one shows the radial profile of the gradient field, for α = 10 −2 and γ = 10 −3 .

Figure 7 .
Figure 7.Long-term evolution of the discrete curl errors, for the non-condensing bubble test case.The left figure corresponds to a simulation with a 512 × 512 grid, while the right figure is for a 128 × 128 grid.The model parameters taken here are the same as in Figure6.The upper horizontal axis indicates the number of time steps, while the lower axis reports the time t.

2 Figure 8 .Figure 9 .
Figure 8.Comparison of the overall shape of the density field ρ (top) and the gradient field component p 1 (bottom) with both a curl-free discretization (left) and without (right).Results are shown for t = 2 on the 512 × 512 grid.

Figure 11 .
Figure 11.Numerical result for the magnitude of the gradient field p for the same simulation as in Figure 10.

Figure 11 . 10 . 2 Figure 12 .
Figure 11.Numerical result for the magnitude of the gradient field p for the same simulation as in Figure 10.

Figure 12 .
Figure 12.Comparison of the density profile with the order parameter η (top) and of the gradient field components p 1 and p 2 with ∂ h x (η) and ∂ h y (η), respectively, computed with finite differences.The comparison is made at t = 5, on two extracted lines from the results of Figures 10 and 11, defined by y = 0.1 (left) and x = −0.05(right).

Table 2 .
(25)ary of the coordinates and radii of the bubbles in the initial condition(25).