Riemann solvers and boundary conditions for two‐dimensional shallow water simulations

Most existing algorithms for two‐dimensional shallow water simulations treat multi‐dimensional waves using wave splitting or time splitting. This often results in anisotropy of the computed flow. Both wave splitting and time splitting are based on a local decomposition of the multi‐dimensional problem into one‐dimensional, orthogonal problems. Therefore, these algorithms handle boundary conditions in a very similar way to classical one‐dimensional algorithms. This should be expected to trigger a dependence of the number of boundary conditions on the direction of the flow at the boundaries. However, most computational codes based on alternate directions do not exhibit such sensitivity, which seems to contradict the theory of existence and uniqueness of the solution. The present paper addresses these issues. A Riemann solver is presented that aims to convert two‐dimensional Riemann problems into a one‐dimensional equivalent Riemann problem (ERP) at the interfaces between the computational cells. The ERP is derived by applying the theory of bicharacteristics at each end of the interface and by performing a linear averaging along the interface. The proposed approach is tested against the traditional one‐dimensional approach on the classical circular dambreak problem. The results show that the proposed solver allows the isotropy of the solution to be better preserved. Use of the two‐dimensional solver with a first‐order scheme may give better results than use of a second‐order scheme with a one‐dimensional solver. The theory of bicharacteristics is also used to discuss the issue of boundary conditions. It is shown that, when the flow is subcritical, the number of boundary conditions affects the accuracy of the solution, but not its existence and uniqueness. When only one boundary condition is to be prescribed, it should not be the velocity in the direction parallel to the boundary. When two boundary conditions are to be prescribed, at least one of them should involve the component of the velocity in the direction parallel to the boundary. Copyright © 2003 John Wiley & Sons, Ltd.     

High resolution Godunov‐type schemes with small stencils

Higher‐order Godunov‐type schemes have to cope with the following two problems: (i) the increase in the size of the stencil that make the scheme computationally expensive, and (ii) the monotony‐preserving treatments (limiters) that must be implemented to avoid oscillations, leading to strong damping of the solution, in particular linear waves (e.g. acoustic waves). When too compressive, limiting procedures may also trigger the instability of oscillatory numerical solutions (e.g. in advection–dispersion phenomena) via the artificial amplification of the shorter modes. The present paper proposes a new approach to carry out the reconstruction. In this approach, the values of the flow variable at the edges of the computational cells are obtained directly from the reconstruction within these cells. This method is applied to the MUSCL and DPM schemes for the solution of the linear advection equation. The modified DPM scheme can capture contact discontinuities within one computational cell, even after millions of time steps at Courant numbers ranging from 1 to values as low as 10^‐4. Linear waves are subject to negligible damping. Application of the method to the DPM for one‐dimensional advection–dispersion problems shows that the numerical instability of oscillatory solutions caused by the over compressive, original DPM limiter is eliminated. One‐ and two‐dimensional shallow water simulations show an improvement over classical methods, in particular for two‐dimensional problems with strongly distorted meshes. The quality of the computational solution in the two‐dimensional case remains acceptable even for mesh aspect ratios Δx/Δy as large as 10. The method can be extend to the discretization of higher‐order PDEs, allowing third‐order space derivatives to be discretized using only two cells in space. Copyright © 2004 John Wiley & Sons, Ltd.     

固定网格下的特征线法求解溃坝问题

溃坝问题是典型的非线性双曲方程的Riemann问题,其数值求解的难点在于对间断面的捕捉以及避免间断面处在数值计算过程中产生数值色散,因而为求解此问题所产生的各种数值计算方法的优劣也体现在这两个方面。本文针对溃坝问题提出一种新的计算方法。该方法基于对偶变量推导的浅水波方程,根据方程的特点,从方程的特征值和黎曼不变量出发,采用高精度的激波捕捉方法计算黎曼不变量的位置随时间的变化,然后映射至不随时间变化的固定网格。根据黎曼不变量的位置,采用保形分段三次Hermite插值将物理量映射至网格节点。计算结果显示,该方法不仅操作简单,计算量小,而且结果准确。     

