An efficient finite difference scheme for free‐surface flows in narrow rivers and estuaries

This paper presents a free‐surface correction (FSC) method for solving laterally averaged, 2‐D momentum and continuity equations. The FSC method is a predictor–corrector scheme, in which an intermediate free surface elevation is first calculated from the vertically integrated continuity equation after an intermediate, longitudinal velocity distribution is determined from the momentum equation. In the finite difference equation for the intermediate velocity, the vertical eddy viscosity term and the bottom‐ and sidewall friction terms are discretized implicitly, while the pressure gradient term, convection terms, and the horizontal eddy viscosity term are discretized explicitly. The intermediate free surface elevation is then adjusted by solving a FSC equation before the intermediate velocity field is corrected. The finite difference scheme is simple and can be easily implemented in existing laterally averaged 2‐D models. It is unconditionally stable with respect to gravitational waves, shear stresses on the bottom and side walls, and the vertical eddy viscosity term. It has been tested and validated with analytical solutions and field data measured in a narrow, riverine estuary in southwest Florida. Model simulations show that this numerical scheme is very efficient and normally can be run with a Courant number larger than 10. It can be used for rivers where the upstream bed elevation is higher than the downstream water surface elevation without any problem. Copyright © 2003 John Wiley & Sons, Ltd.     

Large eddy simulation of high-Reynolds-number atmospheric boundary layer flow with improved near-wall correction

It is highly attractive to develop an efficient and flexible large eddy simulation(LES) technique for high-Reynolds-number atmospheric boundary layer(ABL) simulation using the low-order numerical scheme on a relatively coarse grid, that could reproduce the logarithmic profile of the mean velocity and some key features of large-scale coherent structures in the outer layer. In this study, an improved near-wall correction scheme for the vertical gradient of the resolved streamwise velocity in the s...     

Numerical problems in simulating tidal flows with a frictional-velocity-dependent eddy viscosity and the influence of stratification

The paper deals with the accurate determination of tidal current profiles in both homogeneous and stratified regions when a no-slip condition is used at the seabed with a flow-dependent eddy viscosity related to the depth-mean current or the bed frictional velocity. Calculations show that it is essential to accurately resolve the high-shear region which occurs at the bed and across the pycnocline/thermocline in the case of stratified flow. A computationally accurate and economic method of resolving these regions is demonstrated using the Galerkin method with a set of basis functions designed to accurately reproduce the high-shear layers which occur in these regions. With a flow-independent eddy viscosity a stability analysis can be readily performed and an unconditionally stable algorithm developed. However, with a flow-dependent viscosity, in particular a viscosity computed from the frictional velocity, a non-linear numerical instability can occur. A method of maintaining numerical stability in this case is also described. The importance of near-bed resolution to the computed value of the frictional velocity is demonstrated and its influence on the total tidal velocity profile is illustrated by a number of idealized calculations using various eddy viscosity formulations. The influence of stratification on the computed tidal profiles is shown in the latter part of the paper.     

A solution for the vertical variation of stress,rather than velocity,in a three-dimensional circulation model

A simple technique is presented that allows a numerical solution to be sought for the vertical variation of shear stress as a substitute for the vertical variation of velocity in a three-dimensional hydrodynamic model. In its most general form the direct stress solution (DSS) method depends only upon the validity of an eddy viscosity relation between the shear stress and the vertical gradient of velocity. The rationale for preferring a numerical solution for shear stress to one for velocity is that shear stress tends to vary more slowly over the vertical than velocity, particularly near boundaries. Consequently, a numerical solution can be obtained much more efficiently for shear stress than for velocity. When needed, the velocity profile can be recovered from the stress profile by solving a one-dimensional integral equation over the vertical. For most practical problems this equation can be solved in closed form. Comparisons are presented between the DSS technique, the standard velocity solution technique and analytical solutions for wind-driven circulation in an unstratified, closed, rectangular channel governed by the linear equations of motion. In no case was the computational effort required by the velocity solution competitive with the DSS when a physically realistic boundary layer was included. The DSS technique should be particularly beneficial in numerical models of relatively shallow water bodies in which the bottom and surface boundary layers occupy a significant portion of the water column.     

