Monotonic Upwind Scheme For Conservation Laws Research Articles

MATRICS: The implicit matrix-free Eulerian hydrodynamics solver

Context. There exists a zoo of different astrophysical fluid dynamics solvers, most of which are based on an explicit formulation and hence stability-limited to small time steps dictated by the Courant number expressing the local speed of sound. With this limitation, the modeling of low-Mach-number flows requires small time steps that introduce significant numerical diffusion, and a large amount of computational resources are needed. On the other hand, implicit methods are often developed to exclusively model the fully incompressible or 1D case since they require the construction and solution of one or more large (non)linear systems per time step, for which direct matrix inversion procedures become unacceptably slow in two or more dimensions. Aims. In this work, we present a globally implicit 3D axisymmetric Eulerian solver for the compressible Navier–Stokes equations including the energy equation using conservative formulation and a fully simultaneous approach. We use the second-order-in-time backward differentiation formula for temporal discretization as well as the κ scheme for spatial discretization. We implement different limiter functions to prohibit the occurrence of spurious oscillations in the vicinity of discontinuities. Our method resembles the well-known monotone upwind scheme for conservation laws (MUSCL). We briefly present efficient solution methods for the arising sparse and nonlinear system of equations. Methods. To deal with the nonlinearity of the Navier–Stokes equations we used a Newton iteration procedure in which the required Jacobian matrix-vector product was reconstructed with a first-order finite difference approximation to machine precision in a matrix-free way. The resulting linear system was solved either completely matrix-free with a combination of a sufficient Krylov solver and an approximate Jacobian preconditioner or semi-matrix-free with the incomplete lower upper factorization technique as a preconditioner. The latter was dependent on a standalone approximation of the Jacobian matrix, which was optionally calculated and needed solely for the purpose of preconditioning. Results. We show our method to be capable of damping sound waves and resolving shocks even at Courant numbers larger than one. Furthermore, we prove the method’s ability to solve boundary value problems like the cylindrical Taylor-Couette flow (TC), including viscosity, and to model transition flows. To show the latter, we recover predicted growth rates for the vertical shear instability, while choosing a time step orders of magnitude larger than the explicit one. Finally, we verify that our method is second order in space by simulating a simplistic, stationary solar wind.

Flood Prediction with Two-Dimensional Shallow Water Equations: A Case Study of Tongo-Bassa Watershed in Cameroon

As a result of urbanization, combined with the anthropogenic effects of climate change, natural events such as floods are showing increasingly adverse impacts on human existence. This study proposes a new model, based on shallow water equations, that is able to predict these floods and minimize their impacts. The first-order finite volume method (FVM), the Harten Lax and van Leer (HLL) scheme, and the monotone upwind scheme for conservation laws (MUSCL) are applied in the model. In addition, a virtual boundary cell approach is adopted to achieve a monotonic solution for both interior and boundary cells and flux computations at the boundary cells. The model integrates the infiltration parameters recorded in the area, as well as the Manning coefficient specific to each land-cover type of the catchment region. The results provided were mapped to highlight the potential flood zones and the distribution of water heights throughout the catchment region at any given time, as well as that at the outlet. It has been observed that when standard infiltration and the Manning parameters were selected, the floodable surface increased, as expected, with the increasing rainfall intensity and duration of the simulation. With sufficient infiltration, only a portion of the water tends to stagnate and flow off on the surface toward the outlet. A sensitivity analysis of certain parameters, such as rainfall data and the final infiltration coefficient in the lower watershed of the littoral region, was conducted; the results show that the model simulates well the general character of water flow in the watershed. Finally, the model’s validation using field-collected parameters during the flood of 25 July 2017 and 18 to 22 July 2016 in the Grand Ouaga basin in Burkina reveals Nash–Sutcliffe values of 0.7 and 0.73, respectively.

Assessment of VoF based numerical scheme for bubble rise in isothermal liquid layer, and some new insight in thermally stratified liquid layers