Using fully implicit conservative statements to close open boundaries passing through recirculations

The numerical solution of the fluid flow governing equations requires the implementation of certain boundary conditions at suitable places to make the problem well‐posed. Most of numerical strategies exhibit weak performance and obtain inaccurate solutions if the solution domain boundaries are not placed at adequate locations. Unfortunately, many practical fluid flow problems pose difficulty at their boundaries because the required information for solving the PDE's is not available there. On the other hand, large solution domains with known boundary conditions normally need a higher number of mesh nodes, which can increase the computational cost. Such difficulties have motivated the CFD workers to confine the solution domain and solve it using artificial boundaries with unknown flow conditions prevailing there. In this work, we develop a general strategy, which enables the control‐volume‐based methods to close the outflow boundary at arbitrary locations where the flow conditions are not known prior to the solution. In this regard, we extend suitable conservative statements at the outflow boundary. The derived statements gradually detect the correct boundary conditions at arbitrary boundaries via an implicit procedure using a finite element volume method. The extended statements are validated by solving the truncated benchmark backward‐facing step problem. The investigation shows that the downstream boundary can pass through a recirculation zone without deteriorating the accuracy of the solution either in the domain or at its boundaries. The results indicate that the extended formulation is robust enough to be employed in solution domains with unknown boundary conditions. Copyright © 2006 John Wiley & Sons, Ltd.     

Solution property preserving reconstruction for finite volume scheme: a boundary variation diminishing+multidimensional optimal order detection framework

The purpose of this work is to build a general framework to reconstruct the underlying fields within a finite volume (FV) scheme solving a hyperbolic system of PDEs (Partial Differential Equations). In an FV context, the data are piecewise constants per computational cell and the physical fields are reconstructed taking into account neighbor cell values. These reconstructions are further used to evaluate the physical states usually used to feed a Riemann solver which computes the associated numerical fluxes. The physical field reconstructions must obey some properties linked to the system of PDEs such as the positivity, but also some numerically based ones, like an essentially nonoscillatory behavior. Moreover, the reconstructions should be highly accurate for smooth flows and robust/stable for discontinuous solutions. To ensure a solution property preserving reconstruction, we introduce a methodology to blend high-/low-order polynomials and hyperbolic tangent reconstructions. A boundary variation diminishing algorithm is employed to select the best reconstruction in each cell. An a posteriori MOOD detection procedure is employed to ensure the positivity by recomputing the rare problematic cells using a robust first-order FV scheme. We illustrate the performance of the proposed scheme via the numerical simulations for some benchmark tests which involve vacuum or near vacuum states, strong discontinuities, and also smooth flows. The proposed scheme maintains high accuracy on smooth profile, preserves the positivity and can eliminate the oscillations in the vicinity of discontinuities while maintaining sharper discontinuity with superior solution quality compared to classical highly accurate FV schemes.     

An efficient discontinuous Galerkin method for aeroacoustic propagation

An efficient discontinuous Galerkin formulation is applied to the solution of the linearized Euler equations and the acoustic perturbation equations for the simulation of aeroacoustic propagation in two‐dimensional and axisymmetric problems, with triangular and quadrilateral elements. To improve computational efficiency, a new strategy of variable interpolation order is proposed in addition to a quadrature‐free approach and parallel implementation. Moreover, an accurate wall boundary condition is formulated on the basis of the solution of the Riemann problem for a reflective wall. Time discretization is based on a low dissipation formulation of a fourth‐order, low storage Runge–Kutta scheme. Along the far‐field boundaries a perfectly matched layer boundary condition is used. For the far‐field computations, the integral formulation of Ffowcs Williams and Hawkings is coupled with the near‐field solver. The efficiency and accuracy of the proposed variable order formulation is assessed for realistic geometries, namely sound propagation around a high‐lift airfoil and the Munt problem. Copyright © 2011 John Wiley & Sons, Ltd.     

