A Multiscale Method for Two-Component, Two-Phase Flow with a Neural Network Surrogate

Understanding the dynamics of phase boundaries in fluids requires quantitative knowledge about the microscale processes at the interface. We consider the sharp-interface motion of compressible two-component flow, and propose a heterogeneous multiscale method (HMM) to describe the flow fields accurately. The multiscale approach combines a hyperbolic system of balance laws on the continuum scale with molecular-dynamics simulations on the microscale level. Notably, the multiscale approach is necessary to compute the interface dynamics because there is -- at present -- no closed continuum-scale model. The basic HMM relies on a moving-mesh finite-volume method, and has been introduced recently for compressible one-component flow with phase transitions in [Magiera and Rohde, JCP. 469 (2022)]. To overcome the numerical complexity of the molecular-dynamics microscale model a deep neural network is employed as an efficient surrogate model. The entire approach is finally applied to simulate droplet dynamics for argon-methane mixtures in several space-dimensions. Up to our knowledge such compressible two-phase dynamics accounting for microscale phase-change transfer rates have not yet been computed.


Introduction
In this work we consider the dynamics of a compressible two-component fluid which allows for phase changes between a liquid and a vapor phase state.Keeping temperature constant, we aim to determine continuum-scale quantities like component-wise densities and momenta on a spatial scale where all phase boundaries are explicitly represented, i.e., without further averaging like in e.g.Baer-Nunziato modelling [4,37,43].This choice implies that the phase -liquid or vaporis uniquely determined by the values of the component densities.Having made this decision, we focus on a sharp-interface approach, and refer to [3] and just recently [24] for diffuse-interface models.We start from an extension of the isothermal Euler equations with frictional forces for one-phase flows, as introduced in [6], to the two-phase case, see (1) below.We will focus mainly on the two-component system comprised of the noble gases argon and methane which allows a relatively simple thermodynamical description.These evolution equations are supposed to hold in the liquid-phase domain  − () and the vapor-phase domain  + (), where they form a set of hyperbolic balance laws.The two bulk domains are separated by a co-dimension-1 manifold, the interface  = ().Thus, the entire fluid domain  ⊂ R  ,  ∈ {1, 2, 3}, at any time  ∈ [0,  end ],   end > 0, is partitioned into  =  + () ∪ () ∪  + (), see Figure 1 for illustration.The basic mathematical model is provided in Section 2. To determine the motion of the interface and to ensure the well-posedness of the resulting free boundary value problem one has to prescribe transmission conditions across the interface ().Obviously, the Rankine-Hugoniot conditions have to be imposed; but it is well-known, even in the one-component case, that additional conditions have to be prescribed at the interface [8,17,21,36].This can be achieved by imposing further algebraic relations that prescribe the entropy dissipation rate across the interface [40].However, determining the explicit form of this relation remains a largely unsolved issue for complex flow regimes, including in particular multi-component flow.To overcome this fundamental modelling problem on the continuum scale and to avoid using ad-hoc closures, one can look at smaller spatial scales using a molecular-dynamics (MD) approach.On this ab-initio level the modelling of a multi-component fluid becomes accessible since each molecule of a component has the same physical properties regardless of the fluid phase it is part of (see Section 3).However, MD simulations can only be run for very limited spatial and temporal scales.We utilize a tailored MD algorithm to solve MD Riemann problems across the (discrete) interface only, whereas the dynamics in the bulk phase domains is approximated by a standard finite-volume method for the first-order two-component system (2).With this hybrid ansatz we propose a heterogeneous multiscale model (HMM, in the sense of [10]) that accounts for physically accurate interfacial mass and momentum transfer, as well as the motion of the phase boundary by an MD interface solver, see Section 4. A critical issue in this context is the coupling of the different models, i.e., from MD systems with pairwise interactions to continuum-scale Euler-type equations and vice versa.We rely on classical, statistical averaging to infer continuum-scale quantities from microscale particle simulations.In that way, we transfer the microscale behavior at fluid interfaces to the continuum scale via a data-based approach.Concerning the numerical discretization on the continuum scale we employ a moving-mesh finite-volume method which enables us to include the MD dynamics directly at the phase boundary, while resolving the sharp interface directly within the mesh.The corresponding moving mesh-algorithms is presented in [7,32], see [1] for the corresponding open-source code.Finally, to overcome the computational costs of expensive MD simulations, we apply a surrogate model for the MD interface solver.This surrogate model is based on the constraint-aware neural networks that have been developed in [30], which come up with the same computational efficiency, but ensure mass conservation of the numerical discretization, in particular across the interface.
In Section 5, we present numerical results for the HMM for compressible two-component flow with phase transitions.We focus on a mixture of the noble gases argon and methane using their exact physical properties.In particular this includes a specific equation of state that is used in the bulk domains.The results comprise simulations in one, two, and three space-dimensions and settle, for multiple space-dimensions, around the impingement of pressure waves on droplets inducing oscillation and condensation processes.We summarize our findings in the concluding Section 6.
Up to our knowledge, the continuum-scale phase-transition regimes in Section 5 for realistic argon-methane mixtures have not been computed before using a physically accurate model for the interfacial motion.This contribution is the main novelty of the paper at hand, showing that the approach we have developed for the much simpler single-component cases in [31][32][33] can be readily extended to complex scenarios lacking any knowledge about the proper choice of transmission conditions on the continuum scale.
In the remainder of the introductory section we embed our work into the existing literature for sharp-interface modelling in phase-transition flow.In fact, two-phase fluid flow models with sharp interfaces for a single component are widely investigated, see e.g.[11,12,14,27,38,39] for various modelling approaches in the mathematical realm.As mentioned above, much less is known for multi-component flow.We refer to the works [19,20] which address the closure problem but do not provide a full continuum-mechanical closure.An alternative ansatz [18] builds on the integration of first-order models into phase-field models [9].MD-simulations are readily and since long time used tools for the dynamics of phase boundaries.However, multiscale modelling that connects to the continuum-scale dynamics for compressible liquid-vapor flow has not been done up to our knowledge.The works in e.g.[13,22] indeed support the excellent match of up-scaled quantities if an appropriate closed continuum-scale model exists.

