Introduction

Having established the tools for discretization and methods of analysis of the discretized equation, here we explain methods for solving Navier–Stokes equations of steady and unsteady incompressible viscous flows. Versions of Navier–Stokes equation in primitive and derived variables, expressed in inertial and non-inertial frames have been given in Chapter 2. Various versions of pressure Poisson equation are also derived in Chapter 2, which is to be solved for accurate evaluation of loads and detailed pressure distribution.

Navier–Stokes equation is an evolution equation for vorticity, a primary physical quantity of interest for unsteady laminar and turbulent flows. Primarily, vorticity is generated at physical boundaries for wall-bounded flows, as a consequence of no-slip condition. In free shear layers, vorticity is generated by flow instabilities at interfaces in mixing layers and jets; primary instability mechanisms are attributed to Rayleigh–Taylor, Kelvin–Helmholtz, Görtler mechanisms. Hence, VTE is central for analysis and solution of Navier–Stokes equation. Attendant velocity field can be obtained from the solution of SFE for 2D flows. Poisson equations relating velocity and vorticity field can also be solved for 2D and 3D flows in vorticity–velocity formulation. As stated in Chapter 2, vorticity–stream function formulation is preferred over vorticity–velocity formulation, whenever accurate solution is desired.

Introduction

We have discussed in Chapter 11 about classes of high accuracy compact finite difference schemes for spatial discretization, which along with time discretization schemes, can be used for DNS. It was noted using spectral analysis that the performance of any numerical method is best judged by the spatial and temporal resolutions provided by the combined analysis rather than by the formal order of truncation error terms of individually discretized terms. In this respect, a few compact difference methods are preferred over explicit methods, due to computational efficiency and accuracy. Higher efficiency is due to lesser requirements of points due to higher resolution and accuracy is due to implicit satisfaction of physical dispersion relation with larger time step. It was seen that two second order methods proposed [104, 260, 277], provided higher resolution than sixth order methods [2, 360]. In [104, 277], higher resolution schemes were obtained by minimizing L2-norm of solution error, evaluated as difference of numerical solution from exact solution in the Fourier spectral space for solving 1D convection equation.

There is no unique procedure for the Padè schemes, relating derivative(s) and function values at the nodes. While this can provide flexibility in the choice of a method, if sufficient care is not taken, then it may lead to inconsistency. FVM and FEM can be viewed in this respect to provide logical development in obtaining discrete equations based on physical statement of the problem. Both these methods and other integral methods are said to satisfy conservation principles in a weak sense. In FVM, physical conservation processes are fulfilled in a finite control volume. Physical “densities” of conserved properties are related to fluxes entering and leaving the control volume through control surfaces. These fluxes are created and governed by physical convection and diffusion processes, along with pressure gradient and body forces.

Introduction

In Chapter 2, we have derived most general form of conservation laws which are solved in CFD. For continuum flows, these are given by Navier-Stokes equation. A quick look at them will show that these are non-linear partial differential equations. To be more precise, they are linear with respect to the highest derivative terms and such equations are called quasi-linear PDEs. It is possible to classify such equations based on the behavior of their solutions. In this chapter, we classify quasilinear PDEs, so that we can derive specific numerical methods for each equation Type in subsequent chapters. Despite the observation that different numerical methods are chosen for different classes of PDEs, we will also see here and in later chapters, that there is a generality of approach in treating these PDEs, whose solution and error propagate in time in a unified manner. This will be clearly evident even for time-independent problems, which are solved iteratively, as will be demonstrated through an example in Section 3.3.

Classification of Partial Differential Equations

Consider the moving boundary problem as shown below. This class of problems is also known as the propagation problem.

This book aims at covering the foundations of high accuracy computing methods within the framework of Computational Fluid Dynamics (CFD) in an era of rapidly developing and evolving hardware and software.

From the hardware point of view, huge parallel machines with tens of thousands cores are installed at national facilities and research laboratories giving the practioners of scientific computing tools that they could not have dreamt of a decade ago. The advent of Graphical Processing Units (GPUs) also modifies the course of CFD as everyone tries to strain the computational tools to their last bits and extracts the highest speed-up. This is not surprising as one of the unsolved problems in classical physics is the understanding and control of turbulence in nature and technological applications.

From the software viewpoint, the advent of commercial packages including mesh generators, solvers and graphics tools, provide the numericists with appealing users interfaces and deliver numerical results for extremely different and various problems involving complicated geometries, peculiar boundary conditions and complex physics to be captured. This has had a major impact on the CFD community.

A question that is often raised consists in asking “Why should we not use the simplest schemes and run them on millions (billions) of processors?” The problem as we will discover rapidly is that simple schemes are very often too naive and lead to numerical disaster. We cannot assume that our intellectual indolence will be compensated by the computer’s power. At the end of the day, a bad method will produce inconsistent and poor results.

Introduction

