An Asymptotic-Preserving Method for a Relaxation of the Navier-Stokes-Korteweg Equations

The Navier-Stokes-Korteweg (NSK) equations are a classical diffuse-interface model for compressible two-phase flow. As direct numerical simulations based on the NSK system are quite expensive and in some cases even impossible, we consider a relaxation of the NSK system, for which robust numerical methods can be designed. However, time steps for explicit numerical schemes depend on the relaxation parameter and therefore numerical simulations in the relaxation limit are very inefficient. To overcome this restriction, we propose an implicit-explicit asymptotic-preserving finite volume method. We prove that the new scheme provides a consistent discretization of the NSK system in the relaxation limit and demonstrate that it is capable of accurately and efficiently computing numerical solutions of problems with realistic density ratios and small interfacial widths.


Introduction
There are in general two approaches to describe the behavior of multi-phase fluids, the sharp interface (SI) and the diffuse interface (DI) approach.The first approach represents multiple phases with different sets of equations that are coupled by some interface conditions.The second approach, which used to describe, e.g. the merging process of droplets and bubbles, needs only one set of equations to model the phases and does not require the location of the interface to be tracked explicitly.
In this paper, we consider two DI models for a homogeneous two-phase compressible fluid: The Navier-Stokes-Korteweg system (NSK) and a relaxation system for the NSK system.The NSK system goes back to the work of Korteweg [31] and was formulated in its present form in [18,2].The NSK system uses a Van-der-Waals like pressure function to identify two distinct phases and a third-order term to model phase transitions.Many authors achieved analytical results on the well-posedness of the NSK system and its variants, e.g.[4,7,24,32,37].While some numerical methods have been successfully developed and used for these and related families of problems, see, e.g., [6], [8], [16], [21], [25], [39], there are still many open problems, for which accurate and efficient numerical methods are yet to be designed.Up to our knowledge, the robust computation of realistic density values for liquid and vapor phases has not been suggested for the NSK model.Additionally, numerical methods also fail in cases when very small interface widths close to a sharp interface are to be considered.In both cases, the occurring problems are related to steep density gradients.Another source of difficulty one comes across while numerically solving NSK systems is related to the Van-der-Waals like form of the pressure equation.The latter prevents one from using upwind hyperbolic solvers, which have been successfully applied, e.g., to stabilize computations for Navier-Stokes equations with high Reynolds numbers.
The issue of very small interfaces is especially important because the NSK model can only provide the correct amount of capillary forces if the interface is extremely small [17,25,27].One idea to loosen the strict coupling between interfacial width and capillary forces is to introduce an additional Cahn-Hilliard or Allen-Cahn type equation for a new phase field variable.This was done for example in [1,5,42].Another ansatz to avoid some of the difficulties for the NSK systems suggests to introduce a relaxation of the NSK system, in which the third-order term is replaced by a first-order term and a Poisson equation, that defines a new phase field parameter, see, e.g., [37].This model is parametrized by a so-called Korteweg parameter α.If the Korteweg parameter tends to zero, the relaxation system formally converges to the NSK system.The most important feature of the relaxation system is the fact, that the first-order part is purely hyperbolic for sufficiently small Korteweg parameter.It should also be observed that the addition of the Poisson equation to the system does not increase the computational cost of numerical simulations as test cases demonstrate that the time savings, that come from the fact that one does not have to solve a third order system, are greater than the loss of time that comes from the numerical solution of the Poisson equation.These properties can be exploited to construct robust numerical schemes for the relaxation system.In [34], it has been shown that the overall approach is robust for problems with large density ratios and small interfacial widths.However, the numerical scheme proposed in [34] is an explicit scheme and thus the time steps decrease as the Korteweg parameter α tends to zero.
It is the main purpose of this contribution to construct an asymptoticpreserving (AP) scheme [28] in the Korteweg limit, that is, a scheme for the relaxation system that provides a consistent approximation of the original NSK system as the Korteweg parameter α tends to zero.The AP approach was developed in the framework of linear transport in diffusive regimes [22,29,33] and has been applied to many different areas, e.g.fluid and diffusion limits of kinetic models, relaxation methods for hyperbolic systems and low-mach number limits for compressible flow problems; see, e.g., [11,12,14,15,20,23,30,35,36].
In [23], the time and spatial discretization of the isentropic Euler and Navier-Stokes Equations in the low Mach number limit was investigated.Inspired by this research, we construct here a scheme that captures the Korteweg limit for the relaxation system and prove the AP property of the proposed scheme.The AP property of the new scheme is achieved by splitting the relaxation system into a non-stiff nonlinear, compressible hyperbolic Navier-Stokes like system and a system that can be treated by a Poisson solver, and allows the use of time and spatial steps that are independent of the Korteweg parameter.As the result the proposed numerical scheme is very efficient for small values of the parameter α, which is a significant improvement compared to an explicit scheme from [34].
Beyond that, we expect our scheme to be asymptotic preserving in the sharp interface limit.We cannot give analytical proof to that, but we support this statement by a numerical example.
The paper is organized as follows.In Section 2, we introduce the NSK and relaxation systems and describe their main properties together with the basic thermodynamical framework.We comment on the advantages of the relaxation system and point out why it is necessary to introduce a new scheme in order to solve the relaxation system efficiently.Section 3 contains the basic outcome of this contribution.We propose the AP scheme for the relaxation system and perform an asymptotic analysis to show that the scheme transforms into a scheme for the NSK equations in the Korteweg limit.In Section 4, we demonstrate that the algorithm provides a massive improvement compared to a standard explicit scheme for a number of problems in one and two space dimensions.

