Multivariant finite elements for three-dimensional simulation of viscous incompressible flows

Within multivariant elements, which have restricted degrees of freedom at some nodes, different velocity components have different variations. Shape functions for the multivariant elements Q P_o and R P_o are developed. With such shape functions the value of a velocity component within a multivariant element is shown to depend upon all the independent components of velocity at the nodes of the element. The use of the Q₁ P₀ element to simulate flows with discontinuous boundary conditions generated disturbance throughout the flow domain, giving erroneous pressure and velocity distributions. The Q P_o element restricted the disturbance due to such discontinuities to a small region near the singular points, whereas the P P_o element completely eliminated the fluctuations. Flows with discontinuous boundary conditions were simulated with reasonable accuracy by partially relaxing the no-slip condition on the Q₁ P_o elements near the singular points.     

Transition bifurcation branches in non-linear water waves

We are concerned with the numerical computation of progressive free surface gravity waves on a horizontal bed. They are regarded as families of bifurcation branches (λ,A)_Q of constant discharge Q. Numerically we determine two transition values Q₁ and Q₂ with corresponding transition bifurcation branches that classify waves into three disjoint branch sets B₁, B₂ and B₃. Their members are families of waves (λ,A)_Q satisfying the conditions 0<Q² ?Q, Q <Q² ?Q and Q <Q² <B/27, respectively. The bifurcation patterns are analysed in some detail from the computed bifurcation diagram, which shows that in B₁ bifurcation is to the left and the amplitude A increases as the wavelength λ decreases; in B₂ bifurcation is to the right and turning points are observed nearly at breaking point. In B₃ bifurcation is to the right and A increases monotonically with λ.     

New numerical schemes based on a criterion for constructing essentially stable and accurate numerical schemes for convection-dominated equations

In order to obtain stable and accurate numerical solutions for the convection-dominated steady transport equations, we propose a criterion for constructing numerical schemes for the convection term that the roots of the characteristic equation of the resulting difference equation have poles. By imposing this criterion on the difference coefficients of the convection term, we construct two numerical schemes for the convection-dominated equations. One is based on polynomial differencing and the other on locally exact differencing. The former scheme coincides with the QUICK scheme when the mesh Reynolds number (Rm) is $\mathop \[{\textstyle{{\rm 8} \over {\rm 3}}}\] $, which is the critical value for its stability, while it approaches the second-order upwind scheme as Rm goes to infinity. Hence the former scheme interpolates a stable scheme between the QUICK scheme at Rm = $\mathop \[{\textstyle{{\rm 8} \over {\rm 3}}}\] $ and the second-order upwind scheme at Rm = ∞. Numerical solutions with the present new schemes for the one-dimensional, linear, steady convection-diffusion equations showed good results.     

The behaviour of several stress–velocity–pressure mixed finite elements for Newtonian flows

The Q₂/P₁, P/P₁, P₂/P₀ and Q₁/P₀ velocity–pressure mixed elements are extended to the stress–velocity–pressure formulation, using the same interpolants for stress and velocity, and tested in the 4-to-1 contraction problem for Stokes flow. The comparison shows significant differences among them, which are not present when the velocity–pressure formulation is used. To provide a better understanding of the phenomenon, several variants of the previous elements are introduced, obtained by either changing the pressure space or by enriching the stress space with bubble functions. The formulation exhibits a strong sensitivity to the first alternative, while the second produces only a minor effect. These observations are confirmed by a convergence test effected on a regular problem with the explicit analytical solution. Also, as a result of the whole comparison, the P/P/P₁ element looks promising for three-field calculations.     

Practical evaluation of five partly discontinuous finite element pairs for the non‐conservative shallow water equations