Unlike analytic solutions where we have information everywhere in the domain at all times, in computing we can obtain answers only at discrete points in the domain, called the grid points, at discrete times. This discretization of the independent variables is the first step in computing and a typical section of a discrete grid of points in space are shown in the Cartesian two dimensional plane in Fig. 5.1. Thus, the solution space is discretized preferably with uniformly spaced nodes or discrete points, as shown in the figure.

The approach in obtaining the numerical solution of differential equation rests on representing various differential operators by an equivalent algebraic expressions, i.e., the differential equations are reduced to a set of algebraic relations, which is eventually solved to obtain the unknowns at the grid points. This step of generating algebraic equation from governing equation is known as the discretization process. There is no unique way of doing this. Hence, there is no unique way of solving governing equations. However, there are guiding principles which allow practitioners to choose methods to obtain correct answer. The main guiding principle is the consistency – by which one means that the discretized equations reproduce the original governing equations in the limit of vanishing distance in the grid network shown in Fig. 5.1 and the time step reduced to zero.

Overview on Scientific Computing

The quotation above encapsulates the role scientific computing is expected to play now. Introduction of calculus in seventeenth century provided revolutionary changes in understanding the physical world by analysis and quantification in various branches of science. Prehistory of computing begins with activities in two fronts.

Charles Babbage (1791–1871) laid out detailed plans of mechanical calculating machines: (i) Difference Engines (1821) to construct mathematical tables (without incurring human error!) and (ii) Analytical Engines for the more ambitious task of performing any mathematical task with the help of programmed punched cards. Scaled version of Difference Engine was made during his lifetime. Analytical Engines remained at the conceptual level, though many of its features have found place in modern computers.

On the other hand, progresses were made in developing numerical methods. Basic advances in numerical techniques made in previous centuries encouraged Richardson [225] to seek a solution of system of equations for weather forecasting using a desk calculator. Results were not particularly encouraging from this pioneering effort, but the seeds were sown for what was to become a major initiative!

Introduction

In Chapter 3 we have noted that hyperbolic PDEs are characterized by the presence of real lines along which the highest order derivatives appearing in the differential equations are discontinuous; these lines are the characteristics. Along the characteristics, information propagates for hyperbolic equations. Thus, accurate solution of hyperbolic PDEs is related to the central problem of correctly allowing the signal/disturbance to propagate, as given by the governing PDE. For other classes of differential equations (parabolic or elliptic PDEs) describing equilibrium state of dynamical systems can also support disturbances which can be posed as wave propagation problems. Such wave propagation problems are decided by the dispersion relation relating spatial and temporal scales (as described in Chapter 4 for some mechanical systems) and resultant waves are termed as dispersive waves.

Thus, we discuss about solution methods for hyperbolic PDEs, while focussing on propagating signals and numerical errors. Although, we discuss only about a few handful classical methods of solving hyperbolic PDEs, our attention is equally on the most important problem of signal and error dynamics from a numerical perspective. We would, however, emphasize on signal and error propagation in this chapter, noting that the celebrated error analysis due to von Neumann [41, 53] has been corrected in [259] with respect to a specific example drawn from hyperbolic equation. This forms the second part of this chapter, where the analysis is for identifying mechanisms and sources of numerical error which is independent of which method is being employed.

Introduction

To solve governing equation of motions, we need to resolve all the excited length and time scales. However, even to solve flow past a flat plate, if one takes uniform grid to accommodate the smallest energy carrying length scales, one would be forced to take too many grid points. Numerical solution of the Euler and Navier–Stokes equations for external flow problems requires an outer boundary far away, where some asymptotic boundary conditions apply. This also adds larger requirement on the grid points. Such large problems are poorly convergent and hence, not within the range of the available resources of prevalent high performance computing platforms. This is circumvented by using non-uniform grids. For example, we may decide to take finer grids inside the shear layer (or within the inner layer of a turbulent boundary layer), and take fewer points in the inviscid part of the flow. With nonuniform grids in the physical plane, we usually transform the governing equation in the computational plane, where the spacing is uniform. This also makes writing a code much easier, because of the uniform grid in transformed plane. One needs to write a single code for the transformed plane problem, and for flow past different geometries, one uses the same code provided one can generate a grid mapping from the physical to the transformed plane. Thus, the main solver is grid-independent, and all one needs to do is to generate an appropriate grid transformation either analytically or numerically.

The calculation of rotorcraft performance is largely a matter of determining the power required and power available over a range of flight conditions. The power information can then be translated into quantities such as payload, range, ceiling, speed, and climb rate, which define the operational capabilities of the aircraft. The rotor power required is divided into four parts: the induced power, required to produce the rotor thrust; the profile power, required to turn the rotor through the air; the parasite power, required to move the aircraft through the air; and the climb power, required to change the gravitational potential energy. The aircraft has additional contributions to power required, including accessory and transmission losses and perhaps anti-torque power. In hover there is no parasite power, and the induced power is 65% to 75% of the total. As the forward speed increases, the induced power decreases, the profile power increases slightly, and the parasite power increases until it is dominant at high speed. Thus the total power required is high at hover, because of the induced power with a low but reasonable disk loading. At first the total power decreases significantly with increasing speed, as the induced power decreases; then it increases again at high speed, because of the parasite power. Minimum power required occurs roughly in the middle of the helicopter speed range.