This paper aims at (a) comparative assessment of the different volume of fluid (VoF) based numerical schemes for a rising, single air bubble in isothermal liquid (water) layers and (b) investigation of bubble rise in thermally stratified liquid (therminol) layers. For numerical investigation of bubble rise in isothermal liquid layers, three available numerical schemes, Upwind, Quadratic Upstream Interpolation for Convective Kinematics (QUICK) and Monotonic Upwind Scheme for Conservation Laws (MUSCL) schemes, combined with the pressure-velocity coupling(p-v) approaches, such as, Pressure Implicit with Splitting of Operators (PISO) and Semi Implicit Method for Pressure-Linked Equations (SIMPLE), are considered. The wake analysis revealed that the region of influence grows beyond 10 times the diameter of the bubble. Moreover, based on the comparative assessment, MUSCL scheme with PISO is selected, for investigating the bubble rise in thermally stratified therminol layers. Based on these investigations, a 3D diagram, describing bubble shape as f (Ra, Eo, Ga), is proposed, and a new insight to the micro-convection, inside the rising bubbles is provided. Furthermore, a time-scale analysis is performed to describe the heat transfer mechanisms, (a) inside air bubble, and (b) between air bubble and the external surrounding liquid. Thus, the findings will be useful for the design of heat exchangers or cooling devices, which rely on the heat transfer augmentation with rising air bubble.

Catchment-scale multi-process modeling with local time stepping

Numerical modeling is an important tool for prediction, analysis and understanding of the dynamics of land surface processes. To increase the usage and impact of such tools, it is crucial to decrease runtime by increasing computational efficiency. Dynamic processes such as water flow are typically described by higher-order differential equations. Solving these accurately requires numerical integration over time, where numerical errors depend on the time steps taken. Typically, flow simulation use the smallest required time steps in a model’s domain to simulate flow. In this paper, we analyze the usage of local time stepping, for catchment-scale simulation of land surface processes such as water flow, infiltration, slope stability and landslide runout. In such a scheme, temporal integration is cell specific, allowing for higher numerical efficiency. The implemented scheme works with fully free local time steps that are synchronized only for visualization. We implement this method in a monotonic upwind scheme for conservation laws (MUSCL). We investigate the influence on stability and the resulting changes in computation time and accuracy in a hydrology-coupled, catchment-scale flood simulation. Results show that local time stepping can be implemented in a total variation diminishing (TVD) numerical scheme that is second-order spatially accurate. Simulation results in both 1D dam-break scenarios and catchment-scale flash flood scenarios show insignificant changes in modeling result, while computation time reduces with over 50%. Finally, the method is successfully implemented in a multi-process lands surface model with hydrology, flooding, slope failure, and runout. The implementation of a local time stepping for computation of dynamic land surface processes could be implemented widely for increased computational efficiency without significant loss of accuracy.

Extension of artificial viscosity technique for solving 2D non-hydrostatic shallow water equations

Abstract This study highlights the solution of the 2D non-hydrostatic shallow water equations using an artificial viscosity (AV) technique. The model recently proposed by the first author in Ginting (2017) is extended in this paper by reconstructing the flow variables using the Monotonic Upwind Scheme for Conservation Laws (MUSCL) technique to achieve second-order accuracy and by adding the non-hydrostatic terms, where the Runge–Kutta second-order scheme is used for temporal discretization. To reduce computational time, a novel way is proposed, where the AV technique as well as the non-hydrostatic terms are computed in a hybrid manner, only once per single time step. The solution of the non-hydrostatic terms forms a system of linear equations and the conjugate gradient method is therefore employed to solve it. The results show this technique is robust and accurate for simulating a set of problems of coastal engineering, thus could become a promising method for non-hydrostatic shallow flow applications.

AUSM‐like expression of HLLC and its all‐speed extension