This paper provides a comparison of five finite element pairs for the shallow water equations. We consider continuous, discontinuous and partially discontinuous finite element formulations that are supposed to provide second‐order spatial accuracy. All of them rely on the same weak formulation, using Riemann solver to evaluate interface integrals. We define several asymptotic limit cases of the shallow water equations within their space of parameters. The idea is to develop a comparison of these numerical schemes in several relevant regimes of the subcritical shallow water flow. Finally, a new pair, using non‐conforming linear elements for both velocities and elevation (P?P), is presented, giving optimal rates of convergence in all test cases. P?P₁ and P?P₁ mixed formulations lack convergence for inviscid flows. P?P₂ pair is more expensive but provides accurate results for all benchmarks. P?P provides an efficient option, except for inviscid Coriolis‐dominated flows, where a small lack of convergence is observed. Copyright © 2009 John Wiley & Sons, Ltd.     

Numerical simulations of interfacial instabilities on a rotating miscible droplet in a time‐dependent gap Hele–Shaw cell with significant Coriolis effects

Interfacial instability of a rotating miscible droplet with significant Coriolis force in a Hele–Shaw cell is simulated numerically. The influences of the relevant control parameters are first discussed qualitatively by fingering patterns. More vigorous fingerings are found at higher rotational effects, a lower viscosity contrast and a weaker effective surface tension (Korteweg constant). For a time‐dependent gap Hele–Shaw cell, a higher cell lifting rate makes the rotating droplet bear an inward straining flow, which leads to fingering enhancement. On the contrary, a higher pressing rate provides more stable effects by additional squeezing outward flow. A quantitative analysis between the Coriolis effects and tilting angles of fingers is addressed. For arbitrary combinations of all relevant control parameters, the values of tilting angles follow a nearly linear relationship with the Coriolis effects. We estimate the correlation between the relevant control parameters (dimensionless Coriolis factor Re, viscosity parameter R, cell lifting rate a) and tilting angles (θ) of fingers that can be approximated as for significant Korteweg stresses. Copyright © 2005 John Wiley & Sons, Ltd.     

Symbolic computation of flow in a rotating pipe

A perturbation solution of the fully developed flow through a pipe of circular cross-section, which rotates uniformly around an axis oriented perpendicularly to its own, is considered. The perturbation parameter is given by R = 2Ωa²/ν in terms of the angular velocity Ω, the pipe radius a and the kinematic viscosity ν of the fluid. The two coupled non-linear equations for the axial velocity ω and the streamfunction ? of the transverse (secondary) flow lead to an infinite system of linear equations. This system allows first the computation of a given order ?_n, n ? 1, of the perturbation expansion ? = ∑ Rⁿ?_n in terms of ω_n-1, the (n-1)-th order of the expansion ω = ∑ Rⁿω_n, and of the lower orders ?₁,…,?_{n ? 1}. Then it permits the computation of ω_n from ω₀,…,ω_{n ? 1} and ?₁,…,?;_n. The computation starts from the Hagen–Poiseuille flow ω₀, i.e. the perturbation is around this flow. The computations are performed analytically by computer, with the REDUCE and MAPLE systems. The essential elements for this are the appropriate co-ordinates: in the complex co-ordinates chosen the two-dimensional harmonic (Laplace, Δ) and biharmonic (Δ²) operators are ideally suited for (symbolic) quadratures. Symmetry considerations as well as analysis of the equations for ω_n, ?_n and of the boundary conditions lead to general (polynomial) formulae for these functions, with coeffcients to be determined. Their determination, order by order, implies, in complex co-ordinates, only (symbolic) differentiation and quadratures. The coefficients themselves are polynomials in the Reynolds number c of the (unperturbed) Hagen–Poiseuille flow. They are tabulated in the paper for the orders n ? 6 of the perturbation expansion.     

A simplified v2–f model for near‐wall turbulence

A simplified version of the v²–f model is proposed that accounts for the distinct effects of low‐Reynolds number and near‐wall turbulence. It incorporates modified C_ε(1,2) coefficients to amplify the level of dissipation in non‐equilibrium flow regions, thus reducing the kinetic energy and length scale magnitudes to improve prediction of adverse pressure gradient flows, involving flow separation and reattachment. Unlike the conventional v²–f, it requires one additional equation (i.e. the elliptic equation for the elliptic relaxation parameter f_µ) to be solved in conjunction with the k–ε model. The scaling is evaluated from k in collaboration with an anisotropic coefficient C_v and f_µ. Consequently, the model needs no boundary condition on and avoids free stream sensitivity. The model is validated against a few flow cases, yielding predictions in good agreement with the direct numerical simulation (DNS) and experimental data. Copyright © 2007 John Wiley & Sons, Ltd.     