The Navier-Stokes-Korteweg system
Let an open set Ω ⊆ R d , d ∈ {1, 2, 3} up to the final time T > 0 be given.
The isentropic Navier-Stokes-Korteweg (NSK) equations in arbitrary spatial dimension are given by where ρ = ρ(x, t) is the density of the fluid, v = (v 1 (x, t), ..., v d (x, t)) T ∈ R d is its momentum and p(x, t) is the pressure.Note that ε is the Reynolds number and γε 2 is the capillary number.We refer to [26,41] for a detailed explanation on the physical meaning of the scaling ε → 0 and γ = O(1).The matrix 1) stands for the viscous part of the stress tensor which is given for the viscosity coefficients ν, µ ∈ R with µ ≥ 0 and 3ν + 2µ > 0 by We augment (2.1) with the initial data and boundary conditions that correspond to a bounded box: To describe a two-phase fluid we choose the Van-der-Waals type pressure We observe that the first-order part of (2.1) is not purely hyperbolic for all density values.Indeed, consider, for instance, the one-dimensional (1-D) case, It is easy to check that the eigenvalues of the Jacobian of the first order part of (2.1) with corresponding eigenvectors (2.7) Therefore, the the first order part of (2.1) is hyperbolic if and only if ρ ∈ The lack of hyperbolicity in the first-order part of (2.1) and the presence of the third-order derivative in the momentum balance make the numerical solution of (2.1) to be a challenging task: Explicit schemes suffer from extremely small time steps while implicit discretizations lead to badly conditioned algebraic problems.The non-monotonicity of p prevents the use of most modern shockcapturing schemes.