SummaryAdvection upstream splitting method (AUSM) and Harten‐Lax‐van Leer with contact (HLLC) are two popular families of flux functions. The AUSM is simple and requires no eigenstructure, which facilitates its extensions to general equations of state. Furthermore, one of its variants, simple low‐dissipation AUSM (SLAU), is applicable to all speeds and features removal of parameter setting by the user. HLLC, on the other hand, clearly defines three distinct waves in Riemann problem, namely, left‐running and right‐running acoustic waves, and entropy wave. This paper demonstrates that HLLC can be written in a very similar form with the AUSM family and that the similar manner in extending AUSM family to all speeds is easily incorporated into HLLC in this AUSM‐like form. Then, we combine the strengths of the both flux functions and offer a new inviscid numerical flux function within the framework of monotone upwind scheme for conservation laws (MUSCL) in computational fluid dynamics (CFD) for Euler and Navier‐Stokes equations. The resultant HLLC with low dissipation (HLLCL) numerical flux can compute low Mach number flows and sound propagations at the same time with high accuracy, as demonstrated by one‐dimensional and two‐dimensional numerical examples. Furthermore, the results indicate that its extensions to general fluids such as supercritical fluids are encouraging.

A Hydraulic Model for Multiphase Flow Based on the Drift Flux Model in Managed Pressure Drilling

Managed pressure drilling (MPD) is a drilling technique used to address the narrow density window under complex geological environments. It has widespread applications in the exploration and exploitation of oil and gas, both onshore and offshore. In this study, to achieve effective control of the downhole pressure to ensure safety, a gas–liquid two-phase flow model based on the drift flux model is developed to describe the characteristics of transient multiphase flow in the wellbore. The advection upwind splitting method (AUSM) numerical scheme is used to assist with calculation and analysis, and the monotonic upwind scheme for conservation laws (MUSCLs) technique with second-order precision is adopted in combination with the Van Leer slope limiter to improve precision. Relevant data sourced from prior literature are used to validate the suggested model, the results of which reveal an excellent statistical consistency. Further, the influences of various parameters in a field application, including backpressure, density, and mass flow, are analyzed. Over the course of later-stage drilling, a combination of wellhead backpressure and displacement is recommended to exercise control.

Magnetohydrodynamic simulation code CANS+: Assessments and applications

Abstract We present a new magnetohydrodynamic (MHD) simulation package with the aim of providing accurate numerical solutions to astrophysical phenomena where discontinuities, shock waves, and turbulence are inherently important. The code implements the Harten–Lax–van Leer–discontinuitues (HLLD) approximate Riemann solver, the fifth-order-monotonicity-preserving interpolation (MP5) scheme, and the hyperbolic divergence cleaning method for a magnetic field. This choice of schemes has significantly improved numerical accuracy and stability, and saved computational costs in multidimensional problems. Numerical tests of one- and two-dimensional problems show the advantages of using the high-order scheme by comparing with results from a standard second-order total variation diminishing monotonic upwind scheme for conservation laws (MUSCL) scheme. The present code enables us to explore the long-term evolution of a three-dimensional accretion disk around a black hole, in which compressible MHD turbulence causes continuous mass accretion via nonlinear growth of the magneto-rotational instability (MRI). Numerical tests with various computational cell sizes exhibits a convergent picture of the early nonlinear growth of the MRI in a global model, and indicates that the MP5 scheme has more than twice the resolution of the MUSCL scheme in practical applications.

Comparison of Shallow Water Solvers: Applications for Dam-Break and Tsunami Cases with Reordering Strategy for Efficient Vectorization on Modern Hardware

