Accurate and consistent particle tracking on unstructured grids

A new numerical method for particle tracking (Lagrangian particle advection) on 2‐D unstructured grids with triangular cells is presented and tested. This method combines key attributes of published methods, including streamline closure for steady flows and local mass conservation (uniformity preservation). The subgrid‐scale velocity reconstruction is linear, and this linear velocity field is integrated analytically to obtain particle trajectories. A complete analytic solution to the 2‐D system of ordinary differential equations (ODEs) governing particle trajectories within a grid cell is provided. The analytic solution to the linear system of locally mass‐conserving constraints that must be enforced to obtain the coefficients in the ODEs is also provided. Numerical experiments are performed to demonstrate that the new method has substantial advantages in accuracy over previously published methods and that it does not suffer from unphysical particle clustering. The method can be used not only in particle‐tracking applications but also as part of a semi‐Lagrangian advection scheme.Copyright © 2015 John Wiley & Sons, Ltd.     

A volume-tracking method for incompressible multifluid flows with large density variations

A numerical technique (FGVT) for solving the time-dependent incompressible Navier–Stokes equations in fluid flows with large density variations is presented for staggered grids. Mass conservation is based on a volume tracking method and incorporates a piecewise-linear interface reconstruction on a grid twice as fine as the velocity–pressure grid. It also uses a special flux-corrected transport algorithm for momentum advection, a multigrid algorithm for solving a pressure-correction equation and a surface tension algorithm that is robust and stable. In principle, the method conserves both mass and momentum exactly, and maintains extremely sharp fluid interfaces. Applications of the numerical method to prediction of two-dimensional bubble rise in an inclined channel and a bubble bursting through an interface are presented. © 1998 John Wiley & Sons, Ltd.     

A WAVE EQUATION MODEL TO SOLVE THE MULTIDIMENSIONAL TRANSPORT EQUATION

The wave equation model, originally developed to solve the advection–diffusion equation, is extended to the multidimensional transport equation in which the advection velocities vary in space and time. The size of the advection term with respect to the diffusion term is arbitrary. An operator-splitting method is adopted to solve the transport equation. The advection and diffusion equations are solved separate ly at each time step. During the advection phase the advection equation is solved using the wave equation model. Consistency of the first-order advection equation and the second-order wave equation is established. A finite element method with mass lumping is employed to calculate the three-dimensional advection of both a Gaussian cylinder and sphere in both translational and rotational flow fields. The numerical solutions are accurate in comparison with the exact solutions. The numerical results indicate that (i) the wave equation model introduces minimal numerical oscillation, (ii) mass lumping reduces the computational costs and does not significantly degrade the numerical solutions and (iii) the solution accuracy is relatively independent of the Courant number provided that a stability constraint is satisfied. © 1997 by John Wiley & Sons, Ltd.     

A non‐iterative implicit algorithm for the solution of advection–diffusion equation on a sphere

A numerical algorithm for the solution of advection–diffusion equation on the surface of a sphere is suggested. The velocity field on a sphere is assumed to be known and non‐divergent. The discretization of advection–diffusion equation in space is carried out with the help of the finite volume method, and the Gauss theorem is applied to each grid cell. For the discretization in time, the symmetrized double‐cycle componentwise splitting method and the Crank–Nicolson scheme are used. The numerical scheme is of second order approximation in space and time, correctly describes the balance of mass of substance in the forced and dissipative discrete system and is unconditionally stable. In the absence of external forcing and dissipation, the total mass and L₂‐norm of solution of discrete system is conserved in time. The one‐dimensional periodic problems arising at splitting in the longitudinal direction are solved with Sherman–Morrison's formula and Thomas's algorithm. The one‐dimensional problems arising at splitting in the latitudinal direction are solved by the bordering method that requires a prior determination of the solution at the poles. The resulting linear systems have tridiagonal matrices and are solved by Thomas's algorithm. The suggested method is direct (without iterations) and rapid in realization. It can also be applied to linear and nonlinear diffusion problems, some elliptic problems and adjoint advection–diffusion problems on a sphere. Copyright © 2015 John Wiley & Sons, Ltd.     

A stencil-like volume of fluid （VOF） method for tracking free interface

