This version is available at https://doi.org/10.14279/depositonce-9565
This work is licensed under a CC BY-NC-ND 4.0 License (Creative
Commons Attribution-NonCommercial-NoDerivatives 4.0 International). For
more information see https://creativecommons.org/licenses/by-nc-nd/4.0/.
Terms of Use
Hitz, T., Heinen, M., Vrabec, J., & Munz, C.-D. (2020). Comparison of macro- and microscopic solutions of
the Riemann problem I. Supercritical shock tube and expansion into vacuum. Journal of Computational
Physics, 402, 109077. https://doi.org/10.1016/j.jcp.2019.109077
Timon Hitz, Matthias Heinen, Jadran Vrabec, Claus-Dieter Munz
Comparison of macro- and microscopic
solutions of the Riemann problem I.
Supercritical shock tube and expansion
into vacuum
Accepted manuscript (Postprint)Journal article |
Comparison of Macro- and Microscopic Solutions of the
Riemann Problem I. Supercritical Shock Tube and
Expansion into Vacuum
Timon Hitza, Matthias Heinenb, Jadran Vrabecb, Claus-Dieter Munza
aInstitute of Aerodynamics and Gasdynamics, University of Stuttgart, 70550 Stuttgart,
Germany
bThermodynamics and Process Engineering, Technical University of Berlin, 10623
Berlin, Germany
Abstract
The Riemann problem is a fundamental concept in the development of nu-
merical methods for the macroscopic flow equations. It allows the resolution
of discontinuities in the solution, such as shock waves, and provides a power-
ful tool for the construction of numerical flux functions. A natural extension
of the Riemann problem involves two phases, a liquid and a vapour phase
which undergo phase change at the material boundary. For this problem, we
aim at a comparison with the macroscopic solution from molecular dynamics
simulations. In this work, as a first step, the macroscopic solution of two
important Riemann problem scenarios, the supercritical shock tube and the
expansion into vacuum, were compared to microscopic solutions produced by
molecular dynamics simulations. High fidelity equations of state were used
to accurately approximate the material behaviour of the model fluid. The
results of both scenarios compare almost perfect with each other. During
Preprint submitted to Journal of Computational Physics October 24, 2019
Accepted manuscript of: Hitz, T., Heinen, M., Vrabec, J., & Munz, C.-D. (2020). Comparison of macro-
and microscopic solutions of the Riemann problem I. Supercritical shock tube and expansion into
vacuum. Journal of Computational Physics, 402, 109077. https://doi.org/10.1016/j.jcp.2019.109077
© 2019 This manuscript version is made available under the CC-BY-NC-ND 4.0 license
https://creativecommons.org/licenses/by-nc-nd/4.0/
the vacuum expansion, the fluid obtained a state of non-equilibrium, where
the microscopic and macroscopic solutions start to diverge. A limiting case
was shown, where liquid droplets appeared in the expansion fan, which was
approximated by the macroscopic solution, assuming an undercooled vapour.
Keywords:
Riemann Problem, Real Gas, Supercritical fluid, Vacuum Riemann
Problem, Finite Volume, Molecular Dynamics Simulation
1. Introduction
The Riemann problem poses an initial value problem for a conservation
equation with piecewise constant initial data. In the numerical approxima-
tion of the compressible flow equations by finite volume (FV) schemes, it is
one of the basic building blocks. Because the FV schemes approximate the
integral conservation equations, they allow for discontinuities between the
grid cells. In his pioneering work, Godunov [1] proposed to use the solution
of the Riemann problem to determine the flux between the grid cells. Taking
into account the information about the waves generated by a discontinuity
establishes the robustness of such FV schemes, e.g., at shock waves or in
under-resolved regions of compressible fluid flows. Later, several approxima-
tions of the solution of the Riemann problem have been developed, which
preserve the advantageous properties of the Godunov scheme while reduc-
ing the computational effort. A comprehensive presentation of the Riemann
problem and its approximations is given by Toro in his book [2].
The solution of the Riemann problem in gas dynamics is well-known under
certain convexity constraints for the underlying equation of state (EOS).
2
Here, the solution consists of four constant states that are separated by three
waves. The outer waves are shock waves or rarefaction waves, associated with
the characteristic wave speeds v−cand v+c. The characteristic field of flow
transport vis linearly degenerate and exhibits a contact discontinuity, with
a constant pressure and velocity across it. Menikoff and Plohr [3] describe
wave structures that may occur if the standard convexity conditions are not
satisfied. Beside the classical waves, shock wave, rarefaction wave and contact
discontinuity, anomalous entities such as splitted waves or composite waves
may form. Near a phase transition, an additional wave that corresponds to the
phase interface appears in the wave structure of the Riemann problem. This is
possible because the hyperbolicity of the Euler equations breaks down in the
spinodal region of the state space and imaginary eigenvalues occur. To avoid
such non-physical states in the macroscopic solution, information from the
local thermodynamic behavior is needed and has to be incorporated into the
macroscopic solution of the Riemann problem. Some additional information
about local thermodynamic processes has to be formulated. One approach is
based on a so-called kinetic relation that controls the mass transfer across
the phase interface (see, e.g. LeFloch [4]). But the situation is not fully clear
up to now and even more unclear is the numerical modelling.
In the last decades, molecular dynamics simulations became more pow-
erful due to advanced algorithms and increasing computer efficiency. They
resolve matter down to the atomistic scale by treating every molecule individ-
ually as a mechanical object. A central role is being played by the force field
that describes the intermolecular interactions, like dispersion or repulsion
due to Pauli exclusion. Once assigned to the molecules, these interactions
3
not only determine their individual trajectories under given boundary condi-
tions, but also contain, for a molecular ensemble, the entire set of thermody-
namic properties in a consistent manner. MD has a sound physical basis and
is thus a versatile approach that can be applied to a wide range of problems,
like adsorption, diffusion or evaporation. However, a drawback is the asso-
ciated computational effort. Newton’s equations of motion must be solved
numerically with a time step on the order of femtoseconds and the spatial
resolution is on the ˚
Angstr¨om scale so that studying macroscopic processes
is challenging.
The objective of our research activities is to compare the macroscopic
modelling and the microscopic thermodynamic behavior. Hence, we compare
molecular dynamics (MD) simulations for Riemann problem scenarios with
macroscopic solutions. The most interesting situations are those, of course,
in which phase transitions occur and the macroscopic solutions depend on
the local thermodynamic modelling. To allow for a full correspondence of
the macro- and microscale solutions, the present MD simulations were based
on the Lennard-Jones model fluid with a truncated and shifted potential
(LJTS). For this model fluid, two highly accurate equations of states for the
macroscopic simulations are available: the empirical EOS by Thol et al. [5]
and the semi-empirical EOS by Heier et al. [6], which have similar mathemat-
ical properties as the van der Waals EOS. This gives the chance to directly
compare two topics. First, the macroscopic and the microscopic Riemann
problem solutions directly and then the molecular dynamics Riemann prob-
lem solution with results by a finite volume method on some appropriate
grid. We start here with a validation of this approach and consider Riemann
4
problems in the supercritical regime without phase transition. The expansion
of a fluid into vacuum leads to simple simulation conditions for MD. Hence,
this is our second class of benchmark problems. Here, we see the following: if
the initial conditions are such that due to the strong pressure drop in the rar-
efaction fan, liquid droplets appear in MD simulations, then the macroscopic
simulations become more difficult to obtain due to instabilities. We assume
that we have to change the modelling of the macroscopic Riemann problem
solution and have to take phase transition in our numerical flux calculation
into account.
We restrict ourselves in this first paper to the description of the methods
and the validation of this approach. The interesting topic of the approxi-
mation of phase transition in macroscopic simulations is only touched. Our
research activities are embedded in the DFG collaborative research center
SFB-TRR 75 ”Droplets under extreme ambient conditions”. Hence, we have
the possibility to continue this research in the coming years, especially look-
ing at situations, in which phase transitions do occur. MD results for these
situations will be documented and compared with macroscopic simulations
that have to take phase transitions appropriately into account.
The structure of the paper is as follows. In Section 2, the continuum
equations, the EOS and the Riemann problem are introduced. The numerical
methods are discussed in Section 3. Results for a supercritical shock tube
Riemann problem and several vacuum Riemann problem scenarios are given
in Section 4, followed by a conclusion in Section 5.
5
2. Theory
2.1. Equations of the Continuum System
The compressible Euler equations in one spatial dimension are written in
conservation form as
Ut+F(U)x= 0 ,(1)
with
U=
ρ
ρv
ρe
,F=
ρv
ρv2+p
v(ρe +p)
.(2)
The variables are density ρ, velocity vand specific total energy e=+1
2v2,
which is composed of the internal part and the kinetic part 1
2v2. In addition
to the conservative state vector U, the primitive variables can be written as
a state vector W= (ρ, v, p)>, where pdenotes the thermodynamic pressure.
An EOS is required as a closure relation between the variables pressure,
density and internal energy
=(ρ, p) ; p=p(ρ, ).(3)
In this study, we consider the Lennard-Jones model fluid with a truncated
and shifted potential (LJTS). Two highly accurate EOS are available in this
case: The empirical EOS by Thol et al. [5] and the semi-empirical EOS by
Heier et al. [6]. Both equations are fundamental EOS in terms of the reduced
Helmholtz free energy α. It can be decomposed into the contribution of a
monatomic ideal gas, α0and a residual part, αr. Note that the superscript 0
always refers to a state described by the ideal gas law of a monatomic gas.
The Helmholtz free energy is given as a function of reduced density δ=ρ/ρc
6
and inverse reduced temperature τ=Tc/T, where the subscript cindicates
the critical point,
a(ρ, T)
RT =α(δ, τ) = α0(δ, τ) + αr(δ, τ).(4)
Using derivatives of Eq. (4), any thermodynamic quantity of interest can be
obtained directly. An overview of the necessary relations is provided in [7, 8].
Since these fundamental EOS are given in terms of density and temperature,
the closure condition for the Euler equations, Eq. (3), is not available in
a closed form. For that purpose, a Newton-Raphson method was used to
calculate the temperature iteratively for a given density and internal energy,
provided by the CFD method. With known density and temperature, all
other thermodynamic quantities can be calculated.
The EOS used in this study apply two different approaches expressing
the residual part αr. Thol et al. [5] provided a multiparameter EOS using the
framework of Span [9]. Here, the residual Helmholtz energy was described
by a sum of polynomial, exponential and Gaussian terms. The coefficients
were determined by a fitting procedure on the basis of reference data from
molecular dynamics simulations. This model is referred to as LJTS EOS in
this study. It is in excellent agreement with reference data, but the evaluation
of Eq. (4) is computationally costly. Furthermore, the occurrence of multiple
Van-der-Waals loops in the two-phase region may cause difficulties for root
finding algorithms. The PeTS EOS by Heier et al. [6] considers a framework
based on perturbation theory [10]. The residual Helmholtz energy was split
into the contribution of the reference potential and the contribution due to
perturbation. Compared to the LJTS EOS, its evaluation is computationally
cheaper and the presence of a single Van-der-Waals loop simplifies the usage
7
of root finding algorithms. However, the PeTS EOS is somewhat less accu-
rate. Both EOS provide formulations [7, 8] to calculate the thermodynamic
properties pressure p, internal energy and speed of sound c.
To handle the transition to vacuum, the asymptotic limits of the Helmholtz
energy for vanishing reduced density and reduced temperature, δ→0 and
τ→+∞, have to be considered. In this case, the residual Helmholtz energy
of both EOS approaches zero
lim
δ→0
τ→+∞
αr= 0,(5)
and only the perfect gas contribution remains. Both EOS use the same for-
mulation for this perfect gas contribution α0corresponding to a monatomic
gas with an isobaric heat capacity c0
p/kB= 2.5, where kBis the Boltzman
constant. The perfect gas contribution reads as
α0(δ, τ) = lnδ+ 1.5 lnτ+c1τ+c2.(6)
The coefficients c1, c2were specified by Thol et al. [5] and Heier et al. [6] such
that the perfect gas contributions to enthalpy and entropy vanish h0
0=s0
0= 0
at the reference state (T0, p0) = (0.8,0.001),
c1=−2.5
τ0
and c2= 1.5−ln δ0−1.5 ln τ0.(7)
The asymptotic limits for pressure and internal energy are
lim
δ→0
τ→+∞
p= lim
δ→0
τ→+∞
ρRT (1 + δαr
δ) = 0,(8)
lim
δ→0
τ→+∞
= lim
δ→0
τ→+∞
RT τα0
τ+αr
τ=RTcc1.(9)
8
The speed of sound cis defined by
c2=∂p
∂ρs
=∂p
∂ρT−∂p
∂T ρ∂s
∂ρT∂s
∂T −1
ρ
.(10)
This formulation also falls back to the classical speed of sound of a perfect
gas, which becomes zero for a vanishing temperature
lim
δ→0
τ→+∞
c2= lim
δ→0
τ→+∞
5
3RT = 0.(11)
2.2. Riemann Problem
The Euler equations (1) are a hyperbolic system of conservation laws with
the real eigenvalues of the Jacobian of the flux
λ1=v−c, λ2=v, λ3=v+c. (12)
The initial value problem for the one-dimensional Euler equations, consisting
of constant initial states, separated by a discontinuity at x= 0, poses the
Riemann problem
U(x, t = 0) =
ULfor x < 0
URfor x > 0.
(13)
Its solution is self-similar with respect to x/t and consists of four constant
states UL,U∗
L,U∗
R,UR, separated by elementary waves. The outer waves are
either shock or rarefaction waves and the intermediate wave is a contact
discontinuity. Each of the three wave types is attached to an eigenvalue of
the system of equations. A typical situation is shown in the diagram in Fig.
1. The shock and rarefaction waves are non-linear waves and associated with
λ1,3, while the contact discontinuity is associated with λ2. The latter is a
9
linearly degenerate discontinuity in mechanical equilibrium, i.e. the velocity
and pressure are equal on both sides. Consequently, there is no macroscopic
mass flux across the contact discontinuity so that it can be considered as a
material interface. The shock wave is a non-linear discontinuity where the
Rankine-Hugoniot conditions apply. The rarefaction wave is a continuous
solution where the so-called generalized Riemann invariants remain constant.
As a consequence, the solution in a rarefaction wave is isentropic.
For the perfect gas, a solution has been proposed by Godunov [1] using
a fixed point iteration. An extension to real gas EOS that remain smooth
and convex has been introduced by Colella and Glaz [11] and further devel-
oped by Kamm [12]. The general approach is to find iteratively the velocity
equilibrium across the contact discontinuity as a function of the unique pres-
sure of the inner states via the Newton or the Secant method. Using this
pressure, the inner states U∗
Land U∗
Rare calculated by applying the rela-
tions across the outer waves, i.e., the Rankine Hugoniot conditions or the
isentropic relations.
A special type of the Riemann problem is the shock tube problem, which
is of elementary importance in gas dynamics. It describes two states in a
tube that are separated by a membrane. A dense state under higher pressure
is adjacent to a dilute state under low pressure. Both states are initially at
rest. At time t= 0, the membrane is removed and different waves propagate
through the domain: A shock wave travels into the dilute fluid and a rarefac-
tion wave travels into the dense fluid with the wave propagation speeds SL
and SR, respectively. Between these waves, two constant states appear that
10
are separated by the contact discontinuity that travels with the speed S∗.
The corresponding solution is shown in Fig. 1.
x
t
UR
UL
U∗
R
U∗
L
SLSR
S∗
Figure 1: Solution structure for the shock tube problem where the high pressure is ini-
tially on the left and the low pressure on the right. Three waves divide four
constant states: UL,U∗
L,U∗
R,UR.
Another special type of Riemann problem is the vacuum Riemann prob-
lem. In contrast to the shock tube problem, one of the states is a vacuum, i.e.
U= (0,0,0)>. Because the Euler equations are not valid in the vacuum, this
scenario has to be considered as a free boundary problem, where the moving
boundary represents the expansion front of the continuum phase into vac-
uum. A self-similar exact solution can be found for the perfect gas as shown
by Halter and Martensen [13] and further discussed in Refs. [14, 15, 2].
If the vacuum is considered on the right side, the solution contains a
single (left) rarefaction wave, cf. Fig. 2. The right wave vanishes and the
contact discontinuity attaches the right boundary of the rarefaction wave.
Integration of the generalized Riemann invariant across the rarefaction wave
to the vacuum state yields a front speed of the expansion wave that is fully
11
determined by the initial continuum state on the left
S∗
L=vL+2cL
γ0−1,(14)
where γ0=c0
p/c0
Vis the isentropic coefficient of a monatomic ideal gas and c0
p,
c0
Vare the heat capacities at constant pressure and volume of a monatomic
ideal gas, respectively. Note that the wave speed S∗
Ldescribes the propagation
speed of the outer most right wave of the left rarefaction wave, that lies
adjacent to the star state, cf. Fig. 2. The solution reads as
W(x, t) =
WLif x
t≤vL−cL
Wfan = (ρfan, vfan, pfan)>if vL−cL<x
t< S∗
L
WRif x
t≥S∗
L,
(15)
with
ρfan =ρL2
γ0+ 1 +γ0−1
(γ0+ 1) cLvL−x
t2
γ0−1
,(16)
vfan =2
γ0+ 1 cL+γ0−1
2vL+x
t,(17)
pfan =pL2
γ0+ 1 +γ0−1
(γ0+ 1) cLvL−x
t2γ0
γ0−1
.(18)
Figure 2 shows the solution structure of the vacuum Riemann problem.
The expansion of the continuum to the vacuum may be described as a
free surface problem. For many physical problems, non-equilibrium effects
have to be taken into account to model the physical phenomena correctly. In
this case, a gas kinetic description has to be adopted. For our purpose here,
the gas vacuum problem is an example for a single rarefaction wave, that can
also be handled by molecular dynamics in an appropriate way.
12
x
t
UR= (0,0,0)>
UL
S∗
L
Figure 2: Solution structure of the vacuum Riemann problem, where the vacuum is lo-
cated on the right. The left state is connected to the vacuum by a rarefaction
wave, the contact discontinuity attaches the boundary of the rarefaction wave
adjacent to the vacuum.
3. Numerical Methods
3.1. Fluid-Solver FLEXI
The conservation equations of the fluid equations are solved by a discon-
tinuous Galerkin spectral element method (DGSEM), which is implemented
as open source code FLEXI 1. The numerical method is described by Hin-
denlang et al. [16] in detail. Hence, we only survey the basic building blocks
to get some impression of the numerical method used. As usual in the dis-
continuous Galerkin approach, the solution and the fluxes are approximated
by polynomials in each grid element allowing discontinuities between the el-
ements. In DGSEM, the polynomial basis is a nodal one with Lagrangian
polynomials defined by the Gauss points. We project the physical grid cell
1https://www.flexi-project.org/
13
to a reference element and derive the weak formulation of the conservation
equations in this reference element. The resulting volume and surface inte-
grals are approximated by Gauss quadratures. As for Finite Volume schemes,
numerical flux functions are needed for the coupling to the neighboring ele-
ments. If discontinuities appear in the solution, such as shock waves or mate-
rial boundaries, the method is supplemented with a shock capturing method
[17, 18], in which in elements with oscillations a total variation diminishing
finite volume method on sub-cells is activated. The indicator of Persson and
Peraire [19] is used to detect oscillations. The structure of the DGSEM is
kept - the local finite volume scheme on sub-cells may be interpreted as the
alternative to evaluate the volume integral. The time integration method is
explicit with high order Runge-Kutta schemes as proposed in [20].
3.2. Numerical Flux Calculation
The numerical flux solver for the vacuum Riemann problem is a little bit
subtle not to produce negative internal energies and a break-up of the simu-
lation. We apply here the Suliciu relaxation Riemann solver [21, 22] for the
isentropic Euler equations which handles the vacuum state very well. The
extension to the complete equations of gas dynamics has been described by
Bouchut [23], based on a relaxation system for the Euler equations using a va-
riety of supplementary variables such as pressure, speed of sound and entropy
amongst others. If the underlying EOS is convex, an approximate Riemann
solver can be constructed to solve the relaxation system and, by extension,
the original Euler equations. To use a general EOS, the original method has
14
been modified for this study by defining a local isentropic coefficient via
γloc =p
ρ + 1.(19)
In case of a perfect gas, the local isentropic coefficient recovers the exact one:
γloc =γ0=c0
p
/c0
V=p
/ρ + 1. As shown in [23], the Suliciu solver preserves
positivity of ρand and guarantees the inequalities of entropy as well as
the maximum principle of entropy. We applied the numerical flux in flux
difference form as proposed by Shen et al. [24]
F=1
2(FL+FR)−1
2(|SL|(U∗
L−UL) + |S∗|(U∗
R−U∗
L) + |SR|(UR−U∗
R)) ,
(20)
in which the vacuum state can be directly inserted and the relaxation speeds
remain bounded.
3.3. Molecular Dynamics Simulations
MD simulations were carried out with the open source code ls1 mardyn
[25] that numerically solves Newton’s equation of motion for large particle
ensembles. Molecular interactions were assumed to be pairwise additive and
evaluated explicitly within a specified cutoff radius rc. In the present study,
the LJTS fluid was simulated, being fully consistent with the EOS by Thol
et al. [5] and Heier et al. [6]. This model has only two state-independent
parameters for a given molecular species, i.e. for size σand dispersive energy
, and is very well suited to describe the thermodynamic properties of simple
fluids, like the noble gases and methane [26, 27]. It was selected since it of-
fers additional advantages. First, the cutoff radius rc= 2.5σis comparatively
small, driving down computational costs. Second, long-range corrections can
15
be omitted entirely because of the truncated potential definition. This is par-
ticularly convenient for inhomogeneous systems, where long range correction
techniques are not straightforward.
4. Results
4.1. Supercritical Shocktube Problem
A shock tube scenario was defined for the LJTS fluid in the supercritical
regime. The initial states were chosen such that one state lies in the liquid-like
region and the other in the vapour-like region. A high temperature ensures
that the solution remains in the supercritical regime. The initial data under
consideration were
(ρ, v, T)>=
(0.6,0.0,3.0)>if x<x0
(0.2,0.0,3.0)>if x≥x0.
(21)
For both equations of state used in this study, the initial pressures are calcu-
lated from the density and temperature to avoid errors by conversion com-
paring with the MD simulations , hence they differ in their values. E.g., for
the LJTS EOS, the initial left pressure was p= 3.99 and for the PeTS EOS
it was p= 3.80 while the right pressure was p= 0.657 and p= 0.648, re-
spectively. This caused a shift in the solution for the pressure but achieved
an improvement of the comparison of temperatures between MD and CFD
data.
All physical properties were non-dimensionalized in terms of reference
length σref = 1 ˚
A, reference energy ref /kB= 1 K and reference mass
mref = 1 u. Consequently, the time reference is tref =σref pmref /ref . In these
16
reduced units, the critical point of the LJTS fluid is given by ρc= 0.319 and
Tc= 1.086. The spatial domain was chosen as [0,4000] and the initial discon-
tinuity was located at the midpoint x0= 2000. The problem was simulated
until tend = 300.91. For the CFD code FLEXI , the one-dimensional domain
was discretized into 100 elements and a polynomial degree of N= 3 was cho-
sen that should result in a fourth order approximation. Discontinuities were
detected with the Persson indicator [19] with threshold Indupper =−3. In this
case, the grid element was handled by switching to the local TVD finite vol-
ume scheme on the subgrid. The DG was applied again for Indlower =−5. For
time integration we used the explicit low storage fourth order Runge-Kutta
scheme with five stages [20]. The stability criterion was chosen as CFL = 0.9.
In addition to the numerical simulation, an exact solution of the Riemann
problem was calculated by the solution procedure of Kamm [12].
For the corresponding MD simulation, two homogeneous phases with a
cuboid geometry were equilibrated at a temperature T= 3.0, one with the
higher density ρ1= 0.6 and one with the lower density ρ2= 0.2. They
were merged such that physical contact between the phases was established
through two planar interfaces, yielding a symmetric system, cf. Fig. 3. A few
individual particles overlapping across at the interface were discarded. Peri-
odic boundary conditions were established in all directions to mimic infinite
extended interfaces and to avoid boundary effects in xdirection.
For symmetry reasons, the origin of the spatial coordinate xwas de-
fined to coincide with the geometric center of the system so that only the
right half of the system was considered in the following. Initially, the phases
had an equal thickness of ∆x= 2000 and the initial interface position was
17
Figure 3: Snapshot of the shock tube scenario in a diagonal view to give an impression of
the system on particle scale, and a plan view to clarify the system dimensions.
The system was symmetric around x= 0. Particles of the higher and lower
density phase ρ1= 0.6 and ρ2= 0.2 were colored in red and green, respectively.
located at x0= 2000. The cross-sectional area was comparatively large
Ayz = ∆y∆z= 1002to obtain good statistics while investigating a rapid
process. Consequently, the system consisted of a large number of particles
N= 3.2·107. Starting the simulation from the initial configuration, all in-
dividual particle trajectories were followed without any interventions, i.e.,
thermostating was avoided. Newton’s equations of motion were solved nu-
merically, employing the Leapfrog integrator that achieves a good energy
conservation. Maintaining these conditions, the MD simulation results are
directly comparable to the solution of the Riemann problem. To follow the
temporal evolution of spatially resolved quantities, one-dimensional sampling
18
was conducted to obtain density, temperature and hydrodynamic velocity
profiles. For this purpose, the system was divided into bins with a thickness
of δ= 0.25 and a volume of ∆V=δAyz = 2500. The local density ρin a bin
was calculated by
ρ=N
∆V,(22)
where Nis the number of particles in the bin. The velocity vector of the j-th
particle is named vjand has the three components vx,j, vy,j and vz,j into x- ,
y- , and z-direction, respectively. The hydrodynamic velocity in x-direction
ˆvxwas obtained from the arithmetic mean of the velocity component vxof
all particles in the bin
ˆvx=1
N
N
X
j=1
vx,j.(23)
The temperature was sampled for the three spatial directions individually
Tx,Tyand Tzbecause it cannot unconditionally be treated as a single scalar
quantity under strong non-equilibrium conditions. The temperature is a mea-
sure for the kinetic energy of thermal motion only. Therefore, when accumu-
lating the kinetic energy of the particle collective in a bin, the contribution
due to hydrodynamic velocity has to be subtracted, yielding for the temper-
ature in x- direction
kBTx=m
3N
N
X
j=1
(vx,j −ˆvx)2.(24)
For the other two spatial directions, a vanishing hydrodynamic velocity can
be assumed ˆvy= ˆvz= 0, hence,
kBTy=m
3N
N
X
j=1
v2
y,j,(25)
19
applies, where Tzis analogous. Since Tyand Tzwere identical in the present
scenario, they were summarized to a mean quantity
Tyz = (Ty+Tz)/2.(26)
Quantities were averaged over a period of 5000 time steps, where a single time
step ∆tcorresponds to a physical time of 2 fs for a fluid like Argon. This
was a compromise for the sampling period being long enough to obtain good
statistics, but short enough to ensure that the rapid changing profiles did not
get blurred. The simulation was terminated when the system’s reaction on
the discontinuity was propagated to the system boundary, which was already
the case after 2 ·105time steps.
The results are shown in Fig. 4 for density, temperature, velocity and
pressure. Each quantity was compared to MD simulation data, except for
the pressure, which was not sampled in the MD simulations. The numerical
results indicate that all properties are in excellent agreement with each other.
The wave phenomena and the constant states in between are clearly iden-
tified: A shock wave propagates to the right, a rarefaction wave propagates
to the left and a contact discontinuity is located in between. The approx-
imate CFD solution adds some numerical diffusion, which smears out the
shock and contact discontinuity compared to the exact solution. However,
diffusion is also visible in the MD data, where viscous and diffusive effects
are fully resolved. The shock capturing method was, as expected, only active
around the shock wave, as indicated by the grey rectangle. Here, the method
switches locally to a second order FV method with TVD reconstruction of
the gradients on the sub-cells within the coarse DG grid element. Outside of
this region, the solution was produced by the high order DG method.
20
Density
0.2
0.4
0.6
Shock
Capturing
ρ[−]
MD LJTS exact PeTS exact LJTS DG PeTS DG
Temperature
2
3
4
T[−]
Velocity
0 1,000 2,000 3,000 4,000
0
0.5
1
x[−]
v[−]
Pressure
0 1,000 2,000 3,000 4,000
0
2
4
x[−]
p[−]
Figure 4: Results for the supercritical shocktube problem. The plot shows the results
of the MD simulations, the exact solution of the corresponding macroscopic
Riemann problem and the approximate solution by the DGSEM scheme using
the Suliciu relaxation solver. The Riemann problem was solved resting either
on the LJTS EOS or the PeTS EOS. The shock capturing was active in the
region indicated by the grey rectangle.
Visible differences were observed for velocity and pressure in the com-
parison of the different EOS that were underlying to the CFD simulations.
At the rarefaction wave, the PeTS EOS produces a slightly shifted position,
whereas the LJTS EOS lies directly on the MD data. This effect can be at-
21
tributed to an underestimation of the speed of sound by the PeTS EOS in
the high density region as well as the different pressures of the initial datum.
Similarly, the PeTS EOS underestimates the pressure in the high density re-
gion. Note that the pressure difference of the initial liquid states stems from
the choice of initial values. To find the best basis for the comparison with
MD simulation data, the temperature was chosen to be equal in the initial
datum which caused the different plateaus in pressure.
Our conclusion with respect to this comparison of the numerical results
is as follows. By using these corresponding thermodynamic relations of the
LJTS EOS and the LJTS flow assumption in the MD simulations, we are
able to get the same results by CFD and MD.
4.2. Vacuum Expansion Problem
As a first step to more subtle cases, we choose the so called vacuum
Riemann problem, in which one state is the vacuum. The solution consists
of a rarefaction wave to the vacuum with the attached contact discontinuity.
The vacuum Riemann problem can be handled well by MD simulations. For
the continuum equations, the transition to vacuum is a limit to the range of
the validity of the equations and the approximation needs a flux calculation
that keeps density and internal energy non-negative.
The vacuum expansion problem was considered for the LJTS fluid again
and with the continuum phase on the left and the vacuum on the right.
We investigate in the following four cases with a common temperature of
TL= 3.0. The strength of the rarefaction wave was prescribed by a varying
density. The baseline test case VP1 had an initial density ρL= 0.05. The
other cases had the densities: ρL= 0.01 (VP2), ρL= 0.1 (VP3) and ρL= 0.2
22
(VP4). The full initial data of the vacuum Riemann initial were
(ρ, v, T)>(x, 0) =
(ρL,0.0, TL)>if x<x0
(0.0,0.0,0.0)>if x≥x0.
(27)
The computational domain was chosen to be x∈[0,4000] and the initial
discontinuity was located at x0= 2000. In DG method we switched for
these problems to the second order FV method with TVD limited gradient
reconstruction, which is used as shock-capturing in sub-cells, to guarantee
the positivity near the vacuum. The approximate Suliciu relaxation solver
was used for the numerical flux function and time integration was done via
an Euler explicit time integration with CFL = 0.1. At the vacuum domain
boundary, the boundary values were extrapolated from the flow region as
in a supersonic case. The domain was discretized into 400 elements. Both
the LJTS and PeTS EOS, were alternatively used to obtain approximate
solutions of the Riemann problem. Additionally, the exact solution for a
perfect gas with γ0= 1.667 was computed, which corresponds to the perfect
gas contribution of the LJTS model fluid.
For the MD simulations, exactly the same system geometry was used.
Moreover, the simulation procedure and the sampling were conducted in the
same way as described for the shock tube scenario with the left values of the
density: ρL= 0.05,0.01,0.10,0.2 and the replacement of the lower density
phase by free space ρR= 0 (vacuum). To maintain instant vacuum conditions,
all particles that reached the system boundary at x= 4000 were discarded,
i.e. they left the computational domain.
23
4.2.1. Time Evolution of VP1
Solutions of the time evolution of density, velocity and temperature for
testcase VP1 at three times t= 91.18, t= 182.37 and t= 273.55 are shown
in Fig. 5. The structure of the exact solution is reproduced by our CFD
t=91.18
Density
10−5
10−3
10−1
Knρ<0.02
ρ[−]
Temperature
0
2
4
KnT<0.02
T[−]
MD MD TxMD Tyz Exact LJTS EOS PeTS EOS
Velocity
0 1,000 2,000 3,000 4,000
0
5
10
15
20
x[−]
v[−]
t=182.37
Knρ<0.02
KnT<0.02
0 1,000 2,000 3,000 4,000
x[−]
t=273.55
Knρ<0.02
KnT<0.02
0 1,000 2,000 3,000 4,000
x[−]
Figure 5: Results for the vacuum Riemann problem VP1 at time instances t= 91.18,
t= 182.37 and t= 273.55. The plot shows the approximate CFD solutions
resting either on the LJTS or the PeTS EOS, an exact solution of the Riemann
problem for the perfect gas law with γ0= 1.667 and MD simulation results.
24
approximation as well as by MD simulation. The solution is determined by
the rarefaction wave, which expands the gas isentropically to the vacuum.
Note that the contact wave coincides with the tail of the rarefaction wave and
is therefore not visible. At all time instances, depicted in Fig. 5, the agreement
between all methods is excellent as long as the density is large enough. The
LJTS and PeTS EOS show almost perfect agreement with each other for
every variable, which is attributed to the identical perfect gas contribution
that dominates the fluid behavior at low density.
The MD data were used to calculate the local Knudsen number to indicate
the continuum breakdown. The Knudsen number is the ratio of the mean free
path λand a characteristic length Lof the problem:
Kn = λ
L.(28)
In the present case, the choice of the characteristic length is not as clear.
Therefore, we consider the local Knudsen number introduced by Boyd et al.
[28]. It relies on the gradient of a physical state variable f
W, usually the
density f
W=ρor temperature f
W=T, and reads as
Knf
W=λ|∇f
W|
f
W.(29)
An approximation of the mean free path λwas given for a monatomic perfect
gas in Bird [29] as
λ=1
√2πρσ2,(30)
where σis the molecular diameter and ρthe number density. Following Burt
and Boyd [30], a flow can be considered to be a continuum if Knf
W<0.02.
To calculate the local Knudsen number, the MD data were smoothed using
25
a moving average with a window size of 100. We considered the density-
based local Knudsen number Knρand the temperature-based local Knud-
sen number KnT. The value Knρ,T <0.02 for the density-based and the
temperature-based Knudsen number is plotted as a line parallel to the x
axes in the density and temperature plots, respectively. Both definitions of
the local Knudsen number lead to similar values that indicate the location,
at which the continuum assumption starts to fail.
Notable differences of the results for the different methods occur beyond
the cut-off Knudsen number only. The density of all the methods approach
zero. But the locations differ. The expansion to the vacuum by the DG scheme
is the fastest one, followed by the MD results. The expansion of the exact
perfect gas solution is the slowest one. This behavior is manifested also in
the velocity plots that show strong deviations in the low density region. All
results have the same gradient, but the front speeds differ. Such an artefact
has also been observed by Bouchut [23]. The largest deviation was found
for the temperature. While the exact perfect gas solution tends to a zero
temperature, the DG scheme and the MD method do not. For the MD simu-
lation the expected temperature is zero, as the number of particles strongly
decrease towards the vacuum, which also decreases the interaction probabil-
ity between particles. Consequently, a stronger deviation from the Maxwell
distribution of the particles appears and the definition of a continuum tem-
perature breaks down into an oscillatory behavior. Figure 5 also shows the
temperatures that were calculated from the thermal motion of the particles
in each spatial direction individually, cf. Eqs. (24) for Txand (26) for Tyz.
26
Differences between both temperatures are visible which clearly indicate that
a continuum temperature can not longer be defined.
4.2.2. Convergence for VP1
Because of the subtle behavior of the temperature approximation in the
low density region, a convergence study for the DG scheme has been per-
formed for the VP1 problem. VP1 was simulated based on the Suliciu nu-
merical flux and the EOS of perfect gas with γ0= 1.667 on different grids
ranging from 100 elements to 6400 elements. The left panel of Fig. 6 shows
results for the temperature profile. Even with the EOS of perfect gas, a simi-
lar increase in temperature can be observed, but the temperature approaches
zero as the mesh spacing is reduced. The convergence rate of the local maxi-
mum of the temperature in the higher Knudsen region is plotted in the right
panel of Fig. 6. For an increasing number of elements, convergence is attained
with a low rate of approximately 0.6. The temperature increase obtained by
the Suliciu flux is caused by different slopes of pressure and density. Consid-
ering the perfect gas law, the temperature was calculated by
T=p
ρkB
.(31)
Near the vacuum, for very small values of density and pressure, the slopes
of density and pressure are not identical. Beyond some threshold, the den-
sity gradient is smaller than the pressure gradient which, following Eq. (31),
caused the non-physical temperature increase. Refining the grid this artefact
disappears. This also appears in the simulations resting on the LJTS EOS
and the PeTS EOS. In the high Knudsen regime, both EOS fall back to a
perfect gas formulation. However, the implementation for both EOS still eval-
27
uates Eq. (4) with a remaining contribution of the residual Helmholtz energy
instead of Eq. (31). Consequently, the peak for the increasing temperature is
reproduced differently than for the classical perfect gas law.
Temperature
0 1,000 2,000 3,000 4,000
0
1
2
3
x[−]
T[−]
Exact 100 Elements 800 Elements 6400 Elements
Convergence Rate
102102.5103103.5
10−1
100
-0.6
1
Elements [−]
Tmax [−]
Figure 6: CFD convergence for the vacuum Riemann problem VP1 at the time instance
t= 45.59. The approximate solutions were obtained for the perfect gas with
γ0= 1.667. Three grids are shown, constituted by 100, 800 and 6400 elements.
The data are compared to the exact solution.
Fig. 7 shows a comparison between the temperature profiles of two MD
simulations with a different number of particles. The coarser simulation pro-
duced more scattering, but both simulations agree very well for small Knud-
sen numbers. In the high Knudsen number regime, the fine grid simulation
shows a similar increase of temperature as the coarse grid simulation. Fur-
thermore, the position of the local minimum is identical. Our explanation of
this non-physical increase of the temperature is that non-equilibrium states
affects again the averaging procedure for the mean temperature.
28
Temperature
0 1,000 2,000 3,000 4,000
0
1
2
3
x[−]
T[−]
MD coarse MD fine
Figure 7: MD simulation of the vacuum Riemann problem VP1 at time t= 45.59 for two
different numbers of particles.
4.2.3. Results for VP2 and VP3
For the next two vacuum Riemann problems only results at time t=
273.55 are shown. VP2 has the initial density ρL= 0.01. The results are
shown in Fig. 8, VP3 hast the initial density ρL= 0.1. The results are shown
in Fig. 9. The qualitative behavior is similar to the results for VP1. In the
case of VP2 with the smaller density, the noise in the MD data increased,
which is caused by the smaller number of particles per volume.
4.2.4. Results for VP4
The vacuum Riemann problem VP4 has a larger density of ρ= 0.2. For
this case, the expansion is much stronger than in the previous cases and leads
to a sooner crossing of the isentropic expansion and the vapour binodal curve.
Figure 10 shows the loci of constant entropy for all vacuum Riemann prob-
29
Density
0 1,000 2,000 3,000 4,000
10−6
10−5
10−4
10−3
10−2
10−1
x[−]
ρ[−]
MD Exact LJTS EOS PeTS EOS
Velocity
0 1,000 2,000 3,000 4,000
0
5
10
x[−]
v[−]
Temperature
0 1,000 2,000 3,000 4,000
0
1
2
3
x[−]
T[−]
Figure 8: Results for the vacuum Riemann problem VP2 at time t= 273.55. Shown is
the CFD solution based on the LJTS EOS, the exact solution for the perfect
gas EOS with γ= 1.667 and the MD simulation results.
Density
0 1,000 2,000 3,000 4,000
10−5
10−3
10−1
x[−]
ρ[−]
MD Exact LJTS EOS PeTS EOS
Velocity
0 1,000 2,000 3,000 4,000
0
5
10
x[−]
v[−]
Temperature
0 1,000 2,000 3,000 4,000
0
1
2
3
x[−]
T[−]
Figure 9: Results for the vacuum Riemann problem VP3 at time t= 273.55. Shown is
the CFD solution based on the LJTS EOS, the exact solution for the perfect
gas EOS with γ= 1.667 and the MD simulation results.
lems in state space. Eventually, most of the curves cross the vapour binodal,
but only for the largest density droplets appeared. The results for density,
velocity and temperature are presented in Fig. 11 at time t= 300.91. Sim-
ilar to the other scenarios, the results of the different methods are in very
good agreement. The LJTS and PeTS EOS allow a robust simulation in this
regime, in which small droplets in the expansion fan occur by accounting for
30
Isentrope
100101102103
10−3
10−2
10−1
100
ρ= 0.01
ρ= 0.05
ρ= 0.1
ρ= 0.2
1
ρ[−]
p[−]
Binodal
Isentrope
Isothermal for T= 3.0
Figure 10: Isentropic expansion curves in state space for the vacuum Riemann problems
VP1-VP4.
a undercooled vapour. Due to the fact that the phase change is not explicitly
included into the simulation and the equation of state, the CFD simulations
break down if the initial density is further increased. In the MD simulations,
the droplets grew in size. In this case the strong isentropic expansion moves
into the region of the LJTS fluid, where both fail to provide macroscopic
physical states. This is the region, in which the underlying pressure function
becomes non-convex and the speeds of sound become imaginary. In these
cases, a phase transition model has to be introduced and the multiphase Rie-
mann problem changes its structure taking into account the phase transition
line.
31
Density
0 1,000 2,000 3,000 4,000
10−5
10−3
10−1
x[−]
ρ[−]
MD Exact LJTS EOS PeTS EOS
Velocity
0 1,000 2,000 3,000 4,000
0
2
4
6
x[−]
v[−]
Temperature
0 1,000 2,000 3,000 4,000
0
1
2
3
x[−]
T[−]
Figure 11: Results for the vacuum Riemann problem VP4 at time t= 300.91. Shown is
the CFD solution based on the LJTS EOS, the exact solution for the perfect
gas EOS with γ= 1.667 and the MD simulation results.
5. Conclusion
We compared the macroscopic Riemann solution and approximate solu-
tions, obtained by a CFD solver, with microscopic solutions, obtained by
molecular dynamic simulations. We picked up two Riemann problem scenar-
ios. The numerical results can be directly compared because both methods
are based on a Lennard-Jones model fluid with truncated and shifted po-
tential. The MD simulations use this potential, for the CFD simulation the
corresponding accurate equations of state were used, proposed in Thol et al.
[5] and in Heier et al. [6]. The CFD simulations were based on a hybrid
DG/FV method with the Suliciu relaxation solver as numerical flux func-
tion. The two scenarios under investigation were a super-critical shock tube
problem and a strong expansion of a gas into vacuum. The agreement be-
tween the macroscopic and microscopic solutions is nearly perfect. In the
super-critical regime, the underlying EOS were able to reproduce liquid-like
and vapour-like states and their transition from one to the other regime.
32
In the second problem, the vacuum Riemann problem, the EOS ap-
proaches the regime of a perfect gas. Hence, the exact solution for perfect gas
law coincides very well with the results. Four test cases with varying initial
densities have been considered. The agreement with the microscopic solu-
tions was very good. Beyond the critical Knudsen number, the microscopic
solution developed non-equilibrium states, where a continuum temperature
is not well-defined any more. Hence, in the low density region, the macro-
scopic and microscopic solutions differed. By increasing the strength of the
expansion in the vacuum Riemann problem, condensation starts and small
liquid droplets appeared in the expansion fan simulated by molecular dy-
namics. The real gas EOS could reproduce the solution as long as this effect
could be approximated as undercooled vapour phase. A further increase of
the strength of the expansion wave lead to a failure and a break down of
the continuum simulations. Hence, for stronger rarefaction waves, the phase
transition has to be explicitly modelled in the numerical approach.
These results strongly motivate the use of MD simulations as comparison
in the cases with phase transitions. The MD results may establish thermo-
dynamic consistent solutions for benchmark problems that can be used as
benchmarks for macroscopic multiphase flow approximations. As mentioned
in the introduction, our research activities are embedded in the DFG col-
laborative research center SFB-TRR 75 ”Droplets under extreme ambient
conditions” and we are able to continue this research in the forthcoming
years. MD results for these situations with phase transitions will be well
documented.
33
Acknowledgements
The work was supported by the German Research Foundation (DFG)
through SFB-TRR 75 “Droplet Dynamics under Extreme Ambient Con-
ditions”. The MD simulations were performed on the national supercom-
puter Cray XC40 (Hazel Hen) at the High Performance Computing Center
Stuttgart (HLRS).
References
[1] S. K. Godunov, A difference method for numerical calculation of discon-
tinuous solutions of the equations of hydrodynamics, Mat. Sb. (N.S.) 47
(1959) 271–306.
[2] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a
practical introduction, 2 ed., Springer Science & Business Media, Berlin,
Heidelberg, 2009.
[3] R. Menikoff, B. J. Plohr, The Riemann problem for fluid flow of real
materials, Reviews of Modern Physics 61 (1989) 75–130.
[4] P. G. LeFloch, Hyperbolic Systems of Conservation Laws, Lectures in
Mathematics. ETH Z¨urich, 1 ed., Birkh¨auser, Basel, 2002.
[5] M. Thol, G. Rutkai, R. Span, J. Vrabec, R. Lustig, Equation of State
for the Lennard-Jones Truncated and Shifted Model Fluid, Interna-
tional Journal of Thermophysics 36 (2015) 25–43. doi:10.1007/s10765-
014-1764-4.
34
[6] M. Heier, S. Stephan, J. Liu, W. G. Chapman, H. Hasse, K. Langen-
bach, Equation of state for the Lennard-Jones truncated and shifted
fluid with a cut-off radius of 2.5 sigma based on perturbation theory
and its applications to interfacial thermodynamics, Molecular Physics
116 (2018) 2083–2094.
[7] O. Kunz, W. Wagner, The GERG-2008 wide-range equation of state for
natural gases and other mixtures: an expansion of GERG-2004, Journal
of Chemical & Engineering Data 57 (2012) 3032–3091.
[8] O. Kunz, R. Klimeck, W. Wagner, M. Jaeschke, The GERG-2004 wide-
range equation of state for natural gases and other mixtures, Ber. VDI,
VDI-Verlag, Duesseldorf (2007).
[9] R. Span, Multiparameter Equations of State, Springer-Verlag, Berlin,
Heidelberg, 2000.
[10] J. A. Barker, D. Henderson, Perturbation Theory and Equation of State
for Fluids: The SquareWell Potential, The Journal of Chemical Physics
47 (1967) 2856–2861. doi:10.1063/1.1712308.
[11] P. Colella, H. M. Glaz, Efficient solution algorithms for the Riemann
problem for real gases, Journal of Computational Physics 59 (1985)
264–289. doi:10.1016/0021-9991(85)90146-9.
[12] J. R. Kamm, An Exact, Compressible One-Dimensional Riemann Solver
for General, Convex Equations of State, Technical Report LA-UR-15-
21616, Los Alamos National Laboratory, 2015.
35
[13] E. Halter, E. Martensen, A fast solver for Riemann problems,
Mathematical Methods in the Applied Sciences 7 (1985) 101–107.
doi:10.1002/mma.1670070107.
[14] T.-P. Liu, J. Smoller, On the vacuum state for the isentropic gas dy-
namics equations, Advances in Applied Mathematics 1 (1980) 345–359.
doi:10.1016/0196-8858(80)90016-0.
[15] C.-D. Munz, A tracking method for gas flow into vacuum based on
the vacuum Riemann problem, Mathematical Methods in the Applied
Sciences 17 (1994) 597–612. doi:10.1002/mma.1670170803.
[16] F. Hindenlang, G. J. Gassner, C. Altmann, A. Beck, M. Stau-
denmaier, C.-D. Munz, Explicit discontinuous Galerkin methods
for unsteady problems, Computers & Fluids 61 (2012) 86–93.
doi:10.1016/j.compfluid.2012.03.006.
[17] M. Sonntag, C.-D. Munz, Shock Capturing for Discontinuous Galerkin
Methods using Finite Volume Subcells, in: J. Fuhrmann, M. Ohlberger,
C. Rohde (Eds.), Finite Volumes for Complex Applications VII-Elliptic,
Parabolic and Hyperbolic Problems, Springer International Publishing,
2014, pp. 945–953.
[18] M. Sonntag, C.-D. Munz, Efficient Parallelization of a Shock Cap-
turing for Discontinuous Galerkin Methods using Finite Volume
Sub-cells, Journal of Scientific Computing 70 (2017) 1262–1289.
doi:10.1007/s10915-016-0287-5.
36
[19] P.-O. Persson, J. Peraire, Sub-Cell Shock Capturing for Discontinu-
ous Galerkin Methods, in: Proceedings of the 44th AIAA Aerospace
Sciences Meeting and Exhibit, American Institute of Aeronautics and
Astronautics, Reno, 2006.
[20] C. A. Kennedy, M. H. Carpenter, R. Lewis, Low-storage, explicit
Runge-Kutta schemes for the compressible Navier-Stokes equations, Ap-
plied Numerical Mathematics 35 (2000) 177–219. doi:10.1016/S0168-
9274(99)00141-5.
[21] I. Suliciu, On modelling phase transitions by means of rate-type consti-
tutive equations. Shock wave structure, International Journal of Engi-
neering Science 28 (1990) 829–841. doi:10.1016/0020-7225(90)90028-H.
[22] I. Suliciu, Some stability-instability problems in phase transitions
modelled by piecewise linear elastic or viscoelastic constitutive equa-
tions, International Journal of Engineering Science 30 (1992) 483–494.
doi:10.1016/0020-7225(92)90039-J.
[23] F. Bouchut, Nonlinear Stability of Finite Volume Methods for Hyper-
bolic Conservation Laws, Frontiers in Mathematics, 1 ed., Birkh¨auser,
Basel, 2004.
[24] Z. Shen, W. Yan, G. Yuan, A robust and contact resolving Riemann
solver on unstructured mesh, Part II, ALE method, Journal of Compu-
tational Physics 268 (2014) 456–484. doi:10.1016/j.jcp.2014.03.003.
[25] C. Niethammer, S. Becker, M. Bernreuther, M. Buchholz, W. Eckhardt,
A. Heinecke, S. Werth, H.-J. Bungartz, C. W. Glass, H. Hasse, J. Vrabec,
37
M. Horsch, ls1 mardyn: The Massively Parallel Molecular Dynamics
Code for Large Systems, Journal of Chemical Theory and Computation
10 (2014) 4455–4464. doi:10.1021/ct500169q.
[26] J. Vrabec, G. K. Kedia, G. Fuchs, H. Hasse, Comprehensive study
of the vapour-liquid coexistence of the truncated and shifted Lennard-
Jones fluid including planar and spherical interface properties, Molecular
Physics 104 (2006) 1509–1527. doi:10.1080/00268970600556774.
[27] G. Rutkai, M. Thol, R. Span, J. Vrabec, How well does
the Lennard-Jones potential represent the thermodynamic proper-
ties of noble gases?, Molecular Physics 115 (2017) 1104–1121.
doi:10.1080/00268976.2016.1246760.
[28] I. D. Boyd, G. Chen, G. V. Candler, Predicting failure of the continuum
fluid equations in transitional hypersonic flows, Physics of Fluids 7
(1995) 210–219. doi:10.1063/1.868720.
[29] G. Bird, Molecular Gas Dynamics and The Direct Simulation of Gas
Flow, Oxford University Press, Oxford, 1994.
[30] J. M. Burt, I. D. Boyd, A hybrid particle approach for continuum and
rarefied flow simulation, Journal of Computational Physics 228 (2009)
460–475. doi:10.1016/j.jcp.2008.09.022.
38