Modelling of Two-Component, Two-Phase Flow
In this section, we first review a continuum-scale model from [6] that governs the description of two-phase flow for a fluid that consists of several components.We focus on argon-methane mixtures with comparably well-known fluid properties.Albeit the fact that the governing equations take the form of first-order balance laws, it turns out that the analytical understanding of hyperbolicity and the properties of solutions for e.g.Riemann problems is quite limited, and cannot serve as the basis for numerical simulation.Throughout this work, we consider all variables in combination with their physical units.This is important for converting microscale quantities to continuum-scale quantities later on.Moreover, the simulations in Section 5 refer to realistic fluid regimes for the argon-methane mixture.

An Isothermal Two-Component Flow Model
We consider isothermal, inviscid, non-reactive fluids consisting of two components, which we mark by an index  ∈ {0, 1}.Later the index 0 indicates argon and the index 1 methane.
To model two-component flow, we choose the so-called class-II model for multi-component flow which has been derived in [6] as a thermodynamically consistent ansatz.It includes Maxwell-Stefan diffusion to account for the crucial frictional interaction between components.The model can be written in terms of a first-order balance law.Precisely, we have for each component  ∈ {0, 1} the equations in  ± () for  ∈ (0,  end ).The primary variables are the partial mass densities   [kg m −3 ] and the partial velocities   [m s −1 ].The constant reference temperature of the mixture is denoted by  [K].The function   =   ( 0 ,  1 ;  ) [J kg −1 ] represents the chemical potential of component  ∈ {0, 1} with respect to the partial mass densities and temperature.We will discuss a specific choice for   below.The term  01 =  10 > 0 denotes the friction factor between the components 0, 1 and is proportional to the reciprocal of the Maxwell-Stefan diffusion coefficients Ð 01 [m 2 s −1 ].More specifically,  01 computes as Here,  [mol m −3 ] denotes the total molar concentration, which is defined by  =  0 +  1 , where   =   /  [mol m −3 ] are the molar concentrations, and   [kg mol −1 ] the molar masses of component  ∈ {0, 1}.The ideal gas constant is given by R ≈ 8.314 J mol −1 K −1 .Another relevant quantity in the context of two-component flow are the mole fractions   [-], which are defined as   =   /,  ∈ {0, 1}.They describe the ratio of the number of 0/1-component particles to the total amount of particles.For an in-depth discussion of friction in multi-component systems we refer to e.g.[25].
To close the system (1) we have to provide expressions for the chemical potentials  0 ,  1 .This means in the isothermal case that an equation of state (EOS) has to be defined which describes the functional dependencies for the chemical potentials   =   ( 0 ,  1 ;  ).For two-component flow we choose the PC-SAFT EOS [15].Its underlying theoretical framework provides us in particular with precise parameter values describing argon-methane mixtures, see Section A. This is the first reason for choosing this mixture in Section 5 on realistic numerical simulations (see Section 3 for the second reason related to molecular-dynamical modelling).We do not give the highly complex functional form of the PC-SAFT EOS.Rather, in Figure 2 we display a plot of the mixture pressure for the argon-methane mixture at  = 110 K as a function of the total density  0 +  1 and the mole fraction  0 of argon.We observe a vapor-state domain and, separated by a spinodal region, a liquid-state domain.Although no proof is available we claim that for sufficiently small temperatures the set {( 0 ,  1 ,  ⊤ 0 ∶=  0  ⊤ 0 ,  ⊤ 1 ∶=  1  ⊤ 1 ) ⊤ ∶  0 ,  1 > 0,  0 ,  1 ∈ R  } is separated into a liquid-state region P − and a vapor-state region P − which define the phases in our case.Having defined a liquid and a vapor phase, we assume that the two-component fluid in  partitions the domain in two subdomains  ± =  ± () separated by a single, sharp phase boundary (), see Figure 1.Just as for the topology of the state space and the splitting, there is up to now no proof that the left-hand side of system (1) is hyperbolic in the state space and elliptic in the spinodal region.Motivated by the situation for single-component, two-phase flow this is exactly what we assume from now on [36].In all the numerical simulations we have not encountered contrary evidence.
Then, the rotational invariance of the Euler-type system (1) ensures that there is a matrix A = A(; ) ∈ R 4×4 such that the solution  = (, ) of the -rotated (one-dimensional) problem defines a planar solution of (1).We recover the original states in U by the operator with  1 , … ,   being the entries of .Now, posing the Riemann problem for the system (4) means to choose initial data of the form where  ± are computed from given Riemann states  ± ∈ U by  ± ∶= ( ± ) ∥ .
The Riemann states  ± ∈ P ± should be understood as traces states in  ± () respectively.The Riemann problem for the -rotated system (4) will play a crucial role in the construction of the heterogeneous multiscale method (HMM) in Section 4. The most convenient situation would be to have an exact Riemann solver for (4) 1 .Motivated by the one-phase case one can then expect that the solution of the Riemann problem ( 4), ( 6) includes exactly one isolated phase boundary connecting two constant states in P ± .Its speed  ∈ R and the values of the adjacent states  * ± ∶= ( * 0,± ,  * 1,± ,  * 0,± ,  * 1,± ) ⊤ ∈ P ± would then provide all the necessary information about the future behavior of the solution locally at  ∈ ().Let us denote such a (perfect) interface solver in terms of the original variables for (1), i.e., using the mapping (3) and the back projection   , by Actually, an exact Riemann solver R ∶ P − × P + × S −1 → P − × P + × R for arbitrary two-phase Riemann data  ± is not available for multicomponent flow (but see [18,19] for some works in this direction).This problem is of course related to the lacking information about hyperbolicity of (1).
We aim to substitute the missing Riemann solver by an interface solver based on MD simulations which will be discussed in the next section.

