Vorticity–velocity formulation of the 3D Navier–Stokes equations in cylindrical co‐ordinates

A finite difference method is presented for solving the 3D Navier–Stokes equations in vorticity–velocity form. The method involves solving the vorticity transport equations in ‘curl‐form’ along with a set of Cauchy–Riemann type equations for the velocity. The equations are formulated in cylindrical co‐ordinates and discretized using a staggered grid arrangement. The discretized Cauchy–Riemann type equations are overdetermined and their solution is accomplished by employing a conjugate gradient method on the normal equations. The vorticity transport equations are solved in time using a semi‐implicit Crank–Nicolson/Adams–Bashforth scheme combined with a second‐order accurate spatial discretization scheme. Special emphasis is put on the treatment of the polar singularity. Numerical results of axisymmetric as well as non‐axisymmetric flows in a pipe and in a closed cylinder are presented. Comparison with measurements are carried out for the axisymmetric flow cases. Copyright © 2003 John Wiley & Sons, Ltd.     

Numerical implementation of the Crank–Nicolson/Adams–Bashforth scheme for the time‐dependent Navier–Stokes equations

This article considers numerical implementation of the Crank–Nicolson/Adams–Bashforth scheme for the two‐dimensional non‐stationary Navier–Stokes equations. A finite element method is applied for the spatial approximation of the velocity and pressure. The time discretization is based on the Crank–Nicolson scheme for the linear term and the explicit Adams–Bashforth scheme for the nonlinear term. Comparison with other methods, through a series of numerical experiments, shows that this method is almost unconditionally stable and convergent, i.e. stable and convergent when the time step is smaller than a given constant. Copyright © 2009 John Wiley & Sons, Ltd.     

Efficient discretization of pressure‐correction equations on non‐orthogonal grids

An alternative discretization of pressure‐correction equations within pressure‐correction schemes for the solution of the incompressible Navier–Stokes equations is introduced, which improves the convergence and robustness properties of such schemes for non‐orthogonal grids. As against standard approaches, where the non‐orthogonal terms usually are just neglected, the approach allows for a simplification of the pressure‐correction equation to correspond to 5‐point or 7‐point computational molecules in two or three dimensions, respectively, but still incorporates the effects of non‐orthogonality. As a result a wide range (including rather high values) of underrelaxation factors can be used, resulting in an increased overall performance of the underlying pressure‐correction schemes. Within this context, a second issue of the paper is the investigation of the accuracy to which the pressure‐correction equation should be solved in each pressure‐correction iteration. The scheme is investigated for standard test cases and, in order to show its applicability to practical flow problems, for a more complex configuration of a micro heat exchanger. Copyright © 2003 John Wiley & Sons, Ltd.     

Augmented upwind numerical schemes for the groundwater transport advection‐dispersion equation with local operators

When solute transport is advection‐dominated, the advection‐dispersion equation approximates to a hyperbolic‐type partial differential equation, and finite difference and finite element numerical approximation methods become prone to artificial oscillations. The upwind scheme serves to correct these responses to produce a more realistic solution. The upwind scheme is reviewed and then applied to the advection‐dispersion equation with local operators for the first‐order upwinding numerical approximation scheme. The traditional explicit and implicit schemes, as well as the Crank‐Nicolson scheme, are developed and analyzed for numerical stability to form a comparison base. Two new numerical approximation schemes are then proposed, namely, upwind–Crank‐Nicolson scheme, where only for the advection term is applied, and weighted upwind‐downwind scheme. These newly developed schemes are analyzed for numerical stability and compared to the traditional schemes. It was found that an upwind–Crank‐Nicolson scheme is appropriate if the Crank‐Nicolson scheme is only applied to the advection term of the advection‐dispersion equation. Furthermore, the proposed explicit weighted upwind‐downwind finite difference numerical scheme is an improvement on the traditional explicit first‐order upwind scheme, whereas the implicit weighted first‐order upwind‐downwind finite difference numerical scheme is stable under all assumptions when the appropriate weighting factor (θ) is assigned.     