Effective eddy viscosities in implicit modeling of decaying high Reynolds number turbulence with and without rotation

We assess the applicability of the numerical dissipation as an implicit turbulence model. The nonoscillatory finite volume numerical scheme MPDATA developed for simulations of geophysical flows is employed as an example of a scheme with an implicit turbulence model. A series of low resolution simulations of decaying homogeneous turbulence with and without Coriolis forces in the limit of zero molecular viscosity are performed. To assess the implicit model the long-time evolution of turbulence in the simulations is investigated and the numerical velocity fields are analyzed to determine the effective spectral eddy viscosity that is attributed to the numerical discretization. The detailed qualitative and quantitative comparisons are made between the numerical eddy viscosity and the theoretical results as well as the intrinsic eddy viscosity computed exactly from the velocity fields by introducing an artificial wave number cutoff. We find that the numerical dissipation depends on the time step and exhibits contradictory dependence on rotation: it is overestimated for rapid rotation cases and is underestimated for nonrotating cases. These results indicate that the numerical dissipation may fail to represent the effects of the physical subgrid scale processes unless the parameters of a numerical scheme are carefully chosen.     

Variational parameter estimation for a two-dimensional numerical tidal model

It is shown that the parameters in a two-dimensional (depth-averaged) numerical tidal model can be estimated accurately by assimilation of data from tide gauges. The tidal model considered is a semi-linearized one in which kinematical non-linearities are neglected but non-linear bottom friction is included. The parameters to be estimated (bottom friction coefficient and water depth) are assumed to be position-dependent and are approximated by piecewise linear interpolations between certain nodal values. The numerical scheme consists of a two-level leapfrog method. The adjoint scheme is constructed on the assumption that a certain norm of the difference between computed and observed elevations at the tide gauges should be minimized. It is shown that a satisfactory numerical minimization can be completed using either the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm or Nash's truncated Newton algorithm. On the basis of a number of test problems, it is shown that very effective estimation of the nodal values of the parameters can be achieved provided the number of data stations is sufficiently large in relation to the number of nodes.     

A Hybrid RANS/LES Simulation of Turbulent Channel Flow

Hybrid models combining large eddy simulation (LES) with Reynolds-averaged Navier–Stokes (RANS) simulation are expected to be useful for wall modeling in the LES of high Reynolds number flows. Some hybrid simulations of turbulent channel flow have a common defect; the mean velocity profile has a mismatch between the RANS and LES regions due to a steep velocity gradient at the interface. This mismatch is reproduced and examined using a simple hybrid model; the Smagorinsky model is switched to a RANS model increasing the filter width. It is suggested that a rapid spatial variation in the eddy viscosity is responsible for an underestimate of the grid-scale shear stress and for the steep velocity gradient. To reduce the mean velocity mismatch a new scheme is proposed; additional filtering is introduced to define two kinds of velocity components at the interface between the two regions. The two components are used to remove inconsistency in the velocity equations due to a rapid variation in the filter width. Using the new scheme, simulations of channel flow are carried out with the simple hybrid model. It is shown that the grid-scale shear stress becomes large enough and most of the mean velocity mismatch is removed. Simulations for higher Reynolds numbers are carried out with the k–ε model and the one-equation subgrid-scale model. Although it is necessary to improve the turbulence models and the treatment of the buffer region, the new scheme is shown to be effective for reducing the mismatch and to be useful for developing better hybrid simulations. Received 5 April 2002 and accepted 8 January 2003 Published online 25 March 2003 Communicated by M.Y. Hussaini     

End-to-end differentiable learning of turbulence models from indirect observations

The emerging push of the differentiable programming paradigm in scientific computing is conducive to training deep learning turbulence models using indirect observations. This paper demonstrates the viability of this approach and presents an end-to-end differentiable framework for training deep neural networks to learn eddy viscosity models from indirect observations derived from the velocity and pressure fields. The framework consists of a Reynolds-averaged Navier–Stokes(RANS) solver and a neuralnetwork-represented turbulence model, each accompanied by its derivative computations. For computing the sensitivities of the indirect observations to the Reynolds stress field, we use the continuous adjoint equations for the RANS equations, while the gradient of the neural network is obtained via its built-in automatic differentiation capability. We demonstrate the ability of this approach to learn the true underlying turbulence closure when one exists by training models using synthetic velocity data from linear and nonlinear closures. We also train a linear eddy viscosity model using synthetic velocity measurements from direct numerical simulations of the Navier–Stokes equations for which no true underlying linear closure exists. The trained deep-neural-network turbulence model showed predictive capability on similar flows.     

