CHAPTER 1. INTRODUCTION 2
a blade delivering an optimal performance over all the range of wind speeds,
may perform very poorly in unsteady operating conditions such as those as-
sociated with yawed wind: the condition occurring when the wind speed is
no longer orthogonal to the rotor disc. In order to avoid this, an unsteady
flow analysis must be included in the optimization process, to ensure that the
selected design optimum also fulfills the constraint of acceptable performance
under unsteady flow conditions. In summary, the aeromechanical design opti-
mization has to be multi-disciplinary, multi-objective and constrained. The
computational burden for carrying out this design task using high-fidelity
CFD tools throughout the process would be prohibitive. Hence the various
phases of the design optimization would make use of variable-fidelity com-
putational tools, whereby each design functional should be determined with
a level of fidelity sufficient but not above that required for an acceptable
accuracy. For example, the potential fluid model could be used in the steady
fluid/structure iteration for determining the deformed running shape of the
blade, because viscous forces have a small effect on structural deformations.
On the other hand, three-dimensional Navier-Stokes analyses ought to be
used for the aerodynamic performance prediction.
1.2 Computational Fluid Dynamics
The development of the high-speed digital computer during the twentieth
century has had a great impact on the way fluid mechanics is applied to design
problems in modern engineering practice. Analysis and design tasks based
on high-fidelity flow models would have taken years of computational time
with the resources available just a few decades ago. For this reason, these
problems used to be tackled by means of low-fidelity and low-dimensional
models. The rapid growth of computing power experienced in relatively re-
cent years, however, has enabled designers and researchers to afford complex
problems by means of sophisticated flow models and keeping computational
CHAPTER 1. INTRODUCTION 3
times at the level of hours. More specifically, great improvements of physical
accuracy and computational efficiency have been achieved in the area of Com-
putational Fluid Dynamics. In CFD one starts by considering a set of partial
differential equations representing fundamental conservation laws, and then
discretizing such equations by means of accurate numerical schemes. All
these models are subsequently implemented in complex computer programs,
which are used to obtain a discrete representation of the flow field under as-
sessment. In the past, both experimental and theoretical methods have been
used to develop designs for equipment and vehicles involving fluid flow. With
the advent of the digital computer, CFD has emerged as a valuable third op-
tion. Although experimentation continues to be important, the general trend
is towards greater reliance on computer-based predictions in design. Several
high-fidelity flow models are available to CFD users and developers, and all
formulations start from a set of physical conservation laws. One of the most
popular model is that based on the steady and unsteady Reynolds-Averaged
Navier-Stokes (RANS) equations for turbulent flows. In this approach, the
effects of turbulence are taken into account by introducing semi-empirical al-
gebraic or differential equations based on physical evidence. In RANS model,
averaging over the turbulence scales allows one to use relatively coarse com-
putational grids, which are discrete representations of the physical domain.
Avoiding such averaging would force one to use spatial and temporal res-
olutions sufficient to explicitly capture turbulent flow features. This is the
approach of Direct Numerical Solution (DNS), which is accurate but compu-
tationally unaffordable for most problems of engineering interest. The use of
Large Eddy Simulation (LES) lies between that of RANS and DNS. The LES
aims at simulating directly the larger turbulent scales, whereas averaging of
the small flow structures is still performed.
CHAPTER 1. INTRODUCTION 4
1.3 About the existing code
The code is entirely an innovative homemade Fortran 77/90 program [25]
based on the inviscid Euler equations and developed by Dr. M.S. Cam-
pobasso and his CFD team in the Aerospace Department at the University
of Glasgow [9]. It is an explicit code based on the finite volume with a cell
centered approach. The code flow chart is available in App. D. In the fol-
lowing is provided a little summary about the features already present in the
code before beginning this work. Was already performed the adaption for
external flows of the code included the implementation of a far field bound-
ary condition and other several Wall Boundary Conditions (WBCs), the best
of which is the curvature corrected Wall Boundary Condition (WBC) that
will be used to compare the results reached with the present work. The time
marching procedure, with dual time stepping and Runge-Kutta (RK) algo-
rithms, was already analyzed and implemented as well as a time-accurate
preconditioning algorithm. The grid tools for the airfoil and for the cylinder
were correctly developed. In Fig. 1.1(a) and in Fig. 1.1(b) are showed re-
spectively the cylinder and the airfoil mesh schemes: note that the i index is
used to move in ξ direction, whereas the j is used for the η one. For further
details about the previous code features the reader can refer to [19].
(a) (b)
Figure 1.1: Meshes, (a): cylinder, (b): airfoil C-type
CHAPTER 1. INTRODUCTION 5
1.4 Project objectives and overviews
The underlying thread of this project is the computation of the flow field past
a two dimensional airfoil (NACA0012) in steady and unsteady case (harmonic
pitching motion). This work will be carried out using Euler equations and
moving grid with rigid body deformation. The first part of this project will
be devoted to the development of a time-accurate inviscid wall boundary
condition which keeps the spatial order of accuracy imposed by the scheme
used to solve the set of partial differential equations in the interior of the
domain. To calculate the code order of accuracy, the building of a geometry
which has an analytical solution will be necessary: for this purpose a vortex
test case will be employed. We will also develop a new mesh code for the
cylinder in order to run unsteady analyses. The second part will deal with
the implementation of the fluxes generated by the grid motion, in the upwind
spatial discretization. In order to validate the moving mesh procedure, we
will set up a NACA0012 test case in pitching motion and we will compare
the lift to the analytical solution of a flat plate derived by Theodorsen [44].
To obtain a further validation will be also introduced an oscillating cylinder
test case with and without background flow. The third part of the work
will be focused on the converge acceleration techniques for steady case, by
means of the development of two different Implicit Residual Smoothing (IRS)
methods, and for unsteady case using also an implicit treatment of the RK
scheme.
Chapter 2
Governing equations
2.1 Two dimensional Euler equations
2.1.1 Integral ALE formulation
The arbitrary Lagrangian Eulerian (ALE) [20], [17] formulation has emerged
in recent years as a technique that can alleviate many of the drawbacks of
the traditional Lagrangian and Eulerian formulations [110]. Using ALE, the
computational grid need not adhere to the material nor be fixed in space
but can be moved arbitrarily. The grid is continuously moved to optimize
element shapes independently from material deformation. A proper ALE
formulation should reduce to both the Lagrangian and Eulerian formulations
at any degree of freedom as desired. Combining the merits of both the La-
grangian and Eulerian formulations, ALE can easily describe different types
of boundary conditions and prevent mesh distortion. ALE equations are de-
rived by substituting the relationship between the material time derivative
and grid time derivative into the continuum mechanics governing equations.
This substitution gives rise to convective terms in the ALE equations which
account for the transport of material through the grid. ALE is usually termed
a coupled formulation since material deformation and convective effects are
coupled in the same equations. A survey of the ALE literature [9], however,
6
CHAPTER 2. GOVERNING EQUATIONS 7
shows that the majority of ALE analyses, whether quasi-static or dynamic,
are based on the computationally convenient operator split technique. In this
approach, material deformation and convective effects are treated separately.
Thus each time step may be divided into two steps: a regular Lagrangian
step followed by an Eulerian step. The main advantage of this technique over
the fully coupled approach is the reduction in the cost of implementation of
ALE to current Lagrangian codes as the Lagrangian step is unchanged and
only the Eulerian step algorithm needs to be added. Moreover, the decou-
pling of the Lagrangian and Eulerian steps results in simpler equations to be
solved. However, from the theoretical point of view, the fully coupled ALE
approach represents a true kinematic description in which material deforma-
tion is described relative to a moving reference configuration.
We will now develop the integral formulation of the Euler equations with
moving grid; the extension to viscous flows is straightforward and can be
found in [12]. If the grid moves, but the coordinate system remains fixed
and the Cartesian velocity components are used, the only change in the
conservation equations is the appearance of the relative velocity in convective
terms. From [18], we consider the 1D continuity equation
∂ρ
∂t
+ ∂ (ρv)
∂x
= 0 (2.1)
Integrate this equation over a control volume whose boundaries move with
time, we get
∫ x2(t)
x1(t)
∂ρ
∂t
dx+
∫ x2(t)
x1(t)
∂ (ρv)
∂x
dx = 0 (2.2)
The second term of Eq. 2.2 causes no problem whereas the first requires
the use of Leibnitz’s theorem. By doing so, Eq. 2.2 becomes
d
dt
∫ x2(t)
x1(t)
ρdx−
[
ρ2
dx2
dt
− ρ1
dx1
dt
]
+ ρ2v2 − ρ1v1 = 0 (2.3)
The derivative dxdt represents the velocity with which the grid (integration
CHAPTER 2. GOVERNING EQUATIONS 8
boundary) moves, let us call it vb. The terms in square bracket have therefore
a form similar to the last two terms involving fluid velocity, so we can rewrite
the Eq. 2.2 as
d
dt
∫ x2(t)
x1(t)
ρdx+
∫ x2(t)
x1(t)
∂ [ρ (v − vb)]
∂x
dx = 0 (2.4)
When boundary moves with the fluid velocity (i.e. vb = v) the second
integral of Eq. 2.4 equals zero and we get the Lagrangian mass conservation
equation
dm
dt
= 0
The three dimensional version of Eq. 2.3 obtained by using the three dimen-
sional version of Leibnitz’s theorem gives
d
dt
∫
Ω
ρdΩ +
∫
S
ρ (v − vb) · ndS (2.5)
The integral form of the conservation equation for the ith momentum
component takes the following form when the Control Volume (CV) surface
moves with velocity ub
d
dt
∫
Ω
ρuidΩ +
∫
Ω
ρui (v − vb) · ndS =
=
∫
S
(
τi,j~ij − p~ii
)
· ndS
∫
Ω
bidΩ
Conservation equations for scalar quantities are easily derived from the
corresponding equations for a fixed CV by replacing the velocity vector on
the convective term with the relative velocity v − vb. When the boundary
moves with the same velocity as the fluid, the mass flux through the CV
face is zero. If this holds for all CV faces, then the same fluid remains
within the CV and it becomes a control mass; we then have the Lagrangian
description of the fluid motion. On the other hand, if the CV does not move,
the equations are the standard ones. When the location of the grid is known
as function of time, solution of the Navier Stokes equations poses no new
CHAPTER 2. GOVERNING EQUATIONS 9
problems; we simply calculate the convective fluxes (eg. the mass fluxes)
using the relative velocity components at the cell faces. However, when the
cell faces move, conservation of mass, and all other conserved quantities, are
not necessarily ensured if the grid velocities are use to calculate the mass
fluxes. In conclusion the Euler equations for two dimensional moving grid
problems are
∂
∂t
∫
V
UdV +
∮
S
[(F − Uub)nx + (F − Uvb)ny] dS = 0 (2.6)
with
U =
ρ
ρu
ρv
ρE
F =
ρu
ρu2 + p
ρuv
ρuH
G =
ρv
ρuv
ρv2 + p
ρvH
where
• ρ is the density
• (u, v) are the Cartesian components of the velocity
• E is the total energy
• H is the total enthalpy
• p is the pressure
• V is the deforming control volume
• S is the bounding surface of V
• (ub, vb) are the x and y components of the velocities of the control
surface S.
The pressure is related to the conservative flow variables (U) by the equa-
tion of state
CHAPTER 2. GOVERNING EQUATIONS 10
p = (γ − 1)
(
ρE − ρ(u
2 + v2)
2
)
(2.7)
where γ is the ratio of specific heats, generally taken as 1.4.
Now we would develop the geometric conservation law for moving grid. Start-
ing from [30] with the general case, V (t) is neither fixed nor a material vol-
ume. The surfaces of V move, but not with the local fluid velocity. Consider
the Leibnitz theorem we get
d
dt
∫ b(t)
x=a(t)
F (x, t) dx =
∫ b
a
∂F
∂t
dx+ db
dt
F (b, t)− da
dt
F (a, t) (2.8)
The first term on the RHS in Eq. 2.8 is the integral of ∂F∂t , the second
term is due to the gain of F at the outer boundary moving at a rate dbdt , and
the third term is due to the loss of F at the outer boundary moving at dadt .
Generalizing the Leibnitz theorem we get
d
dt
∫
V (t)
F (x, t) dV =
∫
V (t)
∂F
∂t
dV +
∫
A(t)
dA · uAF (2.9)
where uA is the velocity of the boundary. For fixed volume uA = 0, and
the Eq. 2.9 becomes
d
dt
∫
V
F (x, t) dV =
∫
V
∂F
∂t
dV (2.10)
For a material volume V (t) the surfaces move with the fluid, so that uA =
u, with u being the fluid velocity. Using the Reynolds transport theorem the
Eq. 2.9 becomes
D
Dt
∫
V
F (x, t) dV =
∫
V
∂F
∂t
dV +
∫
A
dA · uAF (2.11)
Another form of the transport theorem is derived by using the mass con-
servation. Using the Gauss theorem the Eq. 2.11 becomes
D
Dt
∫
V
FdV =
∫
V
(
∂F
∂t
+ ∂Fuj
∂xj
)
dV (2.12)
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.