Molecular Dynamics for Multicomponent Fluids
One big advantage of molecular dynamics (MD) are their versatility to simulate fluid mixtures consisting of more than one fluid component.Fluid components in this context denote different types of molecules that make up the composition of a fluid.
In this section, we review the molecular-dynamical modelling for binary mixtures.Then we extend a numerical method that has been developed in [32] for one-component fluids to the case of binary mixtures, i.e., fluid mixtures consisting of two components.We also discuss the compatibility of MD modelling for argon-methane mixtures with the PC-SAFT EOS introduced in Section 2.1.As before, we indicate the two components by the numbers 0 and 1.

Molecular Dynamics for Binary Mixtures
We focus on fluid components that can be modeled by spherical, rotational invariant Lennard-Jones particles, such as noble gases or relatively simple molecules like methane.In the twocomponent setting, each particle belongs to either component 0, or to component 1.We denote the particle type for the -th particle by   .Accordingly, the mass of each particle is On the atomistic scale the dynamics of an ensemble of  =  0 + 1 particles with  0 ∈ N particles of component 0 and  1 ∈ N particles of component 1 are governed by the ordinary initial value problem [2] for the particle positions   () = (  (),   (),   ()) ⊤ ∈ R 3 , velocities   () ∈ R 3 , and accelerations   () ∈ R 3 of the -th particle, with  ∈ [0,  end ].Here  end > 0 is the MD end time.The MD simulations are carried out in the atomistic-scale units [K] for energy, [ Å] for distances, and [u] for mass.The crucial ingredients for the MD model ( 8) are the Lennard-Jones potentials   with The Lennard-Jones parameters  01 =  10 ,  01 =  10 depend on the corresponding parameters  0 ,  0 ,  1 ,  1 for single-component potentials.If both interacting particles belong to the same component  ∈ {0, 1}, the parameters are given by   ∶=   ,   ∶=   .If they belong to different components, the Lorentz-Berthelot combination rules [5,26] suggest for some scaling numbers  [-] and  [-] In Figure 3, the Lennard-Jones potentials for different parameters are plotted.
As discretization for (8) we use the classical Velocity-Verlet algorithm [16] which amounts to performing the following steps with a given time step Δ > 0,  end = ⌈ end /Δ⌉ times: Remark 1 (Reduced units, time scales, and cutoff-potential).
(i) The MD simulations using the Velocity-Verlet algorithm (11) are done in reduced units according to Table 1 below.With this procedure we follow [2].It avoids floating point over-/underflow in performing the numerical algorithm.
The number of particles  , the time step Δ (in reduced units) and the number of time steps  end will enter among others as parameters in Algorithm 2 for solving an MD Riemann problem.The precise values of the method's parameters are listed in Table 4.These parameters are chosen in a way that allows a reliable extraction of Riemann data in the final multiscale method which acts on the continuum scale.Note that we expect the MD Riemann problem to have a self-similar solution such that we can extract data at any point in time.Without any information about numerical errors the ad-hoc choice of these parameters remains a weakness of the entire multiscale approach.
(ii) The application of the algorithm (11) leads to a quadratic complexity in terms of the total particles number  .This is due to the particle-wise interaction in (9).A remedy are cutoffpotentials.For a two-component mixture we substitute   for ,  ∈ {0, 1} by Like all numerical parameters for the MD simulations the actually used value for the cutoff radius  cutoff can be found in Table 4 in the appendix.The error introduced by the cutoff functional will be corrected by the introduction of long-range interaction forces on the discrete level.We refer to [28] for a special choice for two-component flows.The two solid-line graphs display intracomponent interaction, and the dashed one the interaction potential between the two components, using the combination rules (10).