Sensitivity of the 1D shallow water equations with source terms: Solution method for discontinuous flows

A finite volume‐based numerical technique is presented concerning the sensitivity of the solution of the one‐dimensional Shallow Water Equations with scalar transport. An approximate Riemann solver is proposed for direct sensitivity calculation even in the presence of discontinuous solutions. The Shallow Water Sensitivity Equations are first derived as well as the expressions of the sensitivity source terms, initial and boundary conditions. The numerical technique is then detailed and application examples are provided to assess the method's efficiency in estimating the sensitivity to different parameters (friction coefficient and initial and boundary conditions). The application of the dam‐break problem to a trapezoidal channel is also provided. The comparison with the analytical solution and the classical empirical approach illustrates the usefulness of the direct sensitivity calculation. Copyright © 2010 John Wiley & Sons, Ltd.     

High order schemes for the scalar transport equation

A finite volume hybrid scheme for the spatial discretization that combines a fixed stencil and a stencil determined by the classical essentially non‐oscillatory (ENO) scheme is presented. Evolution equations are obtained for the mean values of each cell by means of piecewise interpolation. Time discretization is accomplished by a classical fourth‐order Runge–Kutta. Interpolation polynomials are determined using information of adjacent cells. While smooth regions are interpolated by means of a fixed molecule, discontinuous or sharp regions are interpolated by the classical ENO algorithm. The algorithm estimates the interpolation error at each time step by means of two interpolants of order q and q+1. The main computational load of the resultant scheme is in the interpolation, which is performed by the divided differences table. This table involves O(qN) operations, where q is the interpolation order and N is the number of cells. Finally, linear test cases of continuous and discontinuous initial conditions are integrated to see the goodness of the hybrid scheme. It is well known that, for some particular initial conditions, the classical ENO scheme does not perform properly, not attaining the truncation error of the scheme. It is shown that, for the smooth initial condition, sin⁴(x), the classical ENO scheme does not preserve the character of stability of the initial value problem, giving rise to unstable eigenvalues. The proposed hybrid scheme solves this problem, choosing a fixed stencil over the whole computational domain. The resultant schemes are equivalent to the classical finite difference schemes, which preserve the character of stability. It is also known that the same degeneracy of the error can be encountered for discontinuous solutions. It is shown for the initial discontinuous solution, e^−x, that the classical ENO algorithm does not perform properly due to the conflict between the selection of the stencil to smoother regions (downwind region) and the hyperbolic character of the problem, which obliges us to take information from downwind. The proposed hybrid scheme solves this problem by choosing a fixed stencil over the whole computational domain except at the discontinuity. Copyright © 2001 John Wiley & Sons, Ltd.     

A general approximate‐state Riemann solver for hyperbolic systems of conservation laws with source terms

An approximate‐state Riemann solver for the solution of hyperbolic systems of conservation laws with source terms is proposed. The formulation is developed under the assumption that the solution is made of rarefaction waves. The solution is determined using the Riemann invariants expressed as functions of the components of the flux vector. This allows the flux vector to be computed directly at the interfaces between the computational cells. The contribution of the source term is taken into account in the governing equations for the Riemann invariants. An application to the water hammer equations and the shallow water equations shows that an appropriate expression of the pressure force at the interface allows the balance with the source terms to be preserved, thus ensuring consistency with the equations to be solved as well as a correct computation of steady‐state flow configurations. Owing to the particular structure of the variable and flux vectors, the expressions of the fluxes are shown to coincide partly with those given by the HLL/HLLC solver. Computational examples show that the approximate‐state solver yields more accurate solutions than the HLL solver in the presence of discontinuous solutions and arbitrary geometries. Copyright © 2006 John Wiley & Sons, Ltd.     

A discontinuous Galerkin‐based sharp‐interface method to simulate three‐dimensional compressible two‐phase flow