An implicit three‐dimensional fully non‐hydrostatic model for free‐surface flows

An implicit method is developed for solving the complete three‐dimensional (3D) Navier–Stokes equations. The algorithm is based upon a staggered finite difference Crank‐Nicholson scheme on a Cartesian grid. A new top‐layer pressure treatment and a partial cell bottom treatment are introduced so that the 3D model is fully non‐hydrostatic and is free of any hydrostatic assumption. A domain decomposition method is used to segregate the resulting 3D matrix system into a series of two‐dimensional vertical plane problems, for each of which a block tri‐diagonal system can be directly solved for the unknown horizontal velocity. Numerical tests including linear standing waves, nonlinear sloshing motions, and progressive wave interactions with uneven bottoms are performed. It is found that the model is capable to simulate accurately a range of free‐surface flow problems using a very small number of vertical layers (e.g. two–four layers). The developed model is second‐order accuracy in time and space and is unconditionally stable; and it can be effectively used to model 3D surface wave motions. Copyright © 2004 John Wiley & Sons, Ltd.     

Finite element and implicit Runge–Kutta implementation of an acoustics–convection upstream resolution algorithm for the time‐dependent two‐dimensional Euler equations

The second of a two‐paper series, this paper details a solver for the characteristics‐bias system from the acoustics–convection upstream resolution algorithm for the Euler and Navier–Stokes equations. An integral formulation leads to several surface integrals that allow effective enforcement of boundary conditions. Also presented is a new multi‐dimensional procedure to enforce a pressure boundary condition at a subsonic outlet, a procedure that remains accurate and stable. A classical finite element Galerkin discretization of the integral formulation on any prescribed grid directly yields an optimal discretely conservative upstream approximation for the Euler and Navier–Stokes equations, an approximation that remains multi‐dimensional independently of the orientation of the reference axes and computational cells. The time‐dependent discrete equations are then integrated in time via an implicit Runge–Kutta procedure that in this paper is proven to remain absolutely non‐linearly stable for the spatially‐discrete Euler and Navier–Stokes equations and shown to converge rapidly to steady states, with maximum Courant number exceeding 100 for the linearized version. Even on relatively coarse grids, the acoustics–convection upstream resolution algorithm generates essentially non‐oscillatory solutions for subsonic, transonic and supersonic flows, encompassing oblique‐ and interacting‐shock fields that converge within 40 time steps and reflect reference exact solutions. Copyright © 2005 John Wiley & Sons, Ltd.     

Calculation of curved open channel flow using physical curvilinear non‐orthogonal co‐ordinates

There are two main difficulties in numerical simulation calculations using FD/FV method for the flows in real rivers. Firstly, the boundaries are very complex and secondly, the generated grid is usually very non‐uniform locally. Some numerical models in this field solve the first difficulty by the use of physical curvilinear orthogonal co‐ordinates. However, it is very difficult to generate an orthogonal grid for real rivers and the orthogonal restriction often forces the grid to be over concentrated where high resolution is not required. Recently, more and more models solve the first difficulty by the use of generalized curvilinear co‐ordinates (ξ,η). The governing equations are expressed in a covariant or contra‐variant form in terms of generalized curvilinearco‐ordinates (ξ,η). However, some studies in real rivers indicate that this kind of method has some undesirable mesh sensitivities. Sharp differences in adjacent mesh size may easily lead to a calculation stability problem oreven a false simulation result. Both approaches used presently have their own disadvantages in solving the two difficulties that exist in real rivers. In this paper, the authors present a method for two‐dimensional shallow water flow calculations to solve both of the main difficulties, by formulating the governing equations in a physical form in terms of physical curvilinear non‐orthogonal co‐ordinates (s,n). Derivation of the governing equations is explained, and two numerical examples are employed to demonstrate that the presented method is applicable to non‐orthogonal and significantly non‐uniform grids. Copyright © 2004 John Wiley & Sons, Ltd.     