An implicit three‐dimensional numerical model to simulate transport processes in coastal water bodies

A three‐dimensional baroclinic numerical model has been developed to compute water levels and water particle velocity distributions in coastal waters. The numerical model consists of hydrodynamic, transport and turbulence model components. In the hydrodynamic model component, the Navier–Stokes equations are solved with the hydrostatic pressure distribution assumption and the Boussinesq approximation. The transport model component consists of the pollutant transport model and the water temperature and salinity transport models. In this component, the three‐dimensional convective diffusion equations are solved for each of the three quantities. In the turbulence model, a two‐equation k–ϵ formulation is solved to calculate the kinetic energy of the turbulence and its rate of dissipation, which provides the variable vertical turbulent eddy viscosity. Horizontal eddy viscosities can be simulated by the Smagorinsky algebraic sub grid scale turbulence model. The solution method is a composite finite difference–finite element method. In the horizontal plane, finite difference approximations, and in the vertical plane, finite element shape functions are used. The governing equations are solved implicitly in the Cartesian co‐ordinate system. The horizontal mesh sizes can be variable. To increase the vertical resolution, grid clustering can be applied. In the treatment of coastal land boundaries, the flooding and drying processes can be considered. The developed numerical model predictions are compared with the analytical solutions of the steady wind driven circulatory flow in a closed basin and of the uni‐nodal standing oscillation. Furthermore, model predictions are verified by the experiments performed on the wind driven turbulent flow of an homogeneous fluid and by the hydraulic model studies conducted on the forced flushing of marinas in enclosed seas. Copyright © 2000 John Wiley & Sons, Ltd.     

黄茅海波-流共同作用下的三维悬沙数值模拟

针对黄茅海河口区波流相互作用显著、三维空间结构明显及泥沙运动复杂等特点,联合第3代海浪模式SWAN,建立了波流共同作用下的三维悬沙数学模型.在水流计算中,考虑了波浪产生的辐射应力影响,波浪依赖的表面风应力影响及波浪影响下的底部剪切应力,并采用k-kl紊流闭合模型提供垂向涡粘系数和垂向扩散系数.悬沙扩散方程中的源函数利用切应力法确定,泥沙沉降速度考虑絮凝作用,从而提高黄茅海悬沙场数值模拟精度.通过含沙量实测资料验证,模拟值与实测值符合较好.在模型验证的基础上,讨论了不考虑波浪与考虑波浪两种情况下黄茅海的悬沙分布特征,为工程实际研究提供了依据.     

Parameter estimation for a three‐dimensional numerical barotropic tidal model with adjoint method

The parameters of a three‐dimensional (3‐D) barotropic tidal model are estimated using the adjoint method. The mode splitting technique is employed in both forward and adjoint models. In the external mode, the alternating direction implicit method is used to discretize the two‐dimensional depth‐averaged equations and a semi‐implicit scheme is used for the 3‐D internal mode computations. In this model the bottom friction is expressed in terms of bottom velocity which is different from the previous works. Besides, the bottom friction coefficients (BFCs) are supposed to be spatially varying, i.e. the BFC at some grid points are selected as the independent BFC, while the BFC at the other grid points can be obtained through linear interpolation with these independent BFCs. On the basis of the simulation of M₂ tide in the Bohai and North Yellow Seas (BNYS), twin experiments are carried out to invert the prescribed distributions of model parameters. The parameters inverted are the Fourier coefficients of open boundary conditions (OBCs), the BFC and the vertical eddy viscosity profiles. In these twin experiments, the real topography of BNYS is installed. The ‘observations’ are produced by the tidal model and recorded at the position of TOPEX/Poseidon altimeter data, tidal gauge data and current data. The experiments discuss the influence of initial guesses, model errors and data number. The inversion has obtained satisfactory results and the prescribed distributions have been successfully inverted. The results indicate that the inversion of BFC is more sensitive to data error than that of OBC and the vertical eddy viscosity profiles. Copyright © 2007 John Wiley & Sons, Ltd.     