Laminar separation bubble on an Eppler 61 airfoil

Laminar separation bubble that occurs on the suction side of the Eppler 61 airfoil at Re=46000 is studied. The incompressible flow equations are solved using a stabilized finite element method. No turbulence model is used. The variation of the bubble length and its location, with the angle of attack (α), is studied in detail. An abrupt increase in the lift coefficient is observed at α∼4.5°. It is found to be related to a sudden decrease in the separation bubble length at the trailing edge of the airfoil. Significant differences are observed in the results from the 2D and 3D computations. Stall is observed in 3D simulations, but is found to be absent in 2D. The laminar bubble, which fails to reattach in 3D for α>14°, continues to reattach for α as large as 20° in the 2D computations. Reynolds stress calculations in both 2D and 3D indicate the extent to which the outer flow is affected by the presence of bubble. It is found that the Reynolds stress components ${\over{u{^\prime}}{v{^\prime}}}$ and ${\over{u{^\prime}}{w{^\prime}}}$ are of comparable order of magnitude indicating that spanwise fluctuations are significant. The effect of the time window used to compute the time‐averaged aerodynamic coefficients is studied. The time‐averaged and root mean square (rms) value of the aerodynamic coefficients are calculated for both 2D and 3D computations and compared with the previously published experimental results. The 3D computations show good agreement with the earlier data. The variation of the rms value of the aerodynamic coefficients with angle of attack shows certain peaks. The cause of their appearance is investigated. The effect of Reynolds number is studied. The increase in Re at α=10° is found to reduce the bubble length and cause it to move closer to the leading edge. Copyright © 2009 John Wiley & Sons, Ltd.     

An evolutionary optimization of diffuser shapes based on CFD simulations

An efficient and robust algorithm is presented for the optimum design of plane symmetric diffusers handling incompressible turbulent flow. The indigenously developed algorithm uses the CFD software: Fluent for the hydrodynamic analysis and employs a genetic algorithm (GA) for optimization. For a prescribed inlet velocity and outlet pressure, pressure recovery coefficient C (the objective function) is estimated computationally for various design options. The CFD software and the GA have been combined in a monolithic platform for a fully automated operation using some special control commands. Based on the developed algorithm, an extensive exercise has been made to optimize the diffuser shape. Different methodologies have been adopted to create a large number of design options. Interestingly, not much difference has been noted in the optimum C values obtained through different approaches. However, in all the approaches, a better design has been obtained through a proper selection of the number of design variables. Finally, the effect of diffuser length on the optimum shape has also been studied. Copyright © 2009 John Wiley & Sons, Ltd.     

k‐Version of finite element method in 2‐D polymer flows: Oldroyd‐B constitutive model