A relaxation for the Navier-Stokes-Korteweg system
To overcome some of the shortcomings of the classical NSK system we propose a relaxation for the NSK system [37] Here α > 0 is the Korteweg parameter and c α is a new unknown, that is defined by the additional Poisson equation, and the stress tensor T ε given by (2.2).The relaxation system (2.8) is augmented with the initial conditions and the boundary conditions (2.10) The system (2.8) has a structural advantage that becomes evident when we rewrite the time-dependent equations as Again, if we consider the 1-D case for the sake of simplicity, then the eigenvalues of the Jacobian of the first-order part of (2.11) with corresponding eigenvectors (2.13) A straightforward computation (see [38]) shows that we obtain a purely hyperbolic system for where (α 1 , α 2 ) is the interval of the decreasing pressures, i.e., p (s) < 0, see Figure 2.1.The system (2.8)-(2.10)can be seen as an approximation of the classical NSK system with (ρ α , ρ α v α , c α ) → (ρ, ρv, ρ) for the Korteweg limit α → 0, where (ρ, ρv) is the solution of the corresponding initial boundary value problem (2.1), (2.3), (2.4).We refer to [8,9,13,19,21] for first rigorous results on the Korteweg limit.It is possible to show that (2.8) formally converges to (2.1).We take the asymptotic expansion for small α and look at the balances within the equations of system (2.8).Therefore we compute the Taylor expansion at ρ (0) for the pressure A short computation leads to the following terms for the different powers of α:
We illustrate this result by the following numerical experiment that is taken from [34].
Example 2.1 (Numerical verification of the Korteweg limit).
We consider the 1-D system (2.8) with ε = 0.01 on interval Ω = (−1, 2) subject to the boundary conditions (2.10) and the following initial data From the physical point of view, these initial conditions describe two vapor bubbles surrounded by liquid fluid.A numerical solution of this initial-boundary value problem, denoted by T and a numerical solution of the corresponding classical NSK system, denoted by u h = (ρ h , ρ h v h ) T .Both solutions were computed in [34] using an explicit local discontinuous Galerkin (LDG) method, see, e.g., [3,10,16].
The distance decreases as α does.The fifth line contains the experimental order of convergence (EOC) with respect to α and the third line contains the CPU time.
The numerical results, presented in Table 2.1, indicate that the relaxed model is an O(α 2 )-approximation of the original system.As one can also see from this Table , for decreasing values of α the CPU-time is increasing due to the dependence on α of the eigenvalues (2.12).The maximum wave speed for system (2.8) is λ max = |v α max | + p (ρ α max ) (where ρ α max and v α max are the maximum values of the density and velocity, respectively).For an explicit scheme one for some 0 ≤ c DG < 1 to satisfy the CFL condition for stability.For small α one needs ∆t = O(α∆x).This restriction is a huge drawback for numerical simulations and we are interested in finding a way to circumvent this restriction.

An Asymptotic-Preserving Scheme in the Korteweg Limit
Having in mind the shortcomings of the relaxation system (2.8), we propose a new AP numerical scheme.Note that from here on we suppress the index α of all of the primal variables in order to shorten the notation.

A hyperbolic splitting
The numerical solution of the relaxation system (2.8) requires a resolution of two scales: the (fast) wave scale, coming from the Korteweg part and the (slow) convection scale.In order to obtain an accurate and efficient numerical scheme, which is able to handle both scales, we implement a splitting approach.
To this end we first introduce a parameter a and rewrite the system (2.8) in the following form (by adding and subtracting aρ∇ρ): Here the pressure is defined as where (α 1 , α 2 ) is the interval of the decreasing pressures, i.e., p (s) < 0, see To resolve the different wave length scales, we then split (3.1) into two systems: The slow dynamics are described by the system and for the fast dynamics we have The splitting parameter 0 < β < 1 in (3.4) and (3.5) determines how much of the momentum is seen by each system and the choice of a in (3.3) is motivated by hyperbolicity.In the 1-D case, for instance, the wave speeds are and when we choose a as in (3.3), we ensure that the slow system (3.4) is hyperbolic and that the wave speeds do not depend on α any more.
In what follows, we present numerical methods used to solve each one of the subsystems.We can discretize the slow system (3.3) using any explicit shockcapturing scheme for the compressible Navier-Stokes equations noting that the wave speeds are no longer stiff which avoids time step problems seen in the original system (3.1).A stability analysis shows that we need to fulfil the following CFL condition for the slow system: for some 0 ≤ c AP < 1.For the fast dynamics governed by (3.5), we introduce an implicit scheme that corresponds to the discretization of a linear hyperbolic system with variable coefficients and therefore does not have any time step restrictions.