A stencil-like volume of fluid （VOF） method is proposed for tracking free interface. A stencil on a grid cell is worked out according to the normal direction of the interface, in which only three interface positions are possible in 2D cases, and the interface can be reconstructed by only requiring the known local volume fraction information. On the other hand, the fluid-occupying-length is defined on each side of the stencil, through which a unified fluid-occupying volume model and a unified algorithm can be obtained to solve the interface advection equation. The method is suitable for the arbitrary geometry of the grid cell, and is extendible to 3D cases. Typical numerical examples show that the current method can give ＂sharp＂ results for tracking free interface.     

Central finite volume methods for ideal and shallow water magnetohydrodynamics

We propose two‐dimensional central finite volume methods based on our multidimensional extensions of Nessyahu and Tadmor's one‐dimensional non‐oscillatory central scheme and a constrained transport‐type method to solve ideal magnetohydrodynamic problems (MHD) and shallow water magnetohydrodynamic problems (SMHD). The main numerical scheme is second‐order accurate both in space and time and uses an original Cartesian grid coupled to a Cartesian‐ or diamond‐staggered dual grid to by‐pass the resolution of the Riemann problems at the cell interfaces. To treat the non‐vanishing magnetic field/flux divergence we have constructed an adaptation of Evans and Hawley's constrained transport method specifically designed for central schemes. Our numerical results show the efficiency and the potential of the scheme. Copyright © 2010 John Wiley & Sons, Ltd.     

Simulating surface tension with smoothed particle hydrodynamics

A method for simulating two‐phase flows including surface tension is presented. The approach is based upon smoothed particle hydrodynamics (SPH). The fully Lagrangian nature of SPH maintains sharp fluid–fluid interfaces without employing high‐order advection schemes or explicit interface reconstruction. Several possible implementations of surface tension force are suggested and compared. The numerical stability of the method is investigated and optimal choices for numerical parameters are identified. Comparisons with a grid‐based volume of fluid method for two‐dimensional flows are excellent. The methods presented here apply to problems involving interfaces of arbitrary shape undergoing fragmentation and coalescence within a two‐phase system and readily extend to three‐dimensional problems. Boundary conditions at a solid surface, high viscosity and density ratios, and the simulation of free‐surface flows are not addressed. Copyright © 2000 John Wiley & Sons, Ltd.     

SOLUTION OF THE ADVECTION–DIFFUSION EQUATION USING A COMBINATION OF DISCONTINUOUS AND MIXED FINITE ELEMENTS

When transport is advection-dominated, classical numerical methods introduce excessive artificial diffusion and spurious oscillations. Special methods are required to overcome these phenomena. To solve the advection‒diffusion equation, a numerical method is developed using a discontinuous finite element method for the discretization of the advective terms. At the discontinuities of the approximate solution, numerical advective fluxes are calculated using one-dimensional approximate Riemann solvers. The method is stabilized with a multidimensional slope limiter which introduces small amounts of numerical diffusion when sharp concentration fronts occur. In addition, the diffusive term is discretized using a mixed hybrid finite element method. With this approach, numerical oscillations are completely avoided for a full range of cell Peclet numbers. The combination of discontinuous and mixed finite elements can be easily applied to 2D and 3D models using various types of elements in regular and irregular meshes. Numerical tests show good agreement with 1D and 2D analytical solutions. This approach is compared at the same time with two different numerical methods, a standard mixed finite method and a finite volume approach with high-resolution upwind terms. Regular and irregular meshes are used for the numerical tests to study the mesh effects on the numerical results. Our data show that in all cases this approach performs well. © 1997 by John Wiley & Sons, Ltd.     

Arnoldi and Crank–Nicolson methods for integration in time of the transport equation

A comparison is made between the Arnoldi reduction method and the Crank–Nicolson method for the integration in time of the advection–diffusion equation. This equation is first discretized in space by the classic finite element (FE) approach, leading to an unsymmetric first‐order differential system, which is then solved by the aforementioned methods. Arnoldi reduces the native FE equations to a much smaller set to be efficiently integrated in the Arnoldi vector space by the Crank–Nicolson scheme, with the solution recovered back by a standard Rayleigh–Ritz procedure. Crank–Nicolson implements a time marching scheme directly on the original first‐order differential system. The computational performance of both methods is investigated in two‐ and three‐dimensional sample problems with a size up 30 000. The results show that in advection‐dominated problems less then 100 Arnoldi vectors generally suffice to give results with a 10⁻³–10⁻⁴ difference relative to the direct Crank–Nicolson solution. However, while the CPU time with the Crank–Nicolson starts from zero and increases linearly with the number of time steps used in the simulation, the Arnoldi requires a large initial cost to generate the Arnoldi vectors with subsequently much less expensive dynamics for the time integration. The break‐even point is problem‐dependent at a number of time steps which may be for some problems up to one order of magnitude larger than the number of Arnoldi vectors. A serious limitation of Arnoldi is the requirement of linearity and time independence of the flow field. It is concluded that Arnoldi can be cheaper than Crank–Nicolson in very few instances, i.e. when the solution is needed for a large number of time values, say several hundreds or even 1000, depending on the problem. Copyright © 2001 John Wiley & Sons, Ltd.     