A numerical method for the simulation of compressible two‐phase flows is presented in this paper. The sharp‐interface approach consists of several components: a discontinuous Galerkin solver for compressible fluid flow, a level‐set tracking algorithm to follow the movement of the interface and a coupling of both by a ghost‐fluid approach with use of a local Riemann solver at the interface. There are several novel techniques used: the discontinuous Galerkin scheme allows locally a subcell resolution to enhance the interface resolution and an interior finite volume Total Variation Diminishing (TVD) approximation at the interface. The level‐set equation is solved by the same discontinuous Galerkin scheme. To obtain a very good approximation of the interface curvature, the accuracy of the level‐set field is improved and smoothed by an additional P_NP_M‐reconstruction. The capabilities of the method for the simulation of compressible two‐phase flow are demonstrated for a droplet at equilibrium, an oscillating ellipsoidal droplet, and a shock‐droplet interaction problem at Mach 3. Copyright © 2015 John Wiley & Sons, Ltd.     

Spectrally accurate method for analysis of stationary flows of second‐order fluids in rough micro‐channels

A spectral method for the analysis of stationary flows of second‐order fluids in rough micro‐channels is developed. The algorithm employs a fixed computational domain with the boundaries of the flow domain being located inside the computational domain. The physical boundary conditions are enforced using the immersed boundary conditions concept. The algorithm relies on the Fourier expansions in the flow direction and the Chebyshev expansions in the transverse direction. Various tests confirm spectral accuracy of the algorithm. Copyright © 2010 John Wiley & Sons, Ltd.     

Coupling methods between finite element–based Boussinesq‐type wave and particle‐based free‐surface flow models

We develop one‐way coupling methods between a Boussinesq‐type wave model based on the discontinuous Galerkin finite element method and a free‐surface flow model based on a mesh‐free particle method to strike a balance between accuracy and computational cost. In our proposed model, computation of the wave model in the global domain is conducted first, and the nonconstant velocity profiles in the vertical direction are reproduced by using its results. Computation of the free‐surface flow is performed in a local domain included within the global domain with interface boundaries that move along the reproduced velocity field in a Lagrangian fashion. To represent the moving interfaces, we used a polygon wall boundary model for mesh‐free particle methods. Verification and validation tests of our proposed model are performed, and results obtained by the model are compared with theoretical values and experimental results to show its accuracy and applicability.     

Numerical simulation of compressible multifluid flows using an adaptive positivity-preserving RKDG-GFM approach

The Runge-Kutta discontinuous Galerkin method together with a refined real-ghost fluid method is incorporated into an adaptive mesh refinement environment for solving compressible multifluid flows, where the level set method is used to capture the moving material interface. To ensure that the Riemann problem is exactly along the normal direction of the material interface, a simple and efficient modification is introduced into the original real-ghost fluid method for constructing the interfacial Riemann problem, and the initial conditions of the Riemann problem are obtained directly from the solution polynomials of the discontinuous Galerkin finite element space. In addition, a positivity-preserving limiter is introduced into the Runge-Kutta discontinuous Galerkin method to suppress the failure of preserving positivity of density or pressure for the problems involving strong shock wave or shock interaction with material interface. For interfacial cells in adaptive mesh refinement, the data transfer between different grid levels is achieved by using a L² projection approach along with the least squares fitting. Various numerical cases, including multifluid shock tubes, underwater explosions, and shock-induced collapse of a underwater air bubble, are computed to assess the capability of the present adaptive positivity-preserving RKDG-GFM approach, and the simulated results show that the present approach is quite robust and can provide relatively reasonable results across a wide variety of flow regimes, even for problems involving strong shock wave or shock wave impacting high acoustic impedance mismatch material interface.     

Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluid–structure interaction problems