The task in rotorcraft performance analysis is the calculation of the rotor forces and power. Procedures to perform these calculations have been developed in the preceding chapters.

Beams and Rotor Blades

An adequate blade structural model is essential for the prediction of rotor loads and stability. Rotor blades almost universally have a high structural fineness ratio and thus are well idealized as beams. The complexities of rotation, and now multiple load paths and composite construction, have required extensive and continuing efforts to develop appropriate beam models for the solution of rotor problems. For exposition of beam theory, particularly relevant to rotor blade analyses, see Hodges (2006) and Bauchau (1985).

A beam is a structure that has small cross-section dimensions relative to an axial line. Based on the slender geometry, beam theory develops a one-dimensional model of the three-dimensional structure. The deflection of the structure is described as functions of the axial coordinate, obtained from ordinary differential equations (in the axial coordinate). The equations depend on cross-section properties, including two-dimensional elastic stiffnesses. The three-dimensional stress field is determined from the deflection variables. Beam theory combines kinematic equations relating strain measures to deflection variables, constitutive equations relating stress resultants to strain measures, and equilibrium equations relating stress resultants to applied loads. When inertial loads are included, the motion is described by partial differential equations, in time and the axial coordinate.

Rotor Configuration

The helicopter rotor type is largely determined by the construction of the blade root and its attachment to the hub. The blade root configuration has a fundamental influence on the blade flap and lag motion and hence on the helicopter handling qualities, vibration, loads, and aeroelastic stability. The basic distinction between rotor types is the presence or absence of flap and lag hinges, and thus whether the blade motion involves rigid-body rotation or bending at the blade root. A simple classification of rotor hubs has the categories articulated, teetering, hingeless, and bearingless, as sketched in Figures 8.1 to 8.4. With real designs (see Figure 1.2) the distinctions are not as clear as in these drawings.

An articulated rotor has its blades attached to the hub with both flap and lag hinges (Figure 8.1). The flap hinge is usually offset from the center of rotation because of mechanical constraints and to improve the helicopter handling qualities. The lag hinge must be offset for the shaft to transmit torque to the rotor. The purpose of the flap and lag hinges is to reduce the root blade loads (since the moments must be zero at the hinge) by allowing blade motion to relieve the bending moments that would otherwise arise at the blade root. With a lag hinge a mechanical lag damper is also needed to avoid a mechanical instability called ground resonance, involving the coupled motion of the rotor lag and hub in-plane displacement.

Hover is the operating state in which the lifting rotor has no velocity relative to the air, either vertical or horizontal. General vertical flight involves axial flow with respect to the rotor. Vertical flight implies axial symmetry of the rotor flow field, so the velocities and loads on the rotor blades are independent of the azimuth position. Axial symmetry greatly simplifies the dynamics and aerodynamics of the helicopter rotor, as is evident when forward flight is considered. The basic analyses of a rotor in axial flow originated in the 19th century with the design of marine propellers and were later applied to airplane propellers. The principal objectives of the analysis of the hovering rotor are to predict the forces generated and power required by the rotating blades and to design the most efficient rotor.

Momentum Theory

Momentum theory applies the basic conservation laws of fluid mechanics (conservation of mass, momentum, and energy) to the rotor and flow as a whole to estimate the rotor performance. The theory is a global analysis, relating the overall flow velocities to the total rotor thrust and power. Momentum theory was developed for marine propellers by W.J.M. Rankine in 1865 and R.E. Froude in 1885, and extended in 1920 by A. Betz to include the rotation of the slipstream; see Glauert (1935) for the history.

Vertical flight of the helicopter rotor at speed V includes the operating states of hover (V = 0), climb (V > 0), and descent (V < 0) and the special case of vertical autorotation (power-off descent). Between the hover and autorotation states, the helicopter is descending at reduced power. Beyond autorotation, the rotor is producing power for the helicopter. The principal subject of this chapter is the induced power of the rotor in vertical flight, including descent. The key physics are associated with the flow states of the rotor in axial flight. Axial flight of a rotor also encompasses the propeller in cruise (V > 0) and static (V = 0) operation, and a horizontal axis wind turbine (V < 0).

Induced Power in Vertical Flight

In Chapter 3, momentum theory was used to estimate the rotor induced power Pi for hover and vertical climb. Momentum theory gives a good power estimate if an empirical factor is included to account for additional induced losses, particularly tip losses and losses due to nonuniform inflow. In the present chapter these results are extended to include vertical descent. Momentum theory is not applicable for a range of descent rates because the assumed wake model is not correct. Indeed, the rotor wake in that range is so complex that no simple model is adequate. In autorotation, the operating state for power-off descent, the rotor is producing thrust with no net power absorption. The energy to produce the thrust (the induced power Pi) and turn the rotor (the profile power Po) comes from the change in gravitational potential energy as the helicopter descends. The range of descent rates where momentum theory is not applicable includes autorotation.