Tracking accuracy of a semi‐Lagrangian method for advection–dispersion modelling in rivers

There is an increasing need to improve the computational efficiency of river water quality models because: (1) Monte‐Carlo‐type multi‐simulation methods, that return solutions with statistical distributions or confidence intervals, are becoming the norm, and (2) the systems modelled are increasingly large and complex. So far, most models are based on Eulerian numerical schemes for advection, but these do not meet the requirement of efficiency, being restricted to Courant numbers below unity. The alternative of using semi‐Lagrangian methods, consisting of modelling advection by the method of characteristics, is free from any inherent Courant number restriction. However, it is subject to errors of tracking that result in potential phase errors in the solutions. The aim of this article is primarily to understand and estimate these tracking errors, assuming the use of a cell‐based backward method of characteristics, and considering conditions that would prevail in practical applications in rivers. This is achieved separately for non‐uniform flows and unsteady flows, either via theoretical considerations or using numerical experiments. The main conclusion is that, tracking errors are expected to be negligible in practical applications in both unsteady flows and non‐uniform flows. Also, a very significant computational time saving compared to Eulerian schemes is achievable. Copyright © 2006 John Wiley & Sons, Ltd.     

VOLUME-TRACKING METHODS FOR INTERFACIAL FLOW CALCULATIONS

A new algorithm for volume tracking which is based on the concept of flux-corrected transport (FCT) is introduced. It is applicable to incompressible 2D flow simulations on finite volume and difference meshes. The method requires no explicit interface reconstruction, is direction-split and can be extended to 3D and orthogonal curvilinear meshes in a straightforward manner. A comparison of the new scheme against well-known existing 2D finite volume techniques is undertaken. A series of progressively more difficult advection tests is used to test the accuracy of each scheme and it is seen that simple advection tests are inadequate indicators of the performance of volume-tracking methods. A straightforward methodology is presented that allows more rigorous estimates to be made of the error in volume advection and coupled volume and momentum advection in real flow situations. The volume advection schemes are put to a final test in the case of Rayleigh–Taylor instability. © 1997 by CSIRO.     

A higher-order eulerian scheme for coupled advection-diffusion transport

A new accurate high-order numerical method is presented for the coupled transport of a passive scalar (concentration) by advection and diffusion. Following the method of characteristics, the pure advection problem is first investigated. Interpolation of the concentration and its first derivative at the foot of the characteristic is carried out with a fifth-degree polynomial. The latter is constructed by using as information the concentration and its first and second derivatives at computational points on current time level t in Eulerian co-ordinates. The first derivative involved in the polynomial is transported by advection along the characteristic towards time level t + Δt in the same way as is the concentration itself. Second derivatives are obtained at the new time level t + Δt by solving a system of linear equations defined only by the concentrations and their derivatives at grid nodes, with the assumption that the third-order derivatives are continuous. The approximation of the method is of sixth order. The results are extended to coupled transport by advection and diffusion. Diffusion of the concentration takes place in parallel with advection along the characteristic. The applicability and precision of the method are demonstrated for the case of a Gaussian initial distribution of concentrations as well as for the case of a steep advancing concentration front. The results of the simulations are compared with analytical solutions and some existing methods.     

Consistent discretization for simulations of flows with moving generalized curvilinear coordinates