We investigate in this paper the behaviors of the Riemann solvers (Roe and Harten-Lax-van Leer-Contact (HLLC) schemes) and the Riemann-solver-free method (central-upwind scheme) regarding their accuracy and efficiency for solving the 2D shallow water equations. Our model was devised to be spatially second-order accurate with the Monotonic Upwind Scheme for Conservation Laws (MUSCL) reconstruction for a cell-centered finite volume scheme—and be temporally fourth-order accurate using the Runge–Kutta fourth-order method. Four benchmark cases of dam-break and tsunami events dealing with highly-discontinuous flows and wet–dry problems were simulated. To this end, we applied a reordering strategy for the data structures in our code supporting efficient vectorization and memory access alignment for boosting the performance. Two main features are pointed out here. Firstly, the reordering strategy employed has enabled highly-efficient vectorization for the three solvers investigated on three modern hardware (AVX, AVX2, and AVX-512), where speed-ups of 4.5–6.5× were obtained on the AVX/AVX2 machines for eight data per vector while on the AVX-512 machine we achieved a speed-up of up to 16.7× for 16 data per vector, all with singe-core computation; with parallel simulations, speed-ups of up to 75.7–121.8× and 928.9× were obtained on AVX/AVX2 and AVX-512 machines, respectively. Secondly, we observed that the central-upwind scheme was able to outperform the HLLC and Roe schemes 1.4× and 1.25×, respectively, by exhibiting similar accuracies. This study would be useful for modelers who are interested in developing shallow water codes.

Investigation of Numerical Dissipation in Classical and Implicit Large Eddy Simulations

The quantitative measure of dissipative properties of different numerical schemes is crucial to computational methods in the field of aerospace applications. Therefore, the objective of the present study is to examine the resolving power of Monotonic Upwind Scheme for Conservation Laws (MUSCL) scheme with three different slope limiters: one second-order and two third-order used within the framework of Implicit Large Eddy Simulations (ILES). The performance of the dynamic Smagorinsky subgrid-scale model used in the classical Large Eddy Simulation (LES) approach is examined. The assessment of these schemes is of significant importance to understand the numerical dissipation that could affect the accuracy of the numerical solution. A modified equation analysis has been employed to the convective term of the fully-compressible Navier–Stokes equations to formulate an analytical expression of truncation error for the second-order upwind scheme. The contribution of second-order partial derivatives in the expression of truncation error showed that the effect of this numerical error could not be neglected compared to the total kinetic energy dissipation rate. Transitions from laminar to turbulent flow are visualized considering the inviscid Taylor–Green Vortex (TGV) test-case. The evolution in time of volumetrically-averaged kinetic energy and kinetic energy dissipation rate have been monitored for all numerical schemes and all grid levels. The dissipation mechanism has been compared to Direct Numerical Simulation (DNS) data found in the literature at different Reynolds numbers. We found that the resolving power and the symmetry breaking property are enhanced with finer grid resolutions. The production of vorticity has been observed in terms of enstrophy and effective viscosity. The instantaneous kinetic energy spectrum has been computed using a three-dimensional Fast Fourier Transform (FFT). All combinations of numerical methods produce a k − 4 spectrum at t * = 4 , and near the dissipation peak, all methods were capable of predicting the k − 5 / 3 slope accurately when refining the mesh.

Quantification of numerical diffusivity due to TVD schemes in the advection equation

In this study, the numerical diffusivity νnum inherent to the Roe–MUSCL scheme has been quantified for the scalar advection equation. The Roe–MUSCL scheme employed is a combination of: (1) the standard extension of the original Roeʼs formulation to the advection equation, and (2) van Leerʼs Monotone Upwind Scheme for Conservation Laws (MUSCL) technique that applies a linear variable reconstruction in a cell along with a scaled limiter function. An explicit expression is derived for the numerical diffusivity in terms of the limiter function, the distance between the cell centers on either side of a face, and the face-normal velocity. The numerical diffusivity formulation shows that a scaled limiter function is more appropriate for MUSCL in order to consistently recover the central-differenced flux at the maximum value of the limiter. The significance of the scaling factor is revealed when the Roe–MUSCL scheme, originally developed for 1-D scenarios, is applied to 2-D scalar advection problems. It is seen that without the scaling factor, the MUSCL scheme may not necessarily be monotonic in multi-dimensional scenarios. Numerical diffusivities of the minmod, superbee, van Leer and Barth–Jesperson TVD limiters were quantified for four problems: 1-D advection of a step function profile, and 2-D advection of step, sinusoidal, and double-step profiles. For all the cases, it is shown that the superbee scheme provides the lowest numerical diffusivity that is also most confined to the vicinity of the discontinuity. The minmod scheme is the most diffusive, as well as active in regions away from high gradients. As expected, the grid resolution study demonstrates that the magnitude and the spatial extent of the numerical diffusivity decrease with increasing resolution.

