Parallel 3D computation of unsteady flows around circular cylinders

In this article we present parallel 3D finite element computation of unsteady incompressible flows around circular cylinders. We employ stabilized finite element formulations to solve the Navier-Stokes equations on a thinking machine CM-5 supercomputer. The time integration is based on an implicit method, and the coupled, nonlinear equations generated every time step are solved iteratively, with an element-vector based evaluation technique. This strategy enables us to carry out these computations with millions of coupled, nonlinear equations, and thus resolve the flow features in great detail. At Reynolds number 300 and 800, our results indicate strong 3D features arising from the instability of the columnar vortices forming the Karman street. At Re = 10 000 we employ a large eddy simulation (LES) turbulence model.     

Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method

The development and validation of a parallel unstructured tetrahedral non-nested multigrid (MG) method for simulation of unsteady 3D incompressible viscous flow is presented. The Navier-Stokes solver is based on the artificial compressibility method (ACM) and a higher-order characteristics-based finite-volume scheme on unstructured MG. Unsteady flow is calculated with an implicit dual time stepping scheme. The parallelization of the solver is achieved by a MG domain decomposition approach (MG-DD), using the Single Program Multiple Data (SPMD) programming paradigm. The Message-Passing Interface (MPI) Library is used for communication of data and loop arrays are decomposed using the OpenMP standard. The parallel codes using single grid and MG are used to simulate steady and unsteady incompressible viscous flows for a 3D lid-driven cavity flow for validation and performance evaluation purposes. The speedups and efficiencies obtained by both the parallel single grid and MG solvers are reasonably good for all test cases, using up to 32 processors on the SGI Origin 3400. The parallel results obtained agree well with those of serial solvers and with numerical solutions obtained by other researchers, as well as experimental measurements.     

Numerical simulation of flow past row of square cylinders for various separation ratios

A systematic study for the flow around a row of five square cylinders placed in a side-by-side arrangement and normal to the oncoming flow at a Reynolds number of 150 is carried out through the numerical solution of the two-dimensional unsteady incompressible Navier-Stokes equations. Special attention is paid to investigate the effect of the spacing between the five cylinders on the wake structure and vortex shedding mechanism. The simulations are performed for the separation ratios (spacing to size ratio) of 1.2, 2, 3 and 4. Depending on the separation ratio the following flow patterns are observed: a flip-flopping pattern, in-phase and anti-phase synchronized pattern and non-synchronized pattern. These flow patterns are supposed to be a consequence of the interaction between two types of frequencies viz. the vortex shedding (primary) and the cylinder interaction (secondary) frequencies. At small separation ratio the flow is predominantly characterized by the jet in the gaps between successive cylinders and the secondary frequencies play a role in the resulting chaotic flow. On the contrary, at higher separation ratio the secondary frequencies almost disappear and the resulting flow becomes more synchronized dominated by the primary frequency.     

Block-implicit multigrid calculation of two-dimensional recirculating flows

The performance of a recently developed calculation procedure for steady multidimensional fluid flows is assessed in four laminar two-dimensional recirculating flows representative of those in industrial equipment. The calculation procedure is based on a coupled solution of the momentum and continuity equations by the multigrid technique. The rate of convergence of the algorithm is critically assessed by varying the flow Reynolds number, the number of finite difference nodes, and the numerical underrelaxation factor. The procedure is observed to converge rapidly to an acceptable accuracy in all the flow situations.     

The computation of unsteady non-equilibrium shock tube flows comprising a tailor-made discretization of boundary layer terms

For one-dimensional unsteady non-equilibrium shock tube flows a method of characteristics is developed that treats the boundary layer terms in a specific way. It is shown that a discretization of a boundary layer term tailored to its inverse square root dependence on the distance from the shock front is necessary in order to avoid extremely small step widths. Ionization relaxation in shock tubes is referred to as an example.     

Parallel-multigrid computation of unsteady incompressible viscous flows using a matrix-free implicit method and high-resolution characteristics-based scheme