Time discretizations of the split schemes
In this section, we outline the first-order in time discretizations for the systems (3.4) and (3.5).Note that the MUSCL methodology allows us to turn any first order scheme into a second order one.As mentioned in the previous section, we choose an explicit time discretization for the slow dynamics and the following implicit-explicit discretization for the fast dynamics: The time discretization of the nonlinear terms in the fast dynamics system (3.8) is one of the key ingredients of our scheme.With this choice, the system for the fast dynamics is a linear hyperbolic system with constant coefficients at each time step.We show the benefit of this choice further below.Now we follow the idea of [23] to obtain an efficient numerical solution strategy.We sum up (3.7) and (3.8) to have and rewrite the third equation in (3.9) as We then substitute (3.10) into the momentum equation of (3.9) to get (3.11) Next, we solve (3.11) for ρ n+1 v n+1 and obtain (3.12) Finally, we substitute (3.12) into the density equation of (3.9), obtaining a plate equation for c n+1 : The last equation is now a plate equation with respect to the unknown variable c n+1 and the terms from the previous steps can be pushed to the right hand side as source terms.The crucial part is the discretization of the non-conservative terms.With the discretization we propose, we have to solve a symmetric, sparse linear system for c n+1 .In two space dimensions, the biharmonic operator is a tridiagonal matrix M ∈ R M ×N with bandwidth O(N ), where N (M ) is the number of grid cells in the x-(y-)direction.This system can for example be solved with a conjugate gradient method.Note that we do not compute the inverse matrix M −1 , but solve the linear system directly.As long as N, M are not too large, this strategy showed to be very efficient in our numerical test cases.The updated density ρ n+1 and momentum ρ n+1 v n+1 are obtained from (3.10), (3.12).
Before we depict the spatial discretization, we want to sum up the key ideas of our numerical scheme.We split up the relaxation system into two smaller system to resolve the two wave scales.With our choice of the splitting parameter a, the time steps of the explicit scheme for the slow system are independent of α.Additionally, we proposed a linearization of the fast system, that leads to an implicit scheme.After some small computations we have also showed, that evolution in time corresponds to the solution of a plate equation.

Spatial discretization of the split systems
We assume a two-dimensional (2-D), rectangular grid with uniform spacing ∆x for ease of explanation.Let T = {B i,j |i = 1, . . ., N ; j = 1, . . ., M } be a partition of Ω where B i,j is a square with length ∆x and N (M ) is the number of cells in the x-(y-)direction.We define φ i,j = φ(x i , y j ) where (x i , y j ) = (−∆x/2+ i∆x, −∆y/2 + j∆y).Within this section we define u = (ρ, ρv) T .
First, we treat the convective flux terms and thus recall that the eigenvalues of the Jacobians Df 1 , Df 2 are given by We discretize the fluxes f 1 and f 2 terms using a HLL-solver [40], so that the numerical fluxes are given by with λ i ± defined as for i ∈ {1, 2}.
With the numerical fluxes (3.15) we define the discrete conservative operator for the convective part f 1 (u) x + f 2 (u) y of (3.13) where ). (3.18) The values of u i±1/2,j at the cell interfaces are reconstructed component-wise using the generalized minmod limiter with θ ∈ [1, 2], The values of u i,j±1/2 at the cell interfaces are reconstructed in a similar manner.
The proposed scheme (3.17) is general and may be used in conjunction with one's favorite numerical flux replacing the HLL-flux (3.15).
The viscous part of the stress tensor can be written as Here we set ξ = 2µ + ν where µ, ν denote the coefficients of viscosity.We define the discrete central difference operators and write the discrete version of the stress tensor We combine (3.17), (3.22) together and obtain the following discretization of the system (3.13): .

Discretization of the fast system
The discretization of the fast system (3.8)consists of three parts.First, we discretize the momentum term by using central differences with the operators D x and D y defined in (3.21).Then we discretize the elliptic term as The terms ρ n ∆φ are non-conservative terms.The discretization of these terms is by no means unique and we chose for our scheme We combine (3.24), (3.25), and (3.26) together and obtain the following discretization for the fast dynamics (3.8): (3.27)

Discretization of the NSK system
For the sake of completeness we also provide the discretization of the NSK system (2.1).We use the same ideas and notations as in the previous chapters.
For the slow dynamics the discretization reads: and for the fast dynamics we have (3.29)

Boundary conditions
We have to account for the boundary conditions (2.4), (2.10) and enforce this conditions by using ghost cells.We introduce an associated cell at the edge of each element B i,j which is part of ∂Ω.For scheme (3.23), (3.27) we set and for scheme (3.28), (3.29) we set

Asymptotic preserving property
In the previous section we provided numerical schemes for the Navier-Stokes-Korteweg system and for its relaxation.This allows use to formulate the main theorem of this work.Proof: We start with the asymptotic expansion with respect to a small parameter α: For the O (1) terms we compute for the slow system ], (3.32) and the fast system and for the fast system This is the numerical scheme that we derived in (3.28), (3.29).This means that the scheme (3.23), (3.27) is an AP scheme for α → 0.

Numerical Experiments for the asymptotic-preserving scheme
In this section, we present a series of experiments that compare the performance of the AP scheme presented in the previous section with the results obtained by using the explicit discontinuous Galerkin scheme from [34].The time step ∆T was chosen according to the CFL condition (2.22) for the discontinuous Galerkin scheme with c DG and according to CFL condition (3.6) with c AP = 0.4 for the AP scheme in all numerical computations.