We develop a consistent discretization of conservative momentum and scalar transport for the numerical simulation of flow using a generalized moving curvilinear coordinate system. The formulation guarantees consistency between the discrete transport equation and the discrete mass conservation equation due to grid motion. This enables simulation of conservative transport using generalized curvilinear grids that move arbitrarily in three dimensions while maintaining the desired properties of the discrete transport equation on a stationary grid, such as constancy, conservation, and monotonicity. In addition to guaranteeing consistency for momentum and scalar transport, the formulation ensures geometric conservation and maintains the desired high‐order time accuracy of the discretization on a moving grid. Through numerical examples we show that, when the computation is carried out on a moving grid, consistency between the discretized scalar advection equation and the discretized equation for flow mass conservation due to grid motion is required in order to obtain stable and accurate results. We also demonstrate that significant errors can result when non‐consistent discretizations are employed. Copyright © 2009 John Wiley & Sons, Ltd.     

Semi‐Lagrangian advection on a spherical geodesic grid

A simple and efficient numerical method for solving the advection equation on the spherical surface is presented. To overcome the well‐known ‘pole problem’ related to the polar singularity of spherical coordinates, the space discretization is performed on a geodesic grid derived by a uniform triangulation of the sphere; the time discretization uses a semi‐Lagrangian approach. These two choices, efficiently combined in a substepping procedure, allow us to easily determine the departure points of the characteristic lines, avoiding any computationally expensive tree‐search. Moreover, suitable interpolation procedures on such geodesic grid are presented and compared. The performance of the method in terms of accuracy and efficiency is assessed on two standard test cases: solid‐body rotation and a deformation flow. Copyright © 2007 John Wiley & Sons, Ltd.     

A simple high‐resolution advection scheme

A simple, robust, mass‐conserving numerical scheme for solving the linear advection equation is described. The scheme can estimate peak solution values accurately even in regions where spatial gradients are high. Such situations present a severe challenge to classical numerical algorithms. Attention is restricted to the case of pure advection in one and two dimensions since this is where past numerical problems have arisen. The authors' scheme is of the Godunov type and is second‐order in space and time. The required cell interface fluxes are obtained by MUSCL interpolation and the exact solution of a degenerate Riemann problem. Second‐order accuracy in time is achieved via a Runge–Kutta predictor–corrector sequence. The scheme is explicit and expressed in finite volume form for ease of implementation on a boundary‐conforming grid. Benchmark test problems in one and two dimensions are used to illustrate the high‐spatial accuracy of the method and its applicability to non‐uniform grids. Copyright © 2007 John Wiley & Sons, Ltd.     

Some considerations for high‐order ‘incremental remap’‐based transport schemes: edges,reconstructions, and area integration

The problem of two‐dimensional tracer advection on the sphere is extremely important in modeling of geophysical fluids and has been tackled using a variety of approaches. A class of popular approaches for tracer advection include ‘incremental remap’ or cell‐integrated semi‐Lagrangian‐type schemes. These schemes achieve high‐order accuracy without the need for multistage integration in time, are capable of large time steps, and tend to be more efficient than other high‐order transport schemes when applied to a large number of tracers over a single velocity field. In this paper, the simplified flux‐form implementation of the Conservative Semi‐LAgrangian Multi‐tracer scheme (CSLAM) is reformulated using quadratic curves to approximate the upstream flux volumes and Gaussian quadrature for integrating the edge flux. The high‐order treatment of edge fluxes is motivated because of poor accuracy of the CSLAM scheme in the presence of strong nonlinear shear, such as one might observe in the midlatitudes near an atmospheric jet. Without the quadratic treatment of upstream edges, we observe at most second‐order accuracy under convergence of grid resolution, which is returned to third‐order accuracy under the improved treatment. A shallow‐water barotropic instability also reveals clear evidence of grid imprinting without the quadratic correction. Consequently, these tests reveal a problem that might arise in tracer transport near nonlinearly sheared regions of the real atmosphere, particularly near cubed‐sphere panel edges. Although CSLAM is used as the foundation for this analysis, the conclusions of this paper are applicable to the general class of incremental remap schemes. Copyright © 2012 John Wiley & Sons, Ltd.     

Lattice‐Boltzmann and meshless point collocation solvers for fluid flow and conjugate heat transfer