A three-dimensional parallel unstructured non-nested multigrid solver for solutions of unsteady incompressible viscous flow is developed and validated. The finite-volume Navier–Stokes solver is based on the artificial compressibility approach with a high-resolution method of characteristics-based scheme for handling convection terms. The unsteady flow is calculated with a matrix-free implicit dual time stepping scheme. The parallelization of the multigrid solver is achieved by multigrid domain decomposition approach (MG-DD), using single program multiple data (SPMD) and multiple instruction multiple data (MIMD) programming paradigm. There are two parallelization strategies proposed in this work, first strategy is a one-level parallelization strategy using geometric domain decomposition technique alone, second strategy is a two-level parallelization strategy that consists of a hybrid of both geometric domain decomposition and data decomposition techniques. Message-passing interface (MPI) and OpenMP standard are used to communicate data between processors and decompose loop iterations arrays, respectively. The parallel-multigrid code is used to simulate both steady and unsteady incompressible viscous flows over a circular cylinder and a lid-driven cavity flow. A maximum speedup of 22.5 could be achieved on 32 processors, for instance, the lid-driven cavity flow of Re = 1000. The results obtained agree well with numerical solutions obtained by other researchers as well as experimental measurements. A detailed study of the time step size and number of pseudo-sub-iterations per time step required for simulating unsteady flow are presented in this paper.     

Parallel unstructured multigrid simulation of 3D unsteady flows and fluid-structure interaction in mechanical heart valve using immersed membrane method

In this work, a second-order accurate immersed membrane method (IMM) is adopted to simulate the fluid-structure interaction phenomena in the mechanical heart valves (MHVs). The leaflets of the MHV are immersed in the fluid flows and move on top of the fixed fluid mesh. The blood flow is computed by a 3D parallel unstructured multigrid implicit finite-volume Navier-Stokes solver for incompressible flows. The opening and closing phases of a St. Jude 29 mm MHV are computed under pulsatile inflow to investigate the blood-leaflet interactions. The results show that the moment generated by the fluid pressure is the major cause for the valve motions, while the moment due to the fluid shear stresses is almost negligible. It is also observed that near the end of the opening phase the valve opening speed decelerates, so the valve leaflets have a cushioning effect and avoid a sudden impact on the hinges. For closing phase, jet flows are formed in the central channel and squeeze flows occur in the side channels near the fully closed positions.     

Parallel multigrid finite volume computation of three-dimensional thermal convection

A parallel implementation of the finite volume method for three-dimensional, time-dependent, thermal convective flows is presented. The algebraic equations resulting from the finite volume discretization, including a pressure equation which consumes most of the computation time, are solved by a parallel multigrid method. A flexible parallel code has been implemented on the Intel Paragon, the Cray T3D, and the IBM SP2 by using domain decomposition techniques and the MPI communication software. The code can use 1D, 2D, or 3D partitions as required by different geometries, and is easily ported to other parallel systems. Numerical solutions for air (Prandtl number Pr = 0.733) with various Rayleigh numbers up to 10⁷ are discussed.     

Exact solutions for the unsteady rotating flows of a generalized Maxwell fluid with oscillating pressure gradient between coaxial cylinders

This paper deals with the unsteady rotating flow of a generalized Maxwell fluid with fractional derivative model between two infinite straight circular cylinders, where the flow is due to an infinite straight circular cylinder rotating and oscillating pressure gradient. The velocity field and the adequate shear stress are determined by means of the combine of the sequential fractional derivatives Laplace transform and finite Hankel transform. The exact solutions are presented by integral and series form in terms of the generalized G and Mittag-Leffler functions. The similar solutions can be easily obtained for ordinary Maxwell and Newtonian fluids as limiting cases. Finally, the influence of the relaxation time and the fractional parameter on the fluid dynamic characteristics, as well as a comparison between models, is shown by graphical illustrations.     

Streamwise computation of duct flows