Computational efficiency
In this test case, we consider two different numerical schemes for the relaxation system (3.1).We set d = 1, γ = 0.16, γε 2 = 10 −5 and start with initial conditions that correspond to a two-phase density distribution.We compare the performance of the AP scheme and Discontinuous Galerkin scheme of polynomial order 1, which was introduced in [34].We use the two schemes to approximate solutions of the NSK system (2.1).
Tables 4.1 and 4.2 provide numerical evidence to the properties of the our scheme that we stated in Section 3. We take a column wise look at the CPUtimes in Table 4.1 and notice that they are constant for all values of α and for fixed ∆x.This means that for this test case the CPU time is independent of the parameter α.In contrast, the CPU-time increases as α decreases for fixed ∆x in Table 4.2.For example the AP scheme is faster by a factor 73 for α = 10000 and ∆x = 0.005 and even by a factor 84 for α = 100000 and ∆x = 0.00125.
Now we look at the discrete L 2 -distances for different α and ∆x compared to a reference solution ū, that was computed with an fully implicit DG-scheme of  The smaller bubble shrinks (t=0.02,t=0.04), and emits a shock wave as it collapses (t=0.4,t=1) .At t=4 the material seems to be in equilibrium.
polynomial order 1 for system (2.1) at ∆x = 0.000625.The last line in Table 4.2 displays the discrete L 2 -distance for a fully implicit DG-scheme of polynomial order 1 for system (2.1).Tables 4.1 and 4.2 allow to make three statements on the L 2 -errors for this test case.First, the L 2 -error for fixed α decreases as ∆x does for both schemes.Secondly, the discrete L 2 -distances decreases for fixed ∆x as α does for both schemes.Thirdly, the L 2 -error for each α and ∆x for the DG scheme is slightly better than the error for the AP scheme.However, for large α both errors are close to the error of the implicit scheme that can be regarded as a reference solution for fixed ∆x.
In summary, we observe that the AP scheme (3.23), (3.27) is able to provide solutions with the same magnitude of error as the DG scheme much faster.The AP scheme provides the numerical solution independent of the parameter α and thus is more efficient for small α.  the reference solution ū.The reference solution was computed with an implicit DG scheme.

Large density variations
We already pointed out in Section 1 that explicit numerical methods for the NSK system (2.1) are not able to deal with large density jumps.As we explained above the first-order part of the system is of mixed hyperbolic-elliptic type and one cannot use shock-capturing numerical solvers for computations.Thus, there is no chance to stabilize the numerical scheme for large density gradients.In [34], it was shown that this problem could be overcome if one uses the relaxation system (2.8).
In this 1-D test case, we introduce a modified version of the Van-der-Waals pressure equation (2.5) with the scaling parameter s = 100 to enlarge the elliptic region (cf.[5]).This enables large density jumps for phase boundaries.We apply the AP scheme (3.23)  0.16, γε 2 = 10 −5 and subject to the following initial conditions: Ω)-distance for ∆x = 0.00125 between the solution of the AP scheme and a reference solution that was computed with an implicit DG scheme.The distance decreases as α does.The second line contains the experimental order of convergence (EOC) with respect to α and the third line contains the CPU time.

Sharp interface limit
We again consider the 1-D version of the relaxation system (2.8) with α = 100000, γ = 0.16 and the following initial conditions: The initial datum is a non-equilibrium bubble that will be driven towards a two-phase equilibrium by the evolution of the relaxation system.We want to study if the AP scheme (3.23), (3.27) can deal with tiny interfaces which appear in the sharp interface limit (ε → 0) for fixed ∆x = 0.00125.the AP scheme is able to accurately capture tiny interfaces which appear in the sharp interface limit (ε → 0).This means that we are able to obtain numerical solutions that show the expected behavior even for underresolved meshes.However, we were not able to give analytical proof for this property of the scheme.
Another property of the AP scheme (3.23), (3.27) can be seen in Table 4.4.For decreasing values of ε the CPU time decreases.For γε 2 ≤ 10 5 the CPUtimes are constant.We recall (2.22) and see that the diffusive term dominates for larger ε while the convective term does so for smaller ε.Additionally, the solution of the fast system is very time consuming for the small time steps that occur for larger ε.In fact, it is a nice feature of the AP scheme, that it is faster for smaller and thus more realistic ε.