壁湍流多尺度相干结构复涡黏模型的实验研究

在湍流相干结构动力学方程中,非相干结构成分对相干结构贡献的雷诺应力的模型为涡黏性模型,即涡黏性系数乘以相干结构平均速度变形率的形式.基于非相干结构成分对相干结构贡献的雷诺应力与相干结构速度变形率之间存在相位差的事实,在理论上提出了非相干结构成分对相干结构贡献的雷诺应力复涡黏性模型的假设.应用热线测速技术,在低速风洞中对湍流边界层非相干结构成分对相干结构贡献的雷诺应力与相干结构法向速度变形率之间的相位关系进行了实验测量.通过分析湍流相干结构猝发过程中非相干结构成分对相干结构贡献的雷诺应力与相干结构速度变形率之间的相位关系,研究了相干结构雷诺应力分量与流向速度法向梯度之间的相位差沿湍流边界层法向的变化规律,肯定了湍流相干结构复涡黏性系数模型的合理性.     

A method for estimating wall friction in turbulent wall-bounded flows

We describe a simple method for estimating turbulent boundary layer wall friction using the fit of measured velocity data to a boundary layer model profile that extends the logarithmic profile all the way to the wall. Two models for the boundary layer profile are examined, the power-series interpolation scheme of Spalding and the Musker profile which is based on the eddy viscosity concept. The performance of the method is quantified using recent experimental data in zero pressure gradient flat-plate turbulent boundary layers, and favorable pressure gradient turbulent boundary layers in a pipe, for which independent measurements of wall shear are also available. Between the two model profiles tested, the Musker profile performs much better than the Spalding profile. Results show that the new procedure can provide highly accurate estimates of wall shear with a mean error of about 0.5% in friction velocity, or 1% in shear stress, an accuracy that is comparable to that from independent direct measurements of wall shear stress. An important advantage of the method is its ability to provide accurate estimates of wall shear not only based on many data points in a velocity profile but also very sparse data points in the velocity profile, including only a single data point such as that originating from a near-wall probe.     

Large eddy simulation (2D) using diffusion–velocity method and vortex‐in‐cell

A large eddy Simulation based on the diffusion‐velocity method and the discrete vortex method is presented. The vorticity‐based and eddy viscosity type subgrid scale model simulating the enstrophy transfer between the large and small scale appears as a convective term in the diffusion‐velocity formulation. The methodology has been tested on a spatially growing mixing layer using the two‐dimensional vortex‐in‐cell method and the Smagorinsky subgrid scale model. The effects on the vorticity contours, momemtum thickness, mean streamwise velocity profiles, root‐mean‐square velocity and vorticity fluctuations and negative cross‐stream correlation are discussed. Comparison is made with experiment and numerical work where diffusion is simulated using random walk. Copyright © 2004 John Wiley & Sons, Ltd.     

A semi-implicit lagrangian scheme for 3D shallow water flow with a two-layer turbulence model

A semi-implicit Lagrangian finite difference scheme for 3D shallow water flow has been developed to include an eddy viscosity model for turbulent mixing in the vertical direction. The α-co-ordinate system for the vertical direction has been introduced to give accurate definition of bed and surface boundary conditions. The simple two-layer mixing length model for rough surfaces is used with the standard assumption that the shear stress across the wall region at a given horizontal location is constant. The bed condition is thus defined only by its roughness height (avoiding the need for a friction formula relating to depth-averaged flow, e.g. Chezy, used previously). The method is shown to be efficient and stable with an explicit Lagrangian formulation for convective terms and terms for surface elevation and vertical mixing handled implicitly. The method is applied to current flow around a circular island with gently sloping sides which produce periodic recirculation zones (vortex shedding). Comparisons are made with experimental measurements of velocity using laser Doppler anemometry (time histories at specific points) and surface particle-tracking velocimetry.     

A multi‐level adaptation model of circulation for the western Indian Ocean

