Accuracy and stability of a set of free‐surface time‐domain boundary element models based on B‐splines

An analysis is given for the accuracy and stability of some perturbation‐based time‐domain boundary element models (BEMs) with B‐spline basis functions, solving hydrodynamic free‐surface problems, including forward speed effects. The spatial convergence rate is found as a function of the order of the B‐spline basis. It is shown that for all the models examined the mixed implicit–explicit Euler time integration scheme is correct to second order. Stability diagrams are found for models based on B‐splines of orders third through to sixth for two different time integration schemes. The stability analysis can be regarded as an extension of the analysis by Vada and Nakos [Vada T, Nakos DE. Time marching schemes for ship motion simulations. In Proceedings of the 8th International Workshop on Water Waves and Floating Bodies, St. John's, Newfoundland, Canada, 1993; 155–158] to include B‐splines of orders higher than three (piecewise quadratic polynomials) and to include finite water depth and a current at an oblique angle to the model grid. Copyright © 2000 John Wiley & Sons, Ltd.     

Using a piecewise linear bottom to fit the bed variation in a laterally averaged,z‐co‐ordinate hydrodynamic model

In developing a 3D or laterally averaged 2D model for free‐surface flows using the finite difference method, the water depth is generally discretized either with the z‐co‐ordinate (z‐levels) or a transformed co‐ordinate (e.g. the so‐called σ‐co‐ordinate or σ‐levels). In a z‐level model, the water depth is discretized without any transformation, while in a σ‐level model, the water depth is discretized after a so‐called σ‐transformation that converts the water column to a unit, so that the free surface will be 0 (or 1) and the bottom will be ‐1 (or 0) in the stretched co‐ordinate system. Both discretization methods have their own advantages and drawbacks. It is generally not conclusive that one discretization method always works better than the other. The biggest problem for the z‐level model normally stems from the fact that it cannot fit the topography properly, while a σ‐level model does not have this kind of a topography‐fitting problem. To solve the topography‐fitting problem in a laterally averaged, 2D model using z‐levels, a piecewise linear bottom is proposed in this paper. Since the resulting computational cells are not necessarily rectangular looking at the x–z plane, flux‐based finite difference equations are used in the model to solve the governing equations. In addition to the piecewise linear bottom, the model can also be run with full cells or partial cells (both full cell and partial cell options yield a staircase bottom that does not fit the real bottom topography). Two frictionless wave cases were chosen to evaluate the responses of the model to different treatments of the topography. One wave case is a boundary value problem, while the other is an initial value problem. To verify that the piecewise linear bottom does not cause increased diffusions for areas with steep bottom slopes, a barotropic case in a symmetric triangular basin was tested. The model was also applied to a real estuary using various topography treatments. The model results demonstrate that fitting the topography is important for the initial value problem. For the boundary value problem, topography‐fitting may not be very critical if the vertical spacing is appropriate. Copyright © 2004 John Wiley & Sons, Ltd.     

A three‐dimensional hydrodynamic model for shallow waters using unstructured Cartesian grids

This paper presents a three‐dimensional unstructured Cartesian grid model for simulating shallow water hydrodynamics in lakes, rivers, estuaries, and coastal waters. It is a flux‐based finite difference model that uses a cut‐cell approach to fit the bottom topography and shorelines and, at the same time, has the flexibility of discretizing complex geometries with Cartesian grids that can be arbitrarily downsized in the two horizontal directions simultaneously. Because of the use of Cartesian grids, the grid generation is very simple and does not suffer the grid generation headache often seen in many other unstructured models, as the unstructured Cartesian grid model does not have any requirements on the orthogonality of the grids. The newly developed unstructured Cartesian grid model was validated against analytical solutions for a three‐dimensional seiching case in a rectangular basin, before it was compared with another three‐dimensional model named LESS3D for circulations and salinity transport processes in an idealized embayment that is driven by tides and freshwater inflows. Model tests show that the numerical procedure used in the unstructured Cartesian grid model is robust. Similar to other unstructured models, a variable grid size has resulted in a smaller number of grids required for a reasonable model simulation, which in turn reduces the CPU time used in the model run. Copyright © 2010 John Wiley & Sons, Ltd.     

A well‐balanced explicit/semi‐implicit finite element scheme for shallow water equations in drying–wetting areas

Numerical solutions of the shallow water equations can be used to reproduce flow hydrodynamics occurring in a wide range of regions. In hydraulic engineering, the objectives include the prediction of dam break wave propagation, fluvial floods and other catastrophic flooding phenomena, the modeling of estuarine and coastal circulations, and the design and optimization of hydraulic structures. In this paper, a well‐balanced explicit and semi‐implicit finite element scheme for shallow water equations over complex domains involving wetting and drying is proposed. The governing equations are discretized by a fractional finite element method using a two‐step Taylor–Galerkin scheme. First, the intermediate increment of conserved variable is obtained explicitly neglecting the pressure gradient term. This is then corrected for the effects of pressure once the pressure increment has been obtained from the Poisson equation. In order to maintain the ‘well‐balanced’ property, the pressure gradient term and bed slope terms are incorporated into the Poisson equation. Moreover, a local bed slope modification technique is employed in drying–wetting interface treatments. The proposed model is well validated against several theoretical benchmark tests. Copyright © 2014 John Wiley & Sons, Ltd.     

A mass‐fraction‐based interface‐capturing method for multi‐component flow

An interface‐capturing method based on mass fraction is developed to solve the Riemann problem in multi‐component compressible flow. Equations of mass fraction with modified form, which is derived from conservative equations of mass, are employed here to capture the interface. By introducing mass fraction into Euler equations system, as well as other conservative coefficients, a quasi‐conservative numerical model is created. Numerical examples show that the mass fraction model performs well not only in multi‐component fluids modeled by simple stiffened gas equation of state (EOS) but also in that modeled by complex Mie–Grüneisen EOS. Moreover, the mass fraction model is applied to Riemann problem with piecewise EOS; the expression of which depends on density. It is found that the mass fraction model can well adapt to the analytic change in piecewise EOS and produce accuracy solutions with fewer unknown quantities, and the model can be easily extended to m‐component fluid mixture by using only m + 4 equations with no additional conditions. Copyright © 2013 John Wiley & Sons, Ltd.     

A new well‐balanced Hermite weighted essentially non‐oscillatory scheme for shallow water equations

Hermite weighted essentially non‐oscillatory (HWENO) methods were introduced in the literature, in the context of Euler equations for gas dynamics, to obtain high‐order accuracy schemes characterized by high compactness (e.g. Qiu and Shu, J. Comput. Phys. 2003; 193 :115). For example, classical fifth‐order weighted essentially non‐oscillatory (WENO) reconstructions are based on a five‐cell stencil whereas the corresponding HWENO reconstructions are based on a narrower three‐cell stencil. The compactness of the schemes allows easier treatment of the boundary conditions and of the internal interfaces. To obtain this compactness in HWENO schemes both the conservative variables and their first derivatives are evolved in time, whereas in the original WENO schemes only the conservative variables are evolved. In this work, an HWENO method is applied for the first time to the shallow water equations (SWEs), including the source term due to the bottom slope, to obtain a fourth‐order accurate well‐balanced compact scheme. Time integration is performed by a strong stability preserving the Runge–Kutta method, which is a five‐step and fourth‐order accurate method. Besides the classical SWE, the non‐homogeneous equations describing the time and space evolution of the conservative variable derivatives are considered here. An original, well‐balanced treatment of the source term involved in such equations is developed and tested. Several standard one‐dimensional test cases are used to verify the high‐order accuracy, the C‐property and the good resolution properties of the model. Copyright © 2010 John Wiley & Sons, Ltd.     

The consistency of pressure‐gradient approximations used in multi‐dimensional shock hydrodynamics

This paper presents a study of the consistency properties of the pressure‐gradient approximation used in multi‐dimensional finite‐element shock hydrodynamics codes today. In specific, consideration is given to the so‐called ‘bent‐element blues’ problem associated with the pressure‐gradient approximation when using the Q1Q0 element. On arbitrary grids comprised of distorted elements, the piecewise‐constant representation of the pressure field leads to a low‐order pressure‐gradient approximation at the global (nodal) level. This results in spurious nodal forces that are not aligned with the pressure gradient. There are several side‐effects of this behavior that include (a) incorrectly exciting physical modes in problems that exhibit unstable behavior, e.g. Rayleigh–Taylor problems (both magnetic and hydrodynamic), (b) potentially seeding hourglass modes, and (c) exhibiting non‐stationary behavior for steady‐state problems. A series of commonly used pressure‐gradient approximations are reviewed and evaluated based on linear consistency—the ability of the approximation to annihilate constant terms and exactly reproduce a linear gradient. The deeper theoretical issues associated with the proper selection of function spaces for the finite‐element hydro formulation are not discussed here. There are two gradient approximations that use piecewise‐constant data and deliver a consistent pressure‐gradient approximation on arbitrary grids. The first is the well‐known least‐squares gradient construction, and the second is a corrected gradient approximation that imposes linear consistency at the (global) nodal level. At the time of this writing, the corrected gradient approximation appears to be the most viable candidate for resolving the consistency issues associated with the Q1Q0 element technology. Copyright © 2009 John Wiley & Sons, Ltd.     

Efficient second‐order time integration for single‐species aerosol formation and evolution

The dynamics of a single‐species aerosol composed of droplets in air is described in terms of nucleation, evaporation, condensation, and coagulation processes. We present a comprehensive overview of the Euler–Euler formulation, which gives rise to a model in which fast nucleation that initiates aerosol droplets co‐exists with comparably slow condensation. The latter process is responsible for the subsequent growth of the droplets. To accurately represent the dynamical consequences of the fast nucleation process, while retaining numerical efficiency, a new second‐order time‐integration method for the nucleation, evaporation, and condensation processes is proposed and analyzed. The new time‐integration method takes the form of a ‘corrected Euler forward’ method. It includes rapid nucleation bursts and their possible cessation within a time step Δt. If the current nucleation burst persists for longer than the next time step, it is included fully, whereas cessation of the nucleation burst within the next Δt implies corrections to the effective rates in the algorithm. The identification of these two situations corresponds to the physical mechanism by which nucleation of a supersaturated vapor is halted because of the progressing condensation onto the already formed droplets. The resulting time‐integration method is shown to be second‐order accurate in time, whereas the computational costs per time step were found to be increased by less than 25% compared with the Euler forward method. The new method is also applied in combination with advective transport of the aerosol forming vapor to investigate a front of rapid nucleation. Adopting robust first‐order upwinding for the spatial discretization, we arrive at a flexible method that shows an overall first‐order convergence in Δt. For the full, spatially dependent system motivated by an aerosol of water droplets in air, the computational benefits of the new time‐integration method over the Euler forward scheme, are a factor of about 10 improvements in accuracy at a given Δt and a similar factor in computing time when keeping the same level of accuracy. Copyright © 2013 John Wiley & Sons, Ltd.     

A robust and well‐balanced scheme for the 2D Saint‐Venant system on unstructured meshes with friction source term

In the following lines, we propose a numerical scheme for the shallow‐water system supplemented by topography and friction source terms, in a 2D unstructured context. This work proposes an improved version of the well‐balanced and robust numerical model recently introduced by Duran et al. (J. Comp. Phys., 235 , 565–586, 2013) for the pre‐balanced shallow‐water equations, accounting for varying topography. The present work aims at relaxing the robustness condition and includes a friction term. To this purpose, the scheme is modified using a recent method, entirely based on a modified Riemann solver. This approach preserves the robustness and well‐balanced properties of the original scheme and prevents unstable computations in the presence of low water depths. A series of numerical experiments are devoted to highlighting the performances of the resulting scheme. Simulations involving dry areas, complex geometry and topography are proposed to validate the stability of the numerical model in the neighbourhood of wet/dry transitions. Copyright © 2015 John Wiley & Sons, Ltd.     

Three‐dimensional modeling of non‐hydrostatic free‐surface flows on unstructured grids

A three‐dimensional numerical model is presented for the simulation of unsteady non‐hydrostatic shallow water flows on unstructured grids using the finite volume method. The free surface variations are modeled by a characteristics‐based scheme, which simulates sub‐critical and super‐critical flows. Three‐dimensional velocity components are considered in a collocated arrangement with a σ‐coordinate system. A special treatment of the pressure term is developed to avoid the water surface oscillations. Convective and diffusive terms are approximated explicitly, and an implicit discretization is used for the pressure term to ensure exact mass conservation. The unstructured grid in the horizontal direction and the σ coordinate in the vertical direction facilitate the use of the model in complicated geometries. Solution of the non‐hydrostatic equations enables the model to simulate short‐period waves and vertically circulating flows. Copyright © 2015 John Wiley & Sons, Ltd.     

Linear Rayleigh-Taylor instability analysis of double-shell Kidder＇s self-similar implosion solution

This paper generalizes the single-shell Kidder＇s self-similar solution to the double-shell one with a discontinuity in density across the interface. An isentropic implosion model is constructed to study the Rayleigh-Taylor instability for the implosion compression. A Godunov-type method in the Lagrangian coordinates is used to compute the one-dimensional Euler equation with the initial and boundary conditions for the double-shell Kidder＇s self-similar solution in spherical geometry. Numerical results are obtained to validate the double-shell implosion model. By programming and using the linear perturbation codes, a linear stability analysis on the Rayleigh-Taylor instability for the double-shell isentropic implosion model is performed. It is found that, when the initial perturbation is concentrated much closer to the interface of the two shells, or when the spherical wave number becomes much smaller, the modal radius of the interface grows much faster, i.e., more unstable. In addition, from the spatial point of view for the compressibility effect on the perturbation evolution, the compressibility of the outer shell has a destabilization effect on the Rayleigh-Taylor instability, while the compressibility of the inner shell has a stabilization effect.     

Evaluation of Smagorinsky‐based subgrid‐scale models in a finite‐volume computation

Smagorinsky‐based models are assessed in a turbulent channel flow simulation at Re_b=2800 and Re_b=12500. The Navier–Stokes equations are solved with three different grid resolutions by using a co‐located finite‐volume method. Computations are repeated with Smagorinsky‐based subgrid‐scale models. A traditional Smagorinsky model is implemented with a van Driest damping function. A dynamic model assumes a similarity of the subgrid and the subtest Reynolds stresses and an explicit filtering operation is required. A top‐hat test filter is implemented with a trapezoidal and a Simpson rule. At the low Reynolds number computation none of the tested models improves the results at any grid level compared to the calculations with no model. The effect of the subgrid‐scale model is reduced as the grid is refined. The numerical implementation of the test filter influences on the result. At the higher Reynolds number the subgrid‐scale models stabilize the computation. An analysis of an accurately resolved flow field reveals that the discretization error overwhelms the subgrid term at Re_b=2800 in the most part of the computational domain. Copyright © 2002 John Wiley & Sons, Ltd.     

Lattice Boltzmann modeling of pore‐scale fluid flow through idealized porous media

In order to capture the hydro‐mechanical impacts on the solid skeleton imposed by the fluid flowing through porous media at the pore‐scale, the flow in the pore space has to be modeled at a resolution finer than the pores, and the no‐slip condition needs to be enforced at the grain–fluid interface. In this paper, the lattice Boltzmann method (LBM), a mesoscopic Navier–Stokes solver, is shown to be an appropriate pore‐scale fluid flow model. The accuracy and lattice sensitivity of LBM as a fluid dynamics solver is demonstrated in the Poiseuille channel flow problem (2‐D) and duct flow problem (3‐D). Well‐studied problems of fluid creeping through idealized 2‐D and 3‐D porous media (J. Fluid Mech. 1959; 5 (2):317–328, J. Fluid Mech. 1982; 115 :13–26, Int. J. Multiphase Flow 1982; 8 (4):343–360, Phys. Fluids A 1989; 1 (1):38–46, Int. J. Numer. Anal. Meth. Geomech. 1999; 23 :881–904, Int. J. Numer. Anal. Meth. Geomech. 2010; DOI: 10.1002/nag.898, Int. J. Multiphase Flow 1982; 8 (3):193–206) are then simulated using LBM to measure the friction coefficient for various pore throats. The simulation results agree well with the data reported in the literature. The lattice sensitivity of the frictional coefficient is also investigated. Copyright © 2010 John Wiley & Sons, Ltd.     

SPH computation of plunging waves using a 2‐D sub‐particle scale (SPS) turbulence model

The paper presents a 2‐D large eddy simulation (LES) modelling approach to investigate the properties of the plunging waves. The numerical model is based on the smoothed particle hydrodynamics (SPH) method. SPH is a mesh‐free Lagrangian particle approach which is capable of tracking the free surfaces of large deformation in an easy and accurate way. The Smagorinsky model is used as the turbulence model due to its simplicity and effectiveness. The proposed 2‐D SPH–LES model is applied to a cnoidal wave breaking and plunging over a mild slope. The computations are in good agreement with the documented data. Especially the computed turbulence quantities under the breaking waves agree better with the experiments as compared with the numerical results obtained by using the k–ε model. The sensitivity analyses of the SPH–LES computations indicate that both the turbulence model and the spatial resolution play an important role in the model predictions and the contributions from the sub‐particle scale (SPS) turbulence decrease with the particle size refinement. Copyright © 2006 John Wiley & Sons, Ltd.     

A coupled surface‐subsurface model for hydrostatic flows under saturated and variably saturated conditions

In this paper, the governing differential equations for hydrostatic surface‐subsurface flows are derived from the Richards and from the Navier‐Stokes equations. A vertically integrated continuity equation is formulated to account for both surface and subsurface flows under saturated and variable saturated conditions. Numerically, the horizontal domain is covered by an unstructured orthogonal grid that may include subgrid specifications. Along the vertical direction, a simple z‐layer discretization is adopted. Semi‐implicit finite difference equations for velocities, and a finite volume approximation for the vertically integrated continuity equation, are derived in such a fashion that, after simple manipulation, the resulting discrete pressure equation can be assembled into a single, two‐dimensional, mildly nonlinear system. This system is solved by a nested Newton‐type method, which yields simultaneously the (hydrostatic) pressure and a nonnegative fluid volume throughout the computational grid. The resulting algorithm is relatively simple, extremely efficient, and very accurate. Stability, convergence, and exact mass conservation are assured throughout also in presence of wetting and drying, in variable saturated conditions, and during flow transition through the soil interface. A few examples illustrate the model applicability and demonstrate the effectiveness of the proposed algorithm.     

Numerical flux functions for Reynolds‐averaged Navier–Stokes and kω turbulence model computations with a line‐preconditioned p‐multigrid discontinuous Galerkin solver

We present an eigen‐decomposition of the quasi‐linear convective flux formulation of the completely coupled Reynolds‐averaged Navier–Stokes and kω turbulence model equations. Based on these results, we formulate different approximate Riemann solvers that can be used as numerical flux functions in a DG discretization. The effect of the different strategies on the solution accuracy is investigated with numerical examples. The actual computations are performed using a p‐multigrid algorithm. To this end, we formulate a framework with a backward‐Euler smoother in which the linear systems are solved with a general preconditioned Krylov method. We present matrix‐free implementations and memory‐lean line‐Jacobi preconditioners and compare the effects of some parameter choices. In particular, p‐multigrid is found to be less efficient than might be expected from recent findings by other authors. This might be due to the consideration of turbulent flow. Copyright © 2012 John Wiley & Sons, Ltd.     

A multidimensional HLL‐Riemann solver for non‐linear hyperbolic systems

This article presents a numerical model that enables to solve on unstructured triangular meshes and with a high order of accuracy, Riemann problems that appear when solving hyperbolic systems. For this purpose, we use a MUSCL‐like procedure in a ‘cell‐vertex’ finite‐volume framework. In the first part of this procedure, we devise a four‐state bi‐dimensional HLL solver (HLL‐2D). This solver is based upon the Riemann problem generated at the barycenter of a triangular cell, from the surrounding cell‐averages. A new three‐wave model makes it possible to solve this problem, approximately. A first‐order version of the bi‐dimensional Riemann solver is then generated for discretizing the full compressible Euler equations. In the second part of the MUSCL procedure, we develop a polynomial reconstruction that uses all the surrounding numerical data of a given point, to give at best third‐order accuracy. The resulting over determined system is solved by using a least‐square methodology. To enforce monotonicity conditions into the polynomial interpolation, we use and adapt the monotonicity‐preserving limiter, initially devised by Barth (AIAA Paper 90‐0013, 1990). Numerical tests and comparisons with competing numerical methods enable to identify the salient features of the whole model. Copyright © 2010 John Wiley & Sons, Ltd.     

Further experiences with computing non‐hydrostatic free‐surface flows involving water waves

A semi‐implicit, staggered finite volume technique for non‐hydrostatic, free‐surface flow governed by the incompressible Euler equations is presented that has a proper balance between accuracy, robustness and computing time. The procedure is intended to be used for predicting wave propagation in coastal areas. The splitting of the pressure into hydrostatic and non‐hydrostatic components is utilized. To ease the task of discretization and to enhance the accuracy of the scheme, a vertical boundary‐fitted co‐ordinate system is employed, permitting more resolution near the bottom as well as near the free surface. The issue of the implementation of boundary conditions is addressed. As recently proposed by the present authors, the Keller‐box scheme for accurate approximation of frequency wave dispersion requiring a limited vertical resolution is incorporated. The both locally and globally mass conserved solution is achieved with the aid of a projection method in the discrete sense. An efficient preconditioned Krylov subspace technique to solve the discretized Poisson equation for pressure correction with an unsymmetric matrix is treated. Some numerical experiments to show the accuracy, robustness and efficiency of the proposed method are presented. Copyright © 2004 John Wiley & Sons, Ltd.     

Baroclinic stability for a family of two‐level,semi‐implicit numerical methods for the 3D shallow water equations

The baroclinic stability of a family of two time‐level, semi‐implicit schemes for the 3D hydrostatic, Boussinesq Navier–Stokes equations (i.e. the shallow water equations), which originate from the TRIM model of Casulli and Cheng (Int. J. Numer. Methods Fluids 1992; 15 :629–648), is examined in a simple 2D horizontal–vertical domain. It is demonstrated that existing mass‐conservative low‐dissipation semi‐implicit methods, which are unconditionally stable in the inviscid limit for barotropic flows, are unstable in the same limit for baroclinic flows. Such methods can be made baroclinically stable when the integrated continuity equation is discretized with a barotropically dissipative backwards Euler scheme. A general family of two‐step predictor‐corrector schemes is proposed that have better theoretical characteristics than existing single‐step schemes. Copyright © 2006 John Wiley & Sons, Ltd.     

Evaluation of well‐balanced bore‐capturing schemes for 2D wetting and drying processes

We consider numerical solutions of the two‐dimensional non‐linear shallow water equations with a bed slope source term. These equations are well‐suited for the study of many geophysical phenomena, including coastal engineering where wetting and drying processes are commonly observed. To accurately describe the evolution of moving shorelines over strongly varying topography, we first investigate two well‐balanced methods of Godunov‐type, relying on the resolution of non‐homogeneous Riemann problems. But even if these schemes were previously proved to be efficient in many simulations involving occurrences of dry zones, they fail to compute accurately moving shorelines. From this, we investigate a new model, called SURF_WB, especially designed for the simulation of wave transformations over strongly varying topography. This model relies on a recent reconstruction method for the treatment of the bed‐slope source term and is able to handle strong variations of topography and to preserve the steady states at rest. In addition, the use of the recent VFRoe‐ncv Riemann solver leads to a robust treatment of wetting and drying phenomena. An adapted ‘second order’ reconstruction generates accurate bore‐capturing abilities.This scheme is validated against several analytical solutions, involving varying topography, time dependent moving shorelines and convergences toward steady states. This model should have an impact in the prediction of 2D moving shorelines over strongly irregular topography. Copyright © 2006 John Wiley & Sons, Ltd.