A robust multi‐grid pressure‐based algorithm for multi‐fluid flow at all speeds

This paper reports on the implementation and testing, within a full non‐linear multi‐grid environment, of a new pressure‐based algorithm for the prediction of multi‐fluid flow at all speeds. The algorithm is part of the mass conservation‐based algorithms (MCBA) group in which the pressure correction equation is derived from overall mass conservation. The performance of the new method is assessed by solving a series of two‐dimensional two‐fluid flow test problems varying from turbulent low Mach number to supersonic flows, and from very low to high fluid density ratios. Solutions are generated for several grid sizes using the single grid (SG), the prolongation grid (PG), and the full non‐linear multi‐grid (FMG) methods. The main outcomes of this study are: (i) a clear demonstration of the ability of the FMG method to tackle the added non‐linearity of multi‐fluid flows, which is manifested through the performance jump observed when using the non‐linear multi‐grid approach as compared to the SG and PG methods; (ii) the extension of the FMG method to predict turbulent multi‐fluid flows at all speeds. The convergence history plots and CPU‐times presented indicate that the FMG method is far more efficient than the PG method and accelerates the convergence rate over the SG method, for the problems solved and the grids used, by a factor reaching a value as high as 15. Copyright © 2003 John Wiley & Sons, Ltd.     

A high‐order alternating direction implicit method for the unsteady convection‐dominated diffusion problem

A high‐order alternating direction implicit (ADI) method for solving the unsteady convection‐dominated diffusion equation is developed. The fourth‐order Padé scheme is used for the discretization of the convection terms, while the second‐order Padé scheme is used for the diffusion terms. The Crank–Nicolson scheme and ADI factorization are applied for time integration. After ADI factorization, the two‐dimensional problem becomes a sequence of one‐dimensional problems. The solution procedure consists of multiple use of a one‐dimensional tridiagonal matrix algorithm that produces a computationally cost‐effective solver. Von Neumann stability analysis is performed to show that the method is unconditionally stable. An unsteady two‐dimensional problem concerning convection‐dominated propagation of a Gaussian pulse is studied to test its numerical accuracy and compare it to other high‐order ADI methods. The results show that the overall numerical accuracy can reach third or fourth order for the convection‐dominated diffusion equation depending on the magnitude of diffusivity, while the computational cost is much lower than other high‐order numerical methods. Copyright © 2011 John Wiley & Sons, Ltd.     

Finite‐difference 4th‐order compact scheme for the direct numerical simulation of instabilities of shear layers

An algorithm based on the 4th‐order finite‐difference compact scheme is developed and applied in the direct numerical simulations of instabilities of channel flow. The algorithm is illustrated in the context of stream function formulation that leads to field equation involving 4th‐order spatial derivatives. The finite‐difference discretization in the wall‐normal direction uses five arbitrarily spaced points. The discretization coefficients are determined numerically, providing a large degree of flexibility for grid selection. The Fourier expansions are used in the streamwise direction. A hybrid Runge–Kutta/Crank–Nicholson low‐storage scheme is applied for the time discretization. Accuracy tests demonstrate that the algorithm does deliver the 4th‐order accuracy. The algorithm has been used to simulate the natural instability processes in channel flow as well as processes occurring when the flow is spatially modulated using wall transpiration. Extensions to three‐dimensional situations are suggested. Copyright © 2005 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.     

Two‐dimensional prediction of time dependent,turbulent flow around a square cylinder confined in a channel