In this paper, a new mathematical framework based on h, p, k and variational consistency (VC) of the integral forms is utilized to develop a finite element computational process of two‐dimensional polymer flows utilizing Oldroyd‐B constitutive model. Alternate forms of the choices of dependent variables in the governing differential equations (GDEs) are considered and is concluded that u, v, p, τ choice yielding strong form of the GDEs is meritorious over others. It is shown that: (a) since, the differential operator in the GDEs is non‐linear, Galerkin method and Galerkin method with weak form are variationally inconsistent (VIC). The coefficient matrices in these processes are non‐symmetric and hence may have partial or completely complex basis and thus the resulting computational processes may be spurious. (b) Since the VC of the VIC integral forms cannot be restored through any mathematically justifiable means, the computational processes in these approaches always have possibility of spurious solutions. (c) Least squares process utilizing GDEs in u, v, p, τ (strong form of the GDEs) variables (as well as others) is variationally consistent. The coefficient matrices are always symmetric and positive definite and hence always have a real basis and thus naturally yield computational processes that are free of spurious solutions. (d) The theoretical solution of the GDEs are generally of higher order global differentiability. Numerical simulations of such solutions in which higher order global differentiability characteristics of the theoretical solution are preserved, undoubtedly requires local approximations in higher order scalar product spaces . (e) LSP with local approximations in spaces provide an incomparable mathematical and computational framework in which it is possible to preserve desired characteristics of the theoretical solution in the computational process. Numerical studies are presented for fully developed flow between parallel plates and a lid driven square cavity. M1 fluid is used in all numerical studies. The range of applicability of the Oldroyd‐B model or lack of it is examined for both model problems for increasing De. A mathematical idealization of the corners where stationary wall meets the lid is presented and is shown to simulate the real physics when the local approximations are in higher order spaces and when h_d→0. For both model problems shear rate is examined in the flow domain to establish validity of the Oldroyd‐B constitutive model. Copyright © 2006 John Wiley & Sons, Ltd.     

A numerical study of heat island flows: Stationary solutions

We present two‐dimensional numerical simulations of a natural convection problem in an unbounded domain. The flow circulation is induced by a heat island located on the ground and thermal stratification is applied in the vertical direction. The main effect of this stably stratified environment is to induce the propagation of thermal perturbations in the horizontal direction far from the local thermal source. Numerical stationary solutions at Ra?10⁵ are computed in large elongated computational domains: convergence with respect to the domain sizes is investigated at different resolutions. On fine grids, with mesh size , a thermal sponge layer is added at the vertical boundaries: this local damping technique improves the convergence with respect to the domain length. Boussinesq equations are discretized with a second‐order finite volume scheme on a staggered grid combined with a second‐order projection method for the time integration. Copyright © 2008 John Wiley & Sons, Ltd.     

Comparative study of Euler/Euler and Euler/Lagrange approaches simulating evaporation in a turbulent gas–liquid flow

This study deals with the Reynolds‐averaged Navier–Stokes simulation of evaporation in a turbulent gas–liquid flow in a three‐dimensional duct, focussing on the results obtained by a four‐equation turbulence model within the framework of the Euler/Euler approach for multiphase flow calculations: in addition to the two‐equation k?ε model describing the turbulence of the continuous (C) phase, the computational model employs transport equations for the turbulence kinetic energy of the disperse (D) phase and for the velocity covariance q=〈{u}^D{u}^C〉^D. In the present study, the evaporation model according to Abramzon and Sirignano (Int. J. Heat Mass Transfer 1989; 32 :1605–1618) has been extended by introducing an additional transport equation for a newly defined quantity ā, defined as the phase‐interface surface fraction. This allows the change in the drop diameter to be quantified in terms of a probability density function. The source term in the equation describing the dynamics of the volumetric fraction of the dispersed phase α^D is related to the evaporation time scale τ_Γ. The performance of the new model is evaluated by performing a comparative analysis of the results obtained by simulating a polydispersed spray in a three‐dimensional duct configuration with the results of the Euler/Lagrange calculations performed in parallel. Prior to these calculations, some selected (solid) particle‐laden flow configurations were computationally examined with respect to the validation of the background, four‐equation, eddy‐viscosity‐based turbulence model. Copyright © 2008 John Wiley & Sons, Ltd.     

Reynolds‐stress modelling of M = 2.25 shock‐wave/turbulent boundary‐layer interaction

M = 2.25 shock‐wave/turbulent‐boundary‐layer interactions over a compression ramp for several angles (8, 13 and 18°) at Reynolds‐number Re=7 × 10³ were simulated with three low‐Reynolds second‐moment closures and a linear low‐Reynolds standard k–ε model. A detailed assessment of the turbulence closures by comparison with both mean‐flow and turbulent experimental quantities is presented. The Reynolds‐stress model which is wall‐topology free and which uses an optimized redistribution closure, is in good agreement with experimental data both for wall‐pressure and mean‐velocity profiles. Detailed analysis of three components of the Reynolds‐stress tensor (comparison with measurements and transport‐equation budgets) provides a critical evaluation of full Reynolds‐stress models for the separated supersonic compression ramp. The discrepancy observed in the shock‐wave foot region, between computations and measurements for the Reynolds‐stresses profiles, could be explained by considering the experimental shock‐wave oscillation and directions for future modelling work are indicated. Copyright © 2007 John Wiley & Sons, Ltd.     