The applicability and performance of the lattice‐Boltzmann (LB) and meshless point collocation methods as CFD solvers in flow and conjugate heat transfer processes are investigated in this work. Lid‐driven cavity flow and flow in a slit with an obstacle including heat transfer are considered as case studies. A comparison of the computational efficiency accuracy of the two methods with that of a finite volume method as implemented in a commercial package (ANSYS CFX, ANSYS Inc., Canonsburg, PA) is made. Utilizing the analogy between heat and mass transfer, an advection–diffusion LB model was adopted to simulate the heat transfer part of the slit flow problem followed by a rigorous mapping of the mass transfer variables to the heat transfer quantities of interest, thus circumventing the need for a thermal LB model. Direct comparison among the results of the three methods revealed excellent agreement over a wide range of Reynolds and Prandtl number values. Furthermore, an integrated computational scheme is proposed, utilizing the rapid convergence of the LB model in the flow part of the conjugate heat transfer problem with that of the meshless collocation method for the heat transfer part. The meshless treatment remains sufficiently rapid even for conduction‐controlled processes in contrast to the LB method, which is very rapid in the convection‐controlled case only. A single, common computational grid, composed of regularly distributed nodes is used, saving significant computational and coding time and ensuring convergence of the discrete Laplacian operator in the heat transfer part of the computations. Copyright © 2012 John Wiley & Sons, Ltd.     

A conservative unstructured scheme for rapidly varied flows

This article introduces a new semi‐implicit, staggered finite volume scheme on unstructured meshes for the modelling of rapidly varied shallow water flows. Rapidly varied flows occur in the inundation of dry land during flooding situations. They typically involve bores and hydraulic jumps after obstacles such as road banks. Near such sudden flow transitions, the grid resolution is often low compared with the gradients of the bathymetry. Locally the hydrostatic pressure assumption may become invalid. In these situations, it is crucial to apply the correct conservation properties to obtain accurate results. An important feature of this scheme is therefore its ability to conserve momentum locally or, by choice, preserve constant energy head along a streamline. This is achieved using a special interpolation method and control volumes for momentum. The efficiency of inundation calculations with locally very high velocities, and in the case of unstructured meshes locally very small grid distances, is severely hampered by the Courant condition. This article provides a solution in the form of a locally implicit time integration for the advective terms that allows for an explicit calculation in most of the domain, while maintaining unconditional stability by implicit calculations only where necessary. The complex geometry of flooded urban areas asks for the flexibility of unstructured meshes. The efficient calculation of the pressure gradient in this, and other semi‐implicit staggered schemes, requires, however, an orthogonality condition to be put on the grid. In this article a simple method is introduced to generate unstructured hybrid meshes that fulfil this requirement. Copyright © 2008 John Wiley & Sons, Ltd.     

Capturing the cross-terms in multidimensional advection schemes

Solving the advection equation is an important part of numerically modeling the atmosphere. Both accuracy and efficiency are desirable traits of an advection scheme. For multidimensional flow, forward-in-time advection schemes must properly capture the cross-terms. Failure to capture the cross-terms can result in reduced accuracy and even instabilities. We show how multidimensional forward-in-time schemes successfully capture the cross-terms of two-dimensional (2D) flow. We then introduce a method to improve the efficiency of the forward-in-time schemes for 2D flow. This method stacks the duplicated cross-terms from one flux into the other, creating asymmetrized fluxes. Numerical testing shows that these asymmetrized flux calculation schemes perform to the same accuracy as the original forwards-in-time schemes but with a significant improvement in computational time. Finally, we show extensions of the method to three-dimensional (3D) flow.     

Numerical simulation of free surface flows with the level set method using an extremely high-order accuracy WENO advection scheme

Three-dimensional dynamic gas–liquid flow simulations that accurately track the phase interface are numerically challenging. This article presents a numerical study of the performance of the level-set phase–interface tracking method when combined with extremely high order (7th to 11th) weighted essentially non-oscillatory (WENO) advection schemes for gas–liquid free surface flows. Comparisons between simulation results and prior benchmark results suggest that such a combination of methods can be satisfactorily applied to the level-set and Navier-Stokes equations for free surface flow simulations when volume conservation is enforced at every time step, and minor numerical oscillations are suppressed through use of an artificial viscosity term. In particular, simulations of solid body rotation, the unsteady flow following an ideal dam break, tank sloshing, and the rise of a single bubble all agree with analytical or experimental results to within ± 3.12% when the level-set method is combined with an 11th order WENO scheme. Furthermore, use of an 11th order WENO advection scheme actually has a computational cost advantage because, for the same accuracy, it can be used on a coarser grid when compared with a more-common second-order advection scheme; computational savings of up to 87% are possible.