This paper presents two‐dimensional and unsteady RANS computations of time dependent, periodic, turbulent flow around a square block. Two turbulence models are used: the Launder–Sharma low‐Reynolds number k–ε model and a non‐linear extension sensitive to the anisotropy of turbulence. The Reynolds number based on the free stream velocity and obstacle side is Re=2.2×10⁴. The present numerical results have been obtained using a finite volume code that solves the governing equations in a vertical plane, located at the lateral mid‐point of the channel. The pressure field is obtained with the SIMPLE algorithm. A bounded version of the third‐order QUICK scheme is used for the convective terms. Comparisons of the numerical results with the experimental data indicate that a preliminary steady solution of the governing equations using the linear k–ε does not lead to correct flow field predictions in the wake region downstream of the square cylinder. Consequently, the time derivatives of dependent variables are included in the transport equations and are discretized using the second‐order Crank–Nicolson scheme. The unsteady computations using the linear and non‐linear k–ε models significantly improve the velocity field predictions. However, the linear k–ε shows a number of predictive deficiencies, even in unsteady flow computations, especially in the prediction of the turbulence field. The introduction of a non‐linear k–ε model brings the two‐dimensional unsteady predictions of the time‐averaged velocity and turbulence fields and also the predicted values of the global parameters such as the Strouhal number and the drag coefficient to close agreement with the data. Copyright © 2009 John Wiley & Sons, Ltd.     

Higher‐order upwind leapfrog methods for multi‐dimensional acoustic equations

A non‐dissipative and very accurate one‐dimensional upwind leapfrog method was successfully extended to higher‐order and multi‐dimensional acoustic equations. The governing equations in characteristic form and staggered grid were utilized to preserve the accuracy. Fourier analysis was performed to find the accurate scheme for acoustics and the resultant two‐dimensional methods were successfully applied to several classical test cases. Copyright © 2004 John Wiley & Sons, Ltd.     

The Lagrange–Galerkin method for the two‐dimensional shallow water equations on adaptive grids

The weak Lagrange–Galerkin finite element method for the two‐dimensional shallow water equations on adaptive unstructured grids is presented. The equations are written in conservation form and the domains are discretized using triangular elements. Lagrangian methods integrate the governing equations along the characteristic curves, thus being well suited for resolving the non‐linearities introduced by the advection operator of the fluid dynamics equations. An additional fortuitous consequence of using Lagrangian methods is that the resulting spatial operator is self‐adjoint, thereby justifying the use of a Galerkin formulation; this formulation has been proven to be optimal for such differential operators. The weak Lagrange–Galerkin method automatically takes into account the dilation of the control volume, thereby resulting in a conservative scheme. The use of linear triangular elements permits the construction of accurate (by virtue of the second‐order spatial and temporal accuracies of the scheme) and efficient (by virtue of the less stringent Courant–Friedrich–Lewy (CFL) condition of Lagrangian methods) schemes on adaptive unstructured triangular grids. Lagrangian methods are natural candidates for use with adaptive unstructured grids because the resolution of the grid can be increased without having to decrease the time step in order to satisfy stability. An advancing front adaptive unstructured triangular mesh generator is presented. The highlight of this algorithm is that the weak Lagrange–Galerkin method is used to project the conservation variables from the old mesh onto the newly adapted mesh. In addition, two new schemes for computing the characteristic curves are presented: a composite mid‐point rule and a general family of Runge–Kutta schemes. Results for the two‐dimensional advection equation with and without time‐dependent velocity fields are illustrated to confirm the accuracy of the particle trajectories. Results for the two‐dimensional shallow water equations on a non‐linear soliton wave are presented to illustrate the power and flexibility of this strategy. Copyright © 2000 John Wiley & Sons, Ltd.     

Simulation of non‐hydrostatic,density‐stratified flow in irregular domains

A numerical model has been developed for simulating density‐stratified flow in domains with irregular but simple topography. The model was designed for simulating strong interactions between internal gravity waves and topography, e.g. exchange flows in contracting channels, tidally or convectively driven flow over two‐dimensional sills or waves propagating onto a shoaling bed. The model is based on the non‐hydrostatic, Boussinesq equations of motion for a continuously stratified fluid in a rotating frame, subject to user‐configurable boundary conditions. An orthogonal boundary fitting co‐ordinate system is used for the numerical computations, which rely on a fourth‐order compact differentiation scheme, a third‐order explicit time stepping and a multi‐grid based pressure projection algorithm. The numerical techniques are described and a suite of validation studies are presented. The validation studies include a pointwise comparison of numerical simulations with both analytical solutions and laboratory measurements of non‐linear solitary wave propagation. Simulation results for flows lacking analytical or laboratory data are analysed a posteriori to demonstrate satisfaction of the potential energy balance. Computational results are compared with two‐layer hydraulic predictions in the case of exchange flow through a contracting channel. Finally, a simulation of circulation driven by spatially non‐uniform surface buoyancy flux in an irregular basin is discussed. Copyright © 2000 John Wiley & Sons, Ltd.     