Near‐wall profiles of mean flow and turbulence quantities predicted by eddy‐viscosity turbulence models

This paper presents for the simple flow over a flat plate the near‐wall profiles of mean flow and turbulence quantities determined with seven eddy‐viscosity turbulence models: the one‐equation turbulence models of Menter and Spalart & Allmaras; the k‐ω two‐equation model proposed by Wilcox and its TNT, BSL and SST variants and the $k-\sqrt{k}L$ two‐equation model. The results are obtained at several Reynolds numbers ranging from 10⁷ to 2.5 × 10⁹. Sets of nine geometrically similar Cartesian grids are adopted to demonstrate that the numerical uncertainty of the finest grid predictions is negligible. The profiles obtained numerically have relevance for the application of so‐called ‘wall function’ boundary conditions. Such wall functions refer to assumptions about the flow in the viscous sublayer and the ‘log law’ region. It turns out that these assumptions are not always satisfied by our results, which are obtained by computing the flow with full near‐wall resolution. In particular, the solution in the ‘log‐law’ region is dependent on the turbulence model and on the Reynolds number, which is a disconcerting result for those who apply wall functions. Copyright © 2009 John Wiley & Sons, Ltd.     

A new hybrid turbulence modelling strategy for industrial CFD

This paper presents a new strategy for turbulence model employment with emphasis on the model's applicability for industrial computational fluid dynamics (CFD). In the hybrid modelling strategy proposed here, the Reynolds stress and mean rate of strain tensors are coupled via Boussinesq's formula as in the standard k–εmodel. However, the turbulent kinetic energy is calculated as the sum of the normal Reynolds‐stress components, representing the solutions of the appropriate transport equations. The equations governing the Reynolds‐stress tensor and dissipation rate have been solved in the framework of a ‘background’ second‐moment closure model. Furthermore, the structure parameter C‐µ has been re‐calculated from a newly proposed functional dependency rather than kept constant. This new definition of C‐µ has been assessed by using direct numerical simulation (DNS) results of several generic flow configurations featuring different phenomena such as separation, reattachment and rotation. Comparisons show a large departure of C‐µ from the commonly used value of 0.09. The model proposed is computationally validated in a number of well‐proven fluid flow benchmarks, e.g. backward‐facing step, 180° turn‐around duct, rotating pipe, impinging jet and three‐dimensional (3D) Ahmed body. The obtained results confirm that the present hybrid model delivers a robust solution procedure while preserving most of the physical advantages of the Reynolds‐stress model over simple k–εmodels. A low Reynolds number version of the hybrid model is also proposed and discussed. Copyright © 2003 John Wiley & Sons, Ltd.     

Three‐dimensional numerical simulation of red blood cell motion in Poiseuille flows

An immersed boundary method based on an FEM has been successfully combined with an elastic spring network model for simulating the dynamical behavior of a red blood cell (RBC) in Poiseuille flows. This elastic spring network preserves the biconcave shape of the RBC in the sense that after the removal of the body force for driving the Poiseuille flow, the RBC with its typical parachute shape in a tube does restore its biconcave resting shape. As a benchmark test, the relationship between the deformation index and the capillary number of the RBCs flowing through a narrow cylindrical tube has been validated. For the migration properties of a single cell in a slit Poiseuille flow, a slipper shape accompanied by a cell membrane tank‐treading motion is obtained for Re , and the cell mass center is away from the center line of the channel due to its asymmetric slipper shape. For the lower Re ?0.0137, an RBC with almost undeformed biconcave shape has a tumbling motion. A transition from tumbling to tank‐treading happens at the Reynolds number between 0.0137 and 0.03. In slit Poiseuille flow, the RBC can also exhibit a rolling motion like a wheel during the migration when the cell is released in the fluid flow with φ = π/2 and θ = π/2 (see Figure 12 for the definition of φ and θ). The lower the Reynolds number, the longer the rolling motion lasts; but the equilibrium shape and position are independent from the cell initial position in the channel. Copyright © 2014 John Wiley & Sons, Ltd.     