Accuracy of high‐order density‐based compressible methods in low Mach vortical flows

SUMMARYAt low Mach numbers, Godunov‐type approaches, based on the method of lines, suffer from an accuracy problem. This paper shows the importance of using the low Mach number correction in Godunov‐type methods for simulations involving low Mach numbers by utilising a new, well‐posed, two‐dimensional, two‐mode Kelvin–Helmholtz test case. Four independent codes have been used, enabling the examination of several numerical schemes. The second‐order and fifth‐order accurate Godunov‐type methods show that the vortex‐pairing process can be captured on a low resolution with the low Mach number correction applied down to 0.002. The results are compared without the low Mach number correction and also three other methods, a Lagrange‐remap method, a fifth‐order accurate in space and time finite difference type method based on the wave propagation algorithm, and fifth‐order spatial and third‐order temporal accurate finite volume Monotone Upwind Scheme for Conservation Laws (MUSCL) approach based on the Godunov method and Simple Low Dissipation Advection Upstream Splitting Method (SLAU) numerical flux with low Mach capture property. The ability of the compressible flow solver of the commercial software, ANSYS FLUENT, in solving low Mach flows is also demonstrated for the two time‐stepping methods provided in the compressible flow solver, implicit and explicit. Results demonstrate clearly that a low Mach correction is required for all algorithms except the Lagrange‐remap approach, where dissipation is independent of Mach number. © 2013 Crown copyright. International Journal for Numerical Methods in Fluids. © 2013 John Wiley & Sons, Ltd.

New approaches for computation of low Mach number flows

By using all speed numerical flux schemes, such as SLAU [Simple Low Dissipation AUSM (Advection Upstream Splitting Method)], in MUSCL (Monotone Upwind Scheme for Conservation Laws) approach for compressible CFD, low Mach number flows can be computed without loss of accuracy nor parameter tuning. For an efficient computation, this paper deals with new approaches of implicit time integration method. In this approach, the large sparse matrix system, which consists of flux Jacobian of numerical flux function, has to be solved in each time step. Firstly, a simple Gauss–Seidel iteration method named TC-PGS1(Time Consistent Preconditioned Gauss–Seidel 1) which has flavor of the time derivative preconditioning is introduced. Secondary, we tried to use FGMRES (k) (Flexible Generalized Minimum Residual Method) to solve the non-diagonal dominant linear system arising from Jacobian of flux function SLAU. TC-PGS1 is also used as the matrix preconditioner for FGMRES (k). Optimal parameters for FGMRES (k) is investigated numerically and the performances on computational efficiency of the new methods are compared. It is indicated that FGMRES (k) has apparent advantage on computation of low Mach number flows with sound propagation, however, simpler TC-PGS1 has comparable performance if only flow fields are of interest.

Modeling of Shock Wave Propagation in Large Amplitude Ultrasound