Two‐dimensional modeling for stability analysis of two‐phase stratified flow

The effect of wavelength and relative velocity on the disturbed interface of two‐phase stratified regime is modeled and discussed. To analyze the stability, a small perturbation is imposed on the interface. Growth or decline of the disturbed wave, relative velocity, and surface tension with respect to time will be discussed numerically. Newly developed scheme applied to a two‐dimensional flow field and the governing Navier–Stokes equations in laminar regime are solved. Finite volume method together with non‐staggered curvilinear grid is a very effective approach to capture interface shape with time. Because of the interface shape, for any time advancement, a new grid is performed separately on each stratified field, liquid, and gas regime. The results are compared with the analytical characteristics method and one‐dimensional modeling. This comparison shows that solving the momentum equation including viscosity term leads to physically more realistic results. In addition, the newly developed method is capable of predicting two‐phase stratified flow behavior more precisely than one‐dimensional modeling. It was perceived that the surface tension has an inevitable role in dissipation of interface instability and convergence of the two‐phase flow model. Copyright © 2009 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 comparative study of high‐order variable‐property segregated algorithms for unsteady low Mach number flows

In this paper, five different algorithms are presented for the simulation of low Mach flows with large temperature variations, based on second‐order central‐difference or fourth‐order compact spatial discretization and a pressure projection‐type method. A semi‐implicit three‐step Runge–Kutta/Crank–Nicolson or second‐order iterative scheme is used for time integration. The different algorithms solve the coupled set of governing scalar equations in a decoupled segregate manner. In the first algorithm, a temperature equation is solved and density is calculated from the equation of state, while the second algorithm advances the density using the differential form of the equation of state. The third algorithm solves the continuity equation and the fourth algorithm solves both the continuity and enthalpy equation in conservative form. An iterative decoupled algorithm is also proposed, which allows the computation of the fully coupled solution. All five algorithms solve the momentum equation in conservative form and use a constant‐ or variable‐coefficient Poisson equation for the pressure. The efficiency of the fourth‐order compact scheme and the performances of the decoupling algorithms are demonstrated in three flow problems with large temperature variations: non‐Boussinesq natural convection, channel flow instability, flame–vortex interaction. Copyright © 2010 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.     

Accuracy and stability analysis of a second‐order time‐accurate loosely coupled partitioned algorithm for transient conjugate heat transfer problems

In this paper, a second‐order time‐accurate loosely coupled partitioned algorithm is presented for solving transient thermal coupling of solids and fluids, also referred to by conjugate heat transfer. The Crank–Nicolson scheme is used for time integration. The accuracy and stability of the loosely coupled solution algorithm are analyzed analytically. Based on the accuracy analysis, the design order of the time integration scheme is preserved by following a predictor (implicit)–corrector (explicit) approach. Hence, the need to perform an additional implicit solve (a subiteration) at each time step is avoided. The analytical stability analysis shows that by using the Crank–Nicolson scheme for time integration, the partitioned algorithm is unstable for large Fourier numbers, unlike the monolithic approach. Accordingly, using the stability analysis, a stability criterion is obtained for the Crank–Nicolson scheme that imposes restriction on Δt given the material properties and mesh spacings of the coupled domains. As the ratio of the thermal effusivities of the coupled domains reaches unity, the stability of the algorithm reduces. To demonstrate the applicability of the algorithm, a numerical example is considered (an unsteady conjugate natural convection in an enclosure). Copyright © 2013 John Wiley & Sons, Ltd.