Performance of very‐high‐order upwind schemes for DNS of compressible wall‐turbulence

The purpose of the present paper is to evaluate very‐high‐order upwind schemes for the direct numerical simulation (DNS ) of compressible wall‐turbulence. We study upwind‐biased (UW ) and weighted essentially nonoscillatory (WENO ) schemes of increasingly higher order‐of‐accuracy (J. Comp. Phys. 2000; 160 :405–452), extended up to WENO 17 (AIAA Paper 2009‐1612, 2009). Analysis of the advection–diffusion equation, both as Δx→0 (consistency), and for fixed finite cell‐Reynolds‐number Re_Δx (grid‐resolution), indicates that the very‐high‐order upwind schemes have satisfactory resolution in terms of points‐per‐wavelength (PPW ). Computational results for compressible channel flow (Re∈[180, 230]; M?_CL ∈[0.35, 1.5]) are examined to assess the influence of the spatial order of accuracy and the computational grid‐resolution on predicted turbulence statistics, by comparison with existing compressible and incompressible DNS databases. Despite the use of baseline O(Δt²) time‐integration and O(Δx²) discretization of the viscous terms, comparative studies of various orders‐of‐accuracy for the convective terms demonstrate that very‐high‐order upwind schemes can reproduce all the DNS details obtained by pseudospectral schemes, on computational grids of only slightly higher density. Copyright © 2009 John Wiley & Sons, Ltd.     

Impact of mass lumping on gravity and Rossby waves in 2D finite‐element shallow‐water models

The goal of this study is to evaluate the effect of mass lumping on the dispersion properties of four finite‐element velocity/surface‐elevation pairs that are used to approximate the linear shallow‐water equations. For each pair, the dispersion relation, obtained using the mass lumping technique, is computed and analysed for both gravity and Rossby waves. The dispersion relations are compared with those obtained for the consistent schemes (without lumping) and the continuous case. The P₀?P₁, RT₀ and P?P₁ pairs are shown to preserve good dispersive properties when the mass matrix is lumped. Test problems to simulate fast gravity and slow Rossby waves are in good agreement with the analytical results. Copyright © 2008 John Wiley & Sons, Ltd.     

An efficient and stable finite element solver of higher order in space and time for nonstationary incompressible flow

In this paper, we present fully implicit continuous Galerkin–Petrov (cGP) and discontinuous Galerkin (dG) time‐stepping schemes for incompressible flow problems which are, in contrast to standard approaches like for instance the Crank–Nicolson scheme, of higher order in time. In particular, we analyze numerically the higher order dG(1) and cGP(2) methods, which are super convergent of third, resp., fourth order in time, whereas for the space discretization, the well‐known LBB‐stable finite element pair of third‐order accuracy is used. The discretized systems of nonlinear equations are treated by using the Newton method, and the associated linear subproblems are solved by means of a monolithic (geometrical) multigrid method with a blockwise Vanka‐like smoother treating all components simultaneously. We perform nonstationary simulations (in 2D) for two benchmarking configurations to analyze the temporal accuracy and efficiency of the presented time discretization schemes w.r.t. CPU and numerical costs. As a first test problem, we consider a classical ‘flow around cylinder’ benchmark. Here, we concentrate on the nonstationary behavior of the flow patterns with periodic oscillations and examine the ability of the different time discretization schemes to capture the dynamics of the flow. As a second test case, we consider the nonstationary ‘flow through a Venturi pipe’. The objective of this simulation is to control the instantaneous and mean flux through this device. Copyright © 2013 John Wiley & Sons, Ltd.