A computational method is presented to calculate momentum and energy transport in two-dimensional viscous compressible duct flows. The flow in the duct is partitioned into finite streams. The difference equations are then obtained by applying momentum and energy conservation principles directly to the individual streams. The method is applicable to laminar and turbulent flows.     

Numerical simulation of low-Reynolds number flows past rectangular cylinders based on adaptive finite element and finite volume methods

The flow past rectangular cylinders has been investigated by two different numerical techniques, an adaptive finite-element (AFEM) and a finite-volume method (FVM). A square and a rectangular cylinder with width-to-height equal to 5 are taken into account. 2D computations have been performed for different Reynolds numbers in order to consider different flow regimes, i.e. the stationary, the periodic and the turbulent flow. The comparison between the two methods regarded both the reliability of the computed solutions and the overall resulting efficiency of the methods. Velocity profiles and integral parameters such as Strouhal number, drag coefficient and recirculation length have been compared. A good agreement between the adaptive FEM and the FVM computations, as well as with the available literature results, has been found. The computational effort has been evaluated in terms of used degrees of freedom in space and time and human resources employed to reach the mesh and timestep-length independence of the solutions. Relevant outcomes of this work are the cross validation of an adaptive FE method and a popular open source FV code.     

The numerical computation of turbulent flows

The paper reviews the problem of making numerical predictions of turbulent flow. It advocates that computational economy, range of applicability and physical realism are best served at present by turbulence models in which the magnitudes of two turbulence quantities, the turbulence kinetic energy $\text{k}$ and its dissipation rate ?, are calculated from transport equations solved simultaneously with those governing the mean flow behaviour. The width of applicability of the model is demonstrated by reference to numerical computations of nine substantially different kinds of turbulent flow.     

Streamwise computation of three-dimensional parabolic flows

A new approach to calculate three-dimensional parabolic flows is presented. The flow field is computed by calculating velocity along a set of streamlines. The dependent variables commonly used in the computation of three-dimensional flows are the three velocity components. In contrast, the dependent variables in the present approach are the streamwise velocity and the two coordinates, in the cross-stream plane, of the chosen streamlines. The streamwise velocity is calculated from the finite difference equations obtained by applying Euler's momentum theorem to streamtubes constructed around the chosen streamlines; the streamline coordinates are calculated from the conservation of mass. Results of the calculations, based on the present approach, are compared with the experimental data for flow through rectangular ducts; the agreement is satisfactory.     

Cores of swirling particle motion in unsteady flows

In nature and in flow experiments particles form patterns of swirling motion in certain locations. Existing approaches identify these structures by considering the behavior of stream lines. However, in unsteady flows particle motion is described by path lines which generally gives different swirling patterns than stream lines. We introduce a novel mathematical characterization of swirling motion cores in unsteady flows by generalizing the approach of Sujudi/Haimes to path lines. The cores of swirling particle motion are lines sweeping over time, i.e., surfaces in the space-time domain. They occur at locations where three derived 4D vectors become coplanar. To extract them, we show how to re-formulate the problem using the Parallel Vectors operator. We apply our method to a number of unsteady flow fields.     

Parallel finite-element computation of 3D flows

The authors describe their work on the massively parallel finite-element computation of compressible and incompressible flows with the CM-200 and CM-5 Connection Machines. Their computations are based on implicit methods, and their parallel implementations are based on the assumption that the mesh is unstructured. Computations for flow problems involving moving boundaries and interfaces are achieved by using the deformable-spatial-domain/stabilized-space-time method. Using special mesh update schemes, the frequency of remeshing is minimized to reduce the projection errors involved and also to make parallelizing the computations easier. This method and its implementation on massively parallel supercomputers provide a capability for solving a large class of practical problems involving free surfaces, two-liquid interfaces, and fluid-structure interactions     

Optimal aerodynamic design of airfoils in unsteady viscous flows