Embedded boundary methods for CFD (computational fluid dynamics) simplify a number of issues. These range from meshing the fluid domain, to designing and implementing Eulerian‐based algorithms for fluid–structure applications featuring large structural motions and/or deformations. Unfortunately, embedded boundary methods also complicate other issues such as the treatment of the wall boundary conditions in general, and fluid–structure transmission conditions in particular. This paper focuses on this aspect of the problem in the context of compressible flows, the finite volume method for the fluid, and the finite element method for the structure. First, it presents a numerical method for treating simultaneously the fluid pressure and velocity conditions on static and dynamic embedded interfaces. This method is based on the exact solution of local, one‐dimensional, fluid–structure Riemann problems. Next, it describes two consistent and conservative approaches for computing the flow‐induced loads on rigid and flexible embedded structures. The first approach reconstructs the interfaces within the CFD solver. The second one represents them as zero level sets, and works instead with surrogate fluid/structure interfaces. For example, the surrogate interfaces obtained simply by joining contiguous segments of the boundary surfaces of the fluid control volumes that are the closest to the zero level sets are explored in this work. All numerical algorithms presented in this paper are applicable with any embedding CFD mesh, whether it is structured or unstructured. Their performance is illustrated by their application to the solution of three‐dimensional fluid–structure interaction problems associated with the fields of aeronautics and underwater implosion. Copyright © 2011 John Wiley & Sons, Ltd.     

Numerical simulation of fully nonlinear irregular wave tank in three dimension

A fully nonlinear irregular wave tank has been developed using a three‐dimensional higher‐order boundary element method (HOBEM) in the time domain. The Laplace equation is solved at each time step by an integral equation method. Based on image theory, a new Green function is applied in the whole fluid domain so that only the incident surface and free surface are discretized for the integral equation. The fully nonlinear free surface boundary conditions are integrated with time to update the wave profile and boundary values on it by a semi‐mixed Eulerian–Lagrangian time marching scheme. The incident waves are generated by feeding analytic forms on the input boundary and a ramp function is introduced at the start of simulation to avoid the initial transient disturbance. The outgoing waves are sufficiently dissipated by using a spatially varying artificial damping on the free surface before they reach the downstream boundary. Numerous numerical simulations of linear and nonlinear waves are performed and the simulated results are compared with the theoretical input waves. Copyright © 2006 John Wiley & Sons, Ltd.     

Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing

The solution for the shallow water equations using smoothed particle hydrodynamics is attractive, being a mesh‐free, automatically adaptive method without special treatment for wet–dry interfaces. However, the relatively new method is limited by the variable kernel size or smoothing length being inversely proportional to water depth causing poor resolution at small depths. Boundary conditions at solid walls have also not been well resolved. To solve the resolution problem in small depths, a particle splitting procedure was developed (conveniently into seven particles), which conserves mass and momentum by varying the smoothing length, velocity and acceleration of each refined particle. This improves predictions in the shallowest depths where the error associated with splitting is reduced by one order of magnitude in comparison to other published works. To provide good shock capturing behaviour, particle interactions are treated as a Riemann problem with Monotone Upstream‐centred Scheme for Conservation Laws (MUSCL) reconstruction providing stability. For solid boundaries, the recent modified virtual boundary particle method was developed further to enable the zeroth moment to be accurately conserved where the smoothing length of particles is changing rapidly during particle splitting. The resulting method is applied to the one‐dimensional and the two‐dimensional axisymmetric wet‐bed dam break problems showing close agreement with analytical solutions, demonstrating the need for particle splitting. To demonstrate wetting and drying in a more complex case, the scheme is applied to oscillating water in a two‐dimensional parabolic basin and produces good agreement with the analytical solution. The method is finally applied to the European Concerted Action on DAm break Modelling dam‐break test case representative of realistic conditions and good predictions are made of experimental measurements with a 40% reduction in the computational time when particle splitting is employed. The overall method has thus become quite sophisticated but its generality and versatility will be attractive for various shallow water problems. Copyright © 2011 John Wiley & Sons, Ltd.     

Toward low‐noise synthetic turbulent inflow conditions for aeroacoustic calculations