2D test case: static equilibrium
In this example, we compute a numerical solution to the relaxation system are not able to resolve the interface by the underlying grid due to CPU time constraints, the results show the desired behavior for a two-phase system.The computations for different grid sizes support the statement, that the behavior of the solution is not only qualitatively but also quantitatively correct.However, to be absolutely sure, one has to compare it with a fully resolved solution.In order to run such a simulation, one would need to implement a scheme that incorporates adaptive time and spatial stepping, which is out of the scope of  the current paper.

Conclusion/Outlook
In this work, we presented a new asymptotic-preserving scheme for a relaxation of the Navier-Stokes-Korteweg system.The method has two key ingredients.First, we introduced a modified pressure that guarantees that the Euler part of the system is hyperbolic.Second, we split the relaxation system into a non linear conservative system for the slow dynamics and a non conservative stiff system for the fast dynamics that come from the Korteweg term.
We showed that our scheme is asymptotic preserving, i.e. it converges to a numerical scheme for the Navier-Stokes-Korteweg equations on the fully discrete level.We supported this result with several numerical tests.In all test cases, the solutions showed the expected behavior.We were able to compute large density ratios and small interfacial widths for fixed grid widths and we were  able to obtain approximate solutions of the same quality as those of an explicit scheme, but with significantly smaller computational costs.Additionally, our scheme showed to be asymptotic preserving in the sharp interface limit in our numerical tests.Note that up to now we have not been able to give analytical proof to that.
The numerical proof of the asymptotic preserving property is one goal for future work.This is by no means an easy task, because the system has a jump of the density at the phase boundaries.Therefore one needs to impose special coupling conditions at the interface and to introduce an interface tracking algorithm.Another main goal is the extension of the numerical algorithm to three space dimensions.In order to accomplish this, we will look to introduce an adaptive time and spatial stepping and to use a discontinuous Galerkin formulation instead of the finite volume formulation.These two extensions guarantee that the solution of the linear system remains efficient, because we are able to control the number of cells, even if we need to refine the grid at the interface.
This could prove to be important for the simulation of realistic two-phase flow examples such as underwater explosions.

. 5 )
Thereby, b, b 1 , R are positive constants and T * is the fixed temperature.If T * is chosen small enough, the pressure p is monotone decreasing in some non-empty density interval.This structure allows one to define phases.If the density ρ lies in the interval (0, α 1 ], ((α 1 , α 2 )), {[α 2 , b)} the corresponding fluid state is called vapor (spinodal) {liquid}, see Figure 2.1 for an illustration.

Figure 4 .
1 shows the evolution of the density at different times that was computed with the AP scheme (3.23), (3.27).

Figure 4 . 1 :
Figure 4.1: Density evolution for the relaxation system (3.1) and α = 100 with initial datum (2.21).The initial configuration at t=0 corresponds to two bubbles in a box filled with liquid.

L 2 -
distances between the solution of the discontinuous Galerkin scheme and the reference solution ū.The last line contains the discrete L 2 -distance for an explicit DG discretization of the NSK system.The reference solution was computed with an implicit DG scheme.

Figure 4 . 3 :
Figure 4.3: Approximate density distribution at t = 15 for different ε.The right Figure displays a zoom of the upper corner of the right interface

( 2 . 8 )
Figures 4.4 and 4.5.The density varies between 0.3 (blue) and 1.8 (red).At t = 0.06, t = 0.1 bubbles collapse and emit shock waves.At t = 1, t = 3 the smaller bubbles shrink and the biggest bubble grows.At t = 5 the fluid seems to be in an equilibrium configuration.Note that the position of the bubbles and the coalescence dynamics are the same for both gird widths.Even though we

Table 4 .
1: CPU-times and discrete L 2 -distances between the solution of the AP scheme and

Table 4 .
3 and Figure4.2 indicate that it is possible to simulate large density rations with the AP scheme (3.23),(3.27).The first and second line ofTable 4.3show that the scheme is asymptotic preserving in this test case, as the discrete L 2 -distance decreases as α does.The third line underlines, that, analogously to the previous test case, the CPU time does not depend on the Korteweg

Table 4 .
4: CPU time for different ε.The CPU-time decreases as ε does.