A continuous adjoint formulation is used to determine optimal airfoil shapes in unsteady viscous flows at Re = 1 × 10⁴. The Reynolds number is based on the free-stream speed and the chord length of the airfoil. A finite element method based on streamline-upwind Petrov/Galerkin (SUPG) and pressure-stabilized Petrov/Galerkin (PSPG) stabilizations is used to solve both the flow and adjoint equations. The airfoil is parametrized via a Non-Uniform Rational B-Splines (NURBS) curve. Three different objective functions are used to obtain optimal shapes: maximize lift, minimize drag and minimize ratio of drag to lift. The objective functions are formulated on the basis of time-averaged aerodynamic coefficients. The three objective functions result in diverse airfoil geometries. The resulting airfoils are thin, with the largest thickness to chord ratio being only 5.4%. The shapes obtained are further investigated for their aerodynamic performance. Maximization of time-averaged lift leads to an airfoil that produces more than six times more lift compared to the NACA 0012 airfoil. The excess lift is a consequence of the large peak and extended region of high suction on the upper surface and high pressure on the lower surface. Minimization of drag results in an airfoil with a sharp leading edge. The flow remains attached for close to 70% of the chord length. Minimization of the ratio of drag to lift results in an airfoil with a shallow dimple on the upper surface. It leads to a fairly large value of the time-averaged ratio of lift to drag (～ 17.8). The high value is mostly achieved by a 447% increase in lift and 16% reduction in drag, compared to a NACA 0012 airfoil. Imposition of volume constraint, for the cases studied, is found to result in airfoils that have lower aerodynamic performance.     

On the computation of a matrix inverse square root

The inverse square root of a matrix plays a role in the computation of an optimal symmetric orthogonalization of a set of vectors. We suggest two iterative techniques to compute an inverse square root of a given matrix. The two schemes are analyzed and their numerical stability properties are investigated.     

Numerical computation of transient coaxial entry tube flows

A numerical program was developed to compute transient laminar flows in two dimensions including multicomponent mixing and chemical reaction. The program can compute both incompressible flows and compressible flows at all speeds, and it is applied to describe transient and steady state solutions for low subsonic, coaxial entry, tube flows. Single component, non-reacting flows comprise most of the solutions, but one steady state solution is presented for trace concentration constituents engaging in a second order reaction. Numerical stability was obtained by adding at each calculation point a correction for numerical diffusion errors caused by truncation of the Taylor series used to finite difference the conservation equations. Transient computations were made for fluids initially at rest, then subjected to step velocity inputs that were uniform across each region of the entry plane and were held constant throughout the computation period. For center tube to annulus velocity ratios of 0.5 and 2.0, the bulk fluid in the tube initially moved in plug flow, but strong radial flows developed near the injection plane which moved the fluid into the high shear region between the jets and away from the tube wall. The entrance flows penetrated the bulk flow until steady state was attained. A computation with only the center jet flowing developed a recirculation vortex in the annulus that propagated downstream. The calculation of steady state reacting flows showed formation of a third specie in the mixing zone. Both short and long tube solutions are presented.     

Estimating friction errors in MOC analyses of unsteady pipe flows

Errors arising from the numerical treatment of friction in unsteady flows in small pipe networks are assessed for fixed-grid, method of characteristics analyses (MOC) with no interpolation. Although the results of the study are targeted at general unsteady flows, the underpinning analytical development is based on the behaviour of standing waves. This enables quantitative conclusions to be drawn about the sensitivity of MOC solutions to frequency. The development is undertaken first for a single pipe and then for networks, with and without a loop. The performance of MOC solutions is modelled with the aid of polynomial transfer matrices that are analogous to transfer function matrices available directly from partial differential equations of water hammer. It is found that, in a single pipe, the numerical errors tend to increase damping at all frequencies, but that, in networks, either increased or decreased damping may occur. The errors place both higher and lower limits on the frequencies that can propagate along a pipe in a numerical analysis. This contrasts with the physical case where only a lower limit exists. The work presented is part of a wider project with the long-term aim of automating the selection of numerical grid sizes in MOC analyses of unsteady flows.