The accuracy of boundary conditions for computational aeroacoustics is a well‐known challenge, due in part to the necessity of truncating the flow domain and replacing the analytical boundary conditions at infinity with numerical boundary conditions. In particular, the inflow boundary condition involving turbulent velocity or scalar fields is likely to introduce spurious waves into the domain, therefore degrading the flow behavior and deteriorating the physical acoustic waves. In this work, a method to generate low‐noise, divergence‐free, synthetic turbulence for inflow boundary conditions is proposed. It relies on the classical view of turbulence as a superposition of random eddies convected with the mean flow. Within the proposed model, the vector potential and the requirement that the individual eddies must satisfy the linearized momentum equations about the mean flow are used. The model is tested using isolated eddies convected through the inflow boundary and an experimental benchmark data for spatially decaying isotropic turbulence. Copyright © 2013 John Wiley & Sons, Ltd.     

Numerical simulation of vortical ideal fluid flow through curved channel

A numerical algorithm to study the boundary‐value problem in which the governing equations are the steady Euler equations and the vorticity is given on the inflow parts of the domain boundary is developed. The Euler equations are implemented in terms of the stream function and vorticity. An irregular physical domain is transformed into a rectangle in the computational domain and the Euler equations are rewritten with respect to a curvilinear co‐ordinate system. The convergence of the finite‐difference equations to the exact solution is shown experimentally for the test problems by comparing the computational results with the exact solutions on the sequence of grids. To find the pressure from the known vorticity and stream function, the Euler equations are utilized in the Gromeka–Lamb form. The numerical algorithm is illustrated with several examples of steady flow through a two‐dimensional channel with curved walls. The analysis of calculations shows strong dependence of the pressure field on the vorticity given at the inflow parts of the boundary. Plots of the flow structure and isobars, for different geometries of channel and for different values of vorticity on entrance, are also presented. Copyright © 2003 John Wiley & Sons, Ltd.     

Algorithm for analysis of flows in ribbed annuli

Spectral methods for analyses of steady flows in annuli bounded by walls with either axi‐symmetric or longitudinal ribs are developed. The physical boundary conditions are enforced using the immersed boundary conditions concept. In the former case, the Stokes stream function is used to eliminate pressure and to reduce system of field equations to a single fourth‐order partial differential equation. The ribs are assumed to be periodic in the axial direction and this permits representation of the solution in terms of the Fourier expansion. In the latter case, the problem is reduced to the Laplace equations for the flow modifications that can be expressed in terms of the Fourier expansions. The modal functions, which are functions of the radial coordinate, are discretized using Chebyshev polynomials. The problem formulations are closed using either the fixed volume flow rate constraint or the fixed pressure gradient constraint. Various tests have been carried out in order to demonstrate the spectral accuracy of the discretizations, as well as the spectral accuracy of the enforcement of the flow boundary conditions at the ribbed walls using the immersed boundary conditions concept. Special linear solver that takes advantage of the matrix structure has been implemented in order to reduce computational time and memory requirements. It is shown that the algorithm has superior performance when one is interested in the analysis of a large number of geometries, as part of the coefficient matrix that corresponds to the field equation is always the same and one needs to change only the part of the matrix that corresponds to the boundary relations when changing geometry of the flow domain. Copyright © 2011 John Wiley & Sons, Ltd.     

Shallow flow simulation on dynamically adaptive cut cell quadtree grids

A computationally efficient, high‐resolution numerical model of shallow flow hydrodynamics is described, based on dynamically adaptive quadtree grids. The numerical model solves the two‐dimensional non‐linear shallow water equations by means of an explicit second‐order MUSCL‐Hancock Godunov‐type finite volume scheme. Interface fluxes are evaluated using an HLLC approximate Riemann solver. Cartesian cut cells are used to improve the fit to curved boundaries. A ghost‐cell immersed boundary method is used to update flow information in the smallest cut cells and overcome the time step restriction that would otherwise apply. The numerical model is validated through simulations of reflection of a surge wave at a wall, a low Froude number potential flow past a circular cylinder, and the shock‐like interaction between a bore and a circular cylinder. The computational efficiency is shown to be greatly improved compared with solutions on a uniform structured grid implemented with cut cells. Copyright © 2006 John Wiley & Sons, Ltd.