A three‐dimensional, fully non‐linear semi‐diagnostic (adaptation) model is described. This model is used to compute the climatological mean circulation and to understand the role of local, steady forcing of the wind and thermohaline forcing on the observed circulation in the western tropical Indian Ocean. The model consists of equations of motion and continuity, sea surface topography, equations of state and temperature, and salinity diffusion equations. While the sea surface topography equation is solved by a successive overrelaxation technique, the other model equations are solved by a leap‐frog numerical scheme. Two versions of the model, having 18 and 33 levels in the vertical direction, were prepared to study climatological mean circulation in the western tropical Indian Ocean. The first numerical experiment is carried out with the 18‐level adaptation model to study the sensitivity of the solution to different values of eddy coefficients. The main scientific rationale behind these numerical experiments was to obtain the most appropriate values of the eddy coefficients for the realistic computation of climatological circulation in the western tropical Indian Ocean. Three numerical experiments were conducted for the month of February to understand the sensitivity of the model solution to different eddy coefficients. The model reproduced the circulation features during February, even with low values of horizontal and vertical eddy coefficients. In the second experiment, the adaptation model, with 33 levels in the vertical direction, is applied to study the seasonal mean climatological circulation at selected depths during Spring in the western tropical Indian Ocean. Adapted (steady state) results of currents, sea surface topography, temperature and salinity anomaly fields are presented. Reasonable agreement is obtained between the model results on currents and the observational data. The computed anomaly fields for temperature and salinity at selected depths during Spring show that the observed temperature and salinity data were adapted with surface wind, flow field and bottom relief of the ocean and that the observed data were found to be fully smoothed during the adaptation stage. Copyright © 1999 John Wiley & Sons, Ltd.     

Three-dimensional hydrodynamics on finite elements. Part II: Non-linear time-stepping model

The development and application of a non-linear 3D hydrodynamic model are described. The model is based on the wave equation rearrangement of the primitive 3D shallow water equations with a general eddy viscosity formulation for the vertical shear. A Galerkin procedure is used to discretize these on simple sixnode elements: linear triangles in the horizontal with linear variations in the vertical. Resolution of surface, bottom and interfacial boundary layers is facilitated and total flexibility is preserved for specifying spatial and temporal variations in the vertical viscosity and density fields. A semi-implicit time-stepping algorithm allows the solutions for elevation and velocity to be uncoupled during each time step. The elevation solution is essentially a 2D wave equation calculation with a stationary sparse matrix representing the gravity waves. With nodal quadrature the subsequent velocity calculation is achieved by factoring only a tridiagonal diffusion matrix representing the vertical viscous terms. As a result the overall calculation scales computationally as only a 2D problem but provides the full 3D solution. Application to field-scale problems is illustrated for the English Channel/Southern Bight system and the Lake Maracaibo system.     

Pure steam condensation model with laminar film in a vertical tube

A new physical model for calculating the liquid film thickness and condensation heat transfer coefficient in a vertical condenser tube is proposed by considering the effects of gravity, liquid viscosity, and vapor flow in the core region of the flow. To estimate the velocity profile in the liquid film, the liquid film was assumed to be in Couette flow forced by the interfacial velocity at the liquid–vapor interface. For simplifying the calculation procedures, the interfacial velocity was estimated by introducing an empirical power-law velocity profile. The resulting film thickness and heat transfer coefficient from the model were compared with the experimental data and the results obtained from the other condensation models. The results demonstrated that the proposed model described the liquid film thinning effect by the vapor shear flow and predicted the condensation heat transfer coefficient from experiments reasonably well.     

A stable algorithm for bed friction in three-dimensional shallow sea modal models

The mathematical formulation of a three-dimensional shallow sea model using a modal expansion in the veitical is briefly described. The importance of the time discretization of the vertical diffusion term and bottom friction term is discussed in some detail. Both stability theory and numerical calculations show the importance of time centring or evaluating the modal form of the viscosity term at the higher time step in order to develop a numerically efficient algorithm. Similar analysis and calculations show that in shallow water it is essential to time centre or evaluate bottom friction at the higher time step. In the case of linear bottom friction it is shown that this condition can be readily accomplished. However, using a quadratic friction formulation (a more physically realistic form), this cannot be readily achieved. A new algorithm is presented whereby a stable solution can be obtained even in shallow water using quadratic bottom friction.