The Rankine-Hugoniot relation for shock wave propagation describes the shock speed of a nonlinear wave. This paper investigates time-domain numerical methods that solve the nonlinear parabolic wave equation, or the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation, and the conditions they require to satisfy the Rankine-Hugoniot relation. Two numerical methods commonly used in hyperbolic conservation laws are adapted to solve the KZK equation: Godunov's method and the monotonic upwind scheme for conservation laws (MUSCL). It is shown that they satisfy the Rankine-Hugoniot relation regardless of attenuation. These two methods are compared with the current implicit solution based method. When the attenuation is small, such as in water, the current method requires a degree of grid refinement that is computationally impractical. All three numerical methods are compared in simulations for lithotripters and high intensity focused ultrasound (HIFU) where the attenuation is small compared to the nonlinearity because much of the propagation occurs in water. The simulations are performed on grid sizes that are consistent with present-day computational resources but are not sufficiently refined for the current method to satisfy the Rankine-Hugoniot condition. It is shown that satisfying the Rankine-Hugoniot conditions has a significant impact on metrics relevant to lithotripsy (such as peak pressures) and HIFU (intensity). Because the Godunov and MUSCL schemes satisfy the Rankine-Hugoniot conditions on coarse grids, they are particularly advantageous for three-dimensional simulations.

A third-order explicit central scheme for open channel flow simulations

A third-order explicit time-stepping model for the simulations in open channels flow is presented. The time discretization of one-dimensional Saint- Venant's basic equations is obtained by Taylor series expansion, formulated in multi-stage approach, whereas the spatial discretization is obtained by finite difference central scheme. First- and second-order models have been obtained by using piecewise constant and piecewise linear MUSCL (monotone upwind scheme for conservation laws) approximations. Afterwards, the third-order model proposed is obtained by a piecewise quadratic polynomial approximation. Using Fourier's linear analysis, the stability and the accuracy of the scheme are investigated. The numerical results, for predicting dam-break in a horizontal and frictionless channel and for the representation of the hydraulic jump in prismatic and non-prismatic channels, with non-uniform bottom slope, are evaluated, and compared with the corresponding analytical and measured solutions.

Discretization of Integral Equations Describing Flow in Nonprismatic Channels with Uneven Beds

Application of the finite-volume method in one dimension for open channel flow predictions mandates the direct discretization of integral equations for mass conservation and momentum balance. The integral equations include source terms that account for the forces due to changes in bed elevation and channel width, and an exact expression for these source term integrals is presented for the case of a trapezoidal channel cross section whereby the bed elevation, bottom width, and inverse side slope are defined at cell faces and assumed to vary linearly and uniformly within each cell, consistent with a second-order accurate solution. The expressions may be used in the context of any second-order accurate finite-volume scheme with channel properties defined at cell faces, and it is used here in the context of the Monotone Upwind Scheme for Conservation Laws (MUSCL)-Hancock scheme which has been adopted by many researchers. Using these source term expressions, the MUSCL-Hancock scheme is shown to preserve stationarity, accurately converge to the steady state in a frictionless flow test problem, and perform well in field applications without the need for upwinding procedures previously reported in the literature. For most applications, an approximate, point-wise treatment of the bed slope and nonprismatic source terms can be used instead of the exact expression and, in contrast to reports on other finite-volume-based schemes, will not cause unphysical oscillations in the solution.

High-resolution and non-oscillatory solution of the St. Venant equations in non-rectangular and non-prismatic channels

A scheme to model open channel flow over wet and dry beds in non-rectangular and non-prismatic channels is presented. The scheme is second-order accurate, stable for Courant numbers up to unity, and monotonicity preserving. The scheme solves the St. Venant equations using a Godunov-type finite volume method. Mass and momentum fluxes are computed using a Roe-type Riemann solver, the MUSCL (Monotone Upwind Scheme for Conservation Laws) approach is applied for second-order spatial accuracy, and a treatment is introduced to model the hydrostatic pressure force exerted by the channel walls in the stream wise direction. The treatment permits momentum fluxes and the channel wall force to be balanced to numerical precision, preventing the artificial acceleration of the flow. Comparisons between model results, exact solutions, and experimental data show that the scheme is robust. Accurate and monotone results are obtained in the presence of discontinuities, supercritical flow, subcritical flow, transcritical flow, and dry-bed flow problems without the need for special front-tracking approaches or deforming grids. In addition, the scheme will conserve mass to numerical precision in all applications.