Reduced units SI units
We are interested in running MD-simulations that provide an approximation of the solution of the Riemann problem (4), (6).For this purpose the -rotated Riemann initial data  ± from (6) have to be represented as initial data for (8), i.e., we need an initial particle and velocity distribution in 3D.In [32] we have suggested a linked-cell algorithm to realize this transfer of single-component, two-phase constant Riemann data into a cubic domain Σ that consists of two cuboids Σ ± separated by the atomistic interface.The algorithm includes a thermalization process that makes sure that the MD-particle and particle velocity distributions matches the constant reference temperature  , see ( 14) below for the definition of atomistic-scale temperatures.The linked-cell algorithm can be readily applied to the two-component case and is not further detailed.
In the next step, the system (8) is solved numerically by the scheme (11).We can again utilize a procedure from [32] to determine the approximate position of the atomistic interface and in particular an approximate interface speed  ∈ R like in the exact interface solver R from (7).To obtain all output data in (7) we need to find approximations of adjacent continuum-scale states  * ± ,  * ± , respectively.Thus, we need to define continuum-scale quantities from the MD-simulation data.Let us consider a system of  particles.The component-wise mass concentrations ρ0 (, ), ρ1 (, ) [u Å−3 ], and momentum concentrations m0 (, ), m1 (, ) [K 0.5 u 0.5 Å−3 ] can be computed by the Irving-Kirkwood formulas in the following way, for  ∈ {0, 1}, where   ,   ,   are the mass, position and velocity of the -th particle.The local temperature distribution of the mixture can be formally understood as where  = (, ) [K 0.5 u −0.5 ] denotes the local barycentric average velocity of the mixture.The local average velocity can be realized by e.g.
where   ∶ R → R is a smooth kernel with length-scale parameter  > 0. It cannot be expected that the temperature field T from ( 14) equals the constant temperature  in the isothermal system (1).This can be partially circumvented by applying a thermostat, but we will observe in Section 5 that substantial discrepancies can occur for the vapor-phase region.
Using the formulas ( 13 where vol(Σ * ± ) denotes the volume of the sampling region Σ * ± , and |  ∈ Σ * ± ,   = | the number of component- particles inside it.We can obtain the local interface states for every time step during an MD simulation.For the final outcome we perform time-averaging of these instantaneous states over a fraction  −smpl of the total simulation time  end .Note that the conversion of MD quantities and continuum-scale quantities has to account for the unit conversion given in Table 1.
We summarize all steps for solving an MD-Riemann problem in the following algorithm, see also Algorithm 4 in the context of the multiscale method.

Notes on the Argon-Methane Mixture
Up to now, we considered a generic binary mixture.In the following, we focus on two-component mixtures consisting of argon and methane.In this case, parameters for the Lennard-Jones potential (9) can be found in [41], and even experimental data are available [35].For the Lennard-Jones parameters in (10) we take the values from [41] listed in Table 2.The corresponding combination parameters in (10) are  = 1.001 41 and  = 0.964 00, see also [41].For methodological background on the used approach we refer to [42].
Remark 3 (Temperature parameterization of Lennard-Jones parameters).The chosen Lennard-Jones parameters are given in [41] for temperatures  = 111 K (Argon) and  = 141 K (Methane).The combination parameters refer to the temperature  = 115 K and have also been computed in [41].In this work the quantities have been shown to be in good fit with available experimental data in [35].Up to our knowledge this setting is commonly used for temperatures well below the critical temperature of the mixture.Notably, our HMM method applies only for settings where the phase boundary can be expected to be a sharp interface separating the liquid and the vapor regions.Close to the critical temperature this separation property would fail.
The studies in the literature refer solely to the homogeneous case, i.e., they do not span the entire two-phase state space.In [28] it has been shown by extended MD simulations for homogeneous data covering the whole state space, that the predictions of the PC-SAFT theory for the chemical potentials  0 ,  1 are in good quantitative agreement with the MD simulations using the Lennard-Jones potentials for argon-methane mixtures.This justifies to use the PC-SAFT chemical potentials for the continuum-scale computations in the HMM, that will be proposed in Section 4.  [41] and species masses in reduced units for the argon-methane mixture, cf.Table 1.

A Heterogeneous Multiscale Method for Isothermal Two-Component, Two-Phase Flow
Following [32] for single-component flow with phase transitions we introduce a corresponding heterogeneous multiscale method (HMM) that combines the isothermal two-component flow model on the continuum scale (2) with two-component MD simulations described in Section 3 for the microscale description of the interface motion.
Here, we focus again on a fluid mixture consisting of argon and methane, see Section 2. The component-specific quantities are marked by the indices 0 for argon, and 1 for methane.The two fluid phases are again denoted by the subscripts − and + for the liquid and vapor phase.First, we formulate on the basis of Algorithm 2 the MD-based microscale interface solver  MD which provides an approximate solution for Riemann initial data (Section 4.1).Using an MDbased interface solver directly in place is still too computational expensive.Therefore we suggest the machine-learned surrogate solver R  in Section 4.2, which exploits the fact that the interface solver  MD can be understood as a nonlinear function mapping finite-dimensional Riemann data to a finite-dimensional space including the phase boundary speed and the state values adjacent to the phase boundary.Finally, in Section 4.3 we present the complete HMM for two-component, two-phase flow.For the sake of simplicity we ignore in the entire section boundary conditions on .

The Atomistic-Scale Interface Solver 𝑅 MD
To reduce the input space of the atomistic-scale interface solver from Algorithm 2 we exploit that the continuum-scale system (1) is invariant with respect to velocity shifts.We define the barycentric velocities and choose v ∶=  − as the reference velocity.Furthermore, we introduce the relative velocity for each phase Assuming we know the densities  0,± ,  1,± , we can compute ( 0,± ,  1,± ) from ( ± ,  rel,± ) and vice versa.In the sequel this conversion is done implicitly to simplify the notation.
The following algorithm implements the final MD-interface solver  MD ∶ P − × P + × S −1 → P − × P + × R substituting the exact Riemann solver R.
Result:  MD ( − ,  + ; ) = ( * − ,  * + , ) Due to the velocity shift described before we can always achieve that the component  0,− can be assumed to vanish when Algorithm 4.1 is executed within Algorithm 4.
For the full multiscale model, the surrogate model R  replaces the molecular-scale interface solver  MD in step 4 of Algorithm 4. We note that the implementation of R  relies on the reduceddimensional state spaces in the rotated Riemann problem (4) and using the invariance with respect to the reference velocity  − in (16).The final surrogate solver offers large computational gains: It takes only around 0.10 ms for a single evaluation.The MD simulations (Algorithm 4) on the other hand needs 14 min to 17 min.

Data Set Generation
To prepare the surrogate interface solver, we have to generate a data set  for the input-output relation of the microscale interface solver  MD .The range of the ±-phase densities  0,± ,  1,± for each component is defined by the convex set, that is formed by the points given in Table 3.
The barycentric velocity  + ranges from  min = −750 m s −1 to  max = 750 m s −1 .The relative velocities  rel,± are bounded by  rel,min = −500 m s −1 and  rel,max = 500 m s −1 .Taken all together, the resulting convex set forms the input bounding domain  in .
Using distance-maximizing sampling similar to [34], we generate  data = 12 000 samples in  in , while exploiting the reduced dimensionality, due to rotational invariance and  − ≡ 0. Consequently, we generate the input data set  Note that we run three MD simulations for each input data point.The data set  is visualized in Figure 4.It is digitally archived at [29].With the run time of 14 to 17 min for a single, twocomponent MD simulation, the generation of the whole data set  took approximately 3200 h of computing time, which can be easily split among several of machines to decrease the real time until all data points are sampled.

Neural Network Training
To train the neural-network surrogate solver R  , we use a standard training procedure, albeit with some model-specific parameters as described in the following.The network is comprised of 5 hidden layers with 60 nodes each.The data set  from Section 4.2.1 is split into a training data set  train with 10 800 samples and a validation data set  val with 1200 samples.Finally, we note that we use a CRes-neural network, which has been developed in [30] to incorporate conservation of selected state variables into the surrogate interface Riemann solver.For two-component, two-phase flow we decided to preserve the mass conservation across the interface.

The Complete HMM
Using the machine-learned interface solver R  we present the final heterogeneous multiscale method and the interplay of each of its components.To see the integration of the interface solvers in detail we describe the HMM in the framework of the numerical discretization, a moving-mesh finite-volume method.The method has been introduced for two space-dimensions in [7] and open-source code for two and three space-dimensions can be found in [1].We partition the continuum-scale time interval [0,  end ] according to 0 =  0 < … <   =  end , and start with an initial mesh T(0) = { 0  ∶  ∈  0 } for  that is assumed to consist of simplices   indexed via some index set  0 .Let the mesh parameter ℎ > 0 be given by the maximum of all edges' length.The mesh T(0) is supposed to include a connected set of mesh facets of the simplices, which form the approximation of the initial interface  ℎ (0), and thus discrete liquid and vapor bulk domains  ℎ,± (0).The initial datum  0 from ( 2) is projected to the mesh T(0) giving We assume that we have  0  ∈ P ± for  0  ∈  ℎ,± (0).To introduce the finite-volume time-stepping let us assume that for some discrete time   ,  ∈ {1, …  − 1}, we have the same situation, i.e., a simplicial mesh T(  ) = {   ∶  ∈   } for , a discrete partition  =  ℎ,− (  ) ∪  ℎ (  ) ∪  ℎ,+ (  ) with discrete interface  ℎ (  ) = {   ∶ ,  ∈   ,    ∈  ℎ,− (  ),    ∈  ℎ,− (  )} of connected facets, and a set of cell averages {   } ∈  with the property    ∈ P ± for    ∈  ℎ,± (  ).Furthermore, for any facet from  ℎ (  ) with normal    ∈ S −1 (pointing to the vapor domain  ℎ,+ (  )) we can compute ( *  ,  *  ,    ) = R  (   ,    ;    ).Using the computed speeds    , the moving-mesh finite-volume method from [7] delivers a mesh T( +1 ) = { T +1  ∶  ∈ Ĩ +1 }, a discrete partition  = Ωℎ,− ( +1 ) ∪ Γℎ ( +1 ) ∪ Ωℎ,+ ( +1 ), and the family { Ũ +1  } ∈ +1 with the properties as in time-step   .The new mesh T( +1 ) evolves from T( +1 ) by affine shifts.To avoid small cells that would deteriorate the time-step this mesh deformation is followed by a re-meshing leading finally to the new approximations at time  +1 , i.e., the mesh and the finite-volume approximations A particular point in the moving-mesh finite-volume method is the choice of the numerical flux functions.They can be chosen arbitrarily for facets not in  ℎ (  ) but have to be Godunov fluxes across  ℎ (  ) using the adjacent states  *  ,  *  .This avoids the occurrence of states not belonging to P ± .For more details and the control of the explicit time-stepping we refer to [7].We conclude the section with the complete HMM for two-component, two-phase flow given in algorithmic form.
Remark 6 (Discretization of source and non-conservative terms).The moving-mesh finite-volume method from [7,32] uses explicit Euler time-stepping and deals with first-order systems in divergence form only.The discretization of the 0-th-order Maxwell-Stefan diffusion terms is handled by simple evaluation of the cell averages.More complicated is the discretization of the gradients ∇  in the right-hand-side term of the continuum-scale system (1).
In one space-dimension we simply apply central finite-differences in the bulk phases.At the cells adjacent to the interface we use left-/right-sided finite differences, to avoid computing the gradient over the discrete phase boundary.
In two and three space-dimensions, for the -th time-step, we approximate the gradient ∇  in a simplex    by linear reconstruction.The stencil for the reconstruction includes those neighbors    of    that share a surface with    and belong to the same phase domain as    .In that way, we avoid mixing the phases during the reconstruction.The reconstruction is done by solving the linear least squares system with where   ∈ R  denotes the cell center of    ,  , ∈ R the cell value of the chemical potential, and  ∈ N the number of neighbor cells in the stencil.The solution  ′  of ( 20) is used as an approximation of the gradient ∇  in cell    .

Numerical Simulations for the HMM
In this section we present a series of numerical results to validate the multiscale method that has been introduced in Section 4 to solve the free boundary value problem for the two-component, two-phase flow model (1).
The model and numerical parameters used for the simulations in this section, if not otherwise stated, are found in Appendix A. We choose the Maxwell-Stefan diffusion coefficient to be Ð 01 = 1.0.We refer to [23] for the computation of specific Maxwell-Stefan diffusion coefficients leading to much smaller values.Numerical tests for simplified settings have shown that the influence of Ð 01 is minor for the overall dynamics, which justifies our choice.

One-dimensional Simulation Results
We consider one-dimensional Riemann problems and compare the multiscale simulation results with MD simulations on the entire domain as reference.The section serves as validation for the multi-dimensional droplet simulations in the subsequent Section 5.2 and Section 5.3.
In the multiscale solution it can be observed that the liquid phase expands to the right side, while the total density inside the liquid phase decreases.Near the interface, argon moves from the liquid into the vapor phase, whereas methane accumulates inside the liquid, and decreases in the vapor phase.By comparing the multiscale and the MD solution, we see that the interface position is captured accurately.Furthermore, on the liquid side, the MD simulation and the multiscale model behave qualitatively in the same manner.On the vapor phase side, we see that the argon-component of the multiscale simulation is close to the MD simulation.In contrast to that, substantial deviations in the methane-component can be observed.We assume that these deviations can be attributed to the isothermal continuum model which does not capture the whole range of the MD simulation.The latter is not perfectly isothermal allowing for small temperature variations.We observed the same effect for isothermal flow of one component in [32].Finally, the deviations in the low-density vapor phase might be caused by poor sampling on the MD scale, as few particles of each component are present near the interface.This problem can be solved by increasing the number of particles, as well as the size of the computational domain (see Section 3) for the MD simulations.Example 8 (Colliding vapor wave).In the second example, we simulate a vapor wave that collides with a liquid argon-methane mixture.The corresponding Riemann initial data are  0,− = 440 kg m −3 ,  1,− = 280 kg m −3 ,  0,− = 0 m s −1 ,  1,− = 0 m s −1 for the liquid phase, and  0,+ = 20 kg m −3 ,  1,+ = 2 kg m −3 ,  0,+ = −50 m s −1 ,  1,+ = −50 m s −1 for the vapor phase.The multiscale simulation results are plotted in Figure 6, alongside with their respective MD simulation.It can be seen that the vapor wave transmits into the liquid phase, and increases the liquid density slightly.Furthermore, the wave speeds, as well as the interface speed, are captured very well by the multiscale model.As in Example 7 we observe deviations for the velocity in the vapor phase domain.
The simulation results are depicted in Figure 7 and Figure 8.In the first figure, the componentwise densities and velocities are shown, and in the second figure the argon-mole fraction is plotted.In the beginning, the liquid and vapor phase are not in equilibrium, resulting in small oscillations of the droplet.The oscillations can be clearly observed in the time evolution of the −-averaged quantities, as illustrated in Figure 9. Then the vapor wave hits the droplet, causing a ripple moving through its surface and finally pushing it through the vapor atmosphere.Throughout the simulation, argon accumulates inside the liquid phase, leading to a growth of the droplet -see Figure 9.
For the multiscale simulation results we refer to Figure 10 and Figure 11.Both figures show the three-dimensional solution visualized on the plane through  = 0.In Figure 10 the componentwise densities and velocities are depicted.We display also the (projected) mesh to demonstrate that the interface is resolved on the discrete level.In Figure 11 the phase boundary is shown alongside the barycentric velocity magnitude.We observe that the droplet is hit by a shock wave, which results in a ripple that runs over the liquid surface.Additionally, the momentum from the vapor causes the droplet to move to the right side.Considering the two fluid components, we see that the methane concentration inside the liquid droplet decreases slowly, whereas the amount of liquid argon slowly increases.Next, we provide an overview over the simulation run times.A single time step of the simulation takes on average 2.7 s.This splits approximately into 1.0 s for the moving mesh operations (see Section 4.3), and 1.7 s for the finite volume part of the simulation.The mesh consists of 2.7 • 10 5 cells, with 2.0 • 10 3 interface facets averaged over all 100 000 time steps.In total the whole simulation takes 75 h on a single desktop computer.Note that in this simulation the need for a surrogate solver becomes apparent.If we would not have employed a surrogate, we would need to perform an MD simulation (taking 14 min to 17 min) for every interface facet at every single time step.For a single time step, that would accumulate to approximately 500 h of computational time (albeit parallelizable).Compared to that, the surrogate solver needs 0.102 ms for a single evaluation.The offline phase of the surrogate solver takes approximately 3200 h plus training.The majority of this workload (generating the training data set) is however embarrassingly parallel.Consequently, even if we include the offline phase, using a surrogate solver pays off after a few time steps.

Conclusions
The main result of our study is a proof-of-concept validation for a HMM that has been designed to model the dynamics of compressible two-component flow with liquid-vapor phase transition.The behavior of the numerical solutions on the continuum scale matches the MD simulations, with even quantitatively correct results in the liquid phase domains.Improving the results in the vapor phase seems to require more computational effort in the sampling procedure and the thermalization process.It is noteworthy that we do not consider a generic fluid mixture, but use a real EOS for argon-methane mixtures.Another significant outcome is that the HMM is able to simulate two-and three-dimensional droplets in case of complex interactions.At the same time, both small and large-scale deformations of the (sharp) phase boundary can be rendered by the moving-mesh finite-volume method during the numerical simulations.Up to our knowledge, this is the first time that compressible mixtures of real fluids with resolved sharp interfaces are simulated using an MD-based and neural-network accelerated multiscale model.In particular, we are able to simulate phase transitions of fluid mixtures -without the need to prescribe some (ad-hoc) closure relations at the interface.We think that the HMM in combination with an ML surrogate modelling paves the way to handle more complex flow scenarios.Surface tension has not been taken into account, despite its importance for two-phase evolution.Here, it is possible to account for curvature effects on the continuum-scale by adding a geometric term to the momentum equation [36].However, the surface tension coefficient might then require additional MD simulations like for the EOS.On the expense of more MD simulations we expect that the extension to mixtures with more than two components is straightforward.A conceptual advantage of the underlying moving-mesh finitevolume method is that tailored numerical schemes can be applied in the liquid and the vapor phase.This would allow for benefitting from e.g.implicit time-stepping in the low-Mach liquid domain and -more important for mixtures -to account for large Maxwell-Stefan coefficients.

8 Figure 2 :
Figure 2: Pressure of the PC-SAFT EOS for an argon-methane mixture at  = 110 K.The left-hand (right-hand) region displays vapor (liquid) states.

Figure 5 : 3 density [kg m − 3 ]Figure 6 :
Figure 5: Multiscale simulation for the isothermal, two-component, two-phase flow model in one space-dimension, overlaid over the corresponding MD simulations.

Figure 7 :Figure 8 :
Figure 7: Two-dimensional multiscale simulations.Each sub-figure depicts the densities  0 ,  1 and velocities  0 ,  1 of each component at various time steps.The upper part of each sub-figure shows  1 ,  1 for methane, and the lower part  0 ,  0 for argon.

Figure 9 :
Figure 9: Time evolution of liquid phase-averaged quantities for the twocomponent two-phase flow multiscale simulation in two spacedimensions.

Figure 10 :Figure 11 :
Figure 10: Three-dimensional multiscale simulation of the two-component twophase flow model.The upper part of each sub-figure shows  1 ,  1 for methane, and the lower part  0 ,  0 for argon, at various time steps.

Table 3 :
Density corner points for the two-component (argon-methane) model input data set.