Project topics and materials
Get Complete Project Material File(s) Now! »
Mathematical model
As mentioned in the Introduction, the aim of the present study is to shed light on the physical behavior and stability characteristics of vortex breakdown in an open flow. The current mathematical flow model consists of an approximation of the flow geometry, and governing equations which describe the physics of the flow inside the considered domain.
Flow Configuration
Figure 2.1 shows a schematic diagram of a turbulent jet in a cylindrical domain of outer radius R and axial length x0. Cylindrical coordinates are used, where x, r, and θ denote the axial, radial, and azimuthal directions, respectively.
The flow enters through the inlet (located on the left of the domain in figure 2.1). The characteristic radius of the jet vortex core under investigation is rcore. The characteristic velocity is chosen as the inlet axial velocity, ux0. As the jet develops downstream, ambient fluid is entrained into the flow. The main purpose of the simulation is to leave the bound-aries open to model the behavior of the fluid passing through them. In this case, the flow outside the domain of interest can affect the flow inside. In order for the jet flow to exit the domain, the outlet boundary must also remain open. The open outlet boundary allows for the propagation of disturbances into the inside flow which is essential for the vortex breakdown phenomenon.
Equations of motion
In the previous section details of the flow configuration associated with the boundaries of the cylindrical domain were considered. This section is dedicated to the equations describ-ing the motion of unsteady, axisymmetric, incompressible viscous fluid in the absence of gravity. In order to derive this system of equations, conservation of mass and momentum of fluid must be considered. The use of non-dimensional quantities allows the resulting system to be easily applied to a broad range of similar problems.
Swirling flow is described by several physical and geometrical parameters such as the radius of the vortex core, rcore, the inlet axial ux0, radial ur0 and azimuthal uθ0 velocity components, viscosity and density ρ. These parameters can be combined to form two non-dimensional parameters: the Reynolds number Re and the swirl number S. As mentioned above we consider an incompressible fluid with constant density, ρ = 1. The Reynolds number is here defined as Re = ux0rcore , where ν is the kinematic viscosity of the fluid. The second important non-dimensional parameter is the swirl number, defined as the ratio of the azimuthal velocity at the edge of the core to the axial freestream velocity at the inlet: S = uθ0(rcore) . ux0
We nondimensionalize equations of motions with the following characteristic quantities:
• length [L] = rcore
• density [ρ] = 1
• velocity [U ] = ux0
• time [T ] = rcore .
The definition of the swirl parameter S implies uθ0(r = 1) = S for the azimuthal ve-locity profile at the inlet.
Boundary conditions
At the inlet, the Grabowski profile [37] is used for the radial velocity, where axial and azimuthal velocity components are defined piecewise for the regions inside and outside of a characteristic radius. The Grabowski profile represents a smooth change from solid body rotation inside the characteristic radius to a potential flow farther away, see figure 2.2.
The velocity profile at the inflow boundary is assumed to be axisymmetric and constant over time
ux(0, r) = 1,
ur (0, r) = 0, (2.2)
uθ (0, 0 ≤ r ≤ 1) = Sr(2 − r2),
uθ (0, 1 ≤ r) = S/r.
The inlet condition at x = 0 was chosen in order to approximate the experimentally-measured velocity in vortex cores such as those of trailing vortices, and were used by Mager[68] in his analytical investigation of vortex breakdown.
As previously mentioned, the lateral boundary must remain open to allow entrainment of the external fluid into the domain. To accomplish this goal, a no normal viscous traction boundary condition in the radial direction is used [14] τ • n = 0
where τ represents the viscous stress tensor and n stands for the unit normal vector in the lateral directions. In cylindrical coordinates this equation can be rewritten in component form as
∂ur (x, R) = 0,
∂r
∂ux
∂ur (2.3)
∂x (x, R) + ∂r (x, R) = 0,
∂uθ uθ
(x, R) − (x, R) = 0.
∂r r
In order to allow for mass and momentum to be exchanged across the outlet boundary, the outflow convective boundary conditions used in the numerical computations read
∂ux + C ∂ux = 0, ∂ur + C ∂ur = 0, ∂uθ + C ∂uθ = 0. (2.4)
∂t ∂x ∂t ∂x ∂t ∂x
Following Ruith et al.[79] we choose the convective velocity to be C = 1. In the steady case, equations (2.4) reduce to
∂ux = 0, ∂ur = 0, ∂uθ = 0,
∂x ∂x ∂x
which has the advantage not to depend on the value of C.
Due to the axisymmetry of the flow, we also impose the conditions ur (x, 0) = 0 and uθ (x, 0) = 0 at the centerline.
Direct Numerical Simulation
In order to carry out the research contained in this dissertation, a computer code has been adapted based on the code developed by Nichols [72] (calculation of reacting jets in the low Mach number approximation) to study constant-density viscous axisymmetric rotating flows. The system of partial differential equations is solved by a fourth-order Runge-Kutta scheme in time and a sixth-order discretization scheme in space similar to the one used in Ref.[72, 20]. Due to the complexity of the code, not every detail can be discussed. This section provides an overview of the algorithm used and then focuses on a few numerical issues that arose during the code development.
Mesh
The DNS code solves the Navier-Stokes equations and the continuity equation in cylindri-cal coordinates (2.1) on a staggered mesh, with highly accurate cell-centered sixth-order compact schemes discussed by Lele[62] (with the Fully Included Metrics (FIM), discretiza-tion of a differential equation and evaluation of derivatives directly on a nonuniform grid, which is based on the Taylor series expansion) used to compute spatial derivative in the x and r directions. The benefits of a staggered mesh include a more accurate representation of advection for large scales, but staggering also eliminates numerical modes with upstream propagating group velocity which otherwise would influence the stability properties of the overall flow.
In order to achieve a high resolution in the radial direction an algebraically mapped mesh was used[43]. Following Ref.[96], it is important to control the number of points in a specified subsection of the physical domain in order to achieve high accuracy. In the current investigation the flow dynamics dominates the area close to the centerline, leading to a local flow deceleration and the rise of recirculation region. The dynamics of the flow at the lateral boundary is important as well, since it allows for the entrainment of the external flow into the computational domain. Therefore the following mapping was used to cluster grid points near the lateral boundary r = R and put half the grid points in the interior (most important) region 0 ≤ r ≤ Rint
r = a 1 + rˆ , −1 ≤ rˆ ≤ 1,
b − rˆ
where RintR 2a
a = , b = 1 + .
R − 2Rint R
The grid point separation for this mapping is presented in the figure 2.3 as a function of the radius. A uniform mesh is used in the axial direction.
Boundary conditions
’Boundary conditions are sometimes said to be the ’bane’ of computational fluid dynam-ics’, see Ref. [71]. As was mentioned in the previous section, open boundaries (such as the lateral and outlet boundaries in the current problem) are especially challenging due to the fact that an effective boundary condition must model the motion of the fluid outside the computational domain, for which no information exists. In addition, while cylindrical coor-dinates have many advantages for the simulation of round jets, they introduce an apparent singularity in the governing equation at the centerline r = 0. Even if the centerline is not a true physical boundary, it must be treated carefully in order to avoid numerical difficulties associated with the apparent singularity.
Lateral boundary
As emphasized above, the lateral boundary must remain open to allow entrainment of fluid into the jet. To accomplish this goal, the zero normal viscous traction boundary condition (2.3) is used. To apply this boundary condition we must approximate derivatives in the radial and axial directions on the lateral boundary.
A compact scheme is used to compute the x derivatives, since this technique does not involve a significant amount of extra computations. The computations of r− derivatives via compact schemes, however, would entail the involvement of data from the entire domain; for this reason one-sided differences are employed. Furthermore, the compact schemes them-selves involve one-sided differences near the boundary and thus induce a loss accuracy near the boundary. Therefore, for additional accuracy at the lateral boundary a higher order one-sided finite difference scheme may be selected.
As mentioned before, the convective boundary condition for each velocity flux com-ponent (2.4) is imposed at the outlet. This boundary condition, also known as a non-reflective, radiative, absorbing or advective condition is spatially and temporally local and thus computationally very efficient to implement. In order to do so, a simple one-sided finite difference is used to approximate the x− derivative at the outlet in (2.4). It was found that this naive discretization is sufficient and numerically efficient.
We define r as a unit vector pointing away from the centerline. This vector changes direction across the centerline, which means that certain flow variables change sign to com-pensate and thus are multi-valued at the centerline. This occurs due to the definition of the coordinate system and does not depend on any sort of physical discontinuity. Considering the flow convecting perpendicularly to the centerline, the radial velocity ur is negative as the flow approaches the centerline from one side and then changes to being positive as soon as the flow crosses the centerline and starts moving away from it; from a physical point of view, however, it remains continuous.
A key property in the relation between single-valued and multi-valued variables in cylindrical coordinates is that taking the radial derivative of a single-valued variable will result in a variable that can be multi-valued at the centerline. Therefore, the multi-valued variables must be antisymmetric in the vicinity of the centerline; the converse is also true, meaning that taking the r-derivative of a multi-valued variable will yield a single-valued variable. Multiplying or dividing by r has the same effect, changing between single and multi-valued variables.
Table of contents :
1 Introduction
1.1 Swirling flows: introduction
1.2 Vortex breakdown and its applications
1.2.1 Definition and classification
1.2.2 Vortex breakdown in tornadoes
1.2.3 Vortex breakdown at the tip of delta wings
1.2.4 Combustion and cyclone chambers
1.3 Previous investigations and state of the art
1.3.1 Experimental investigations
1.3.2 Theoretical and numerical studies
1.4 Outline of the thesis
2 Numerical Methodology
2.1 Mathematical model
2.1.1 Flow Configuration
2.1.2 Equations of motion
2.1.3 Boundary conditions
2.2 Direct Numerical Simulation
2.2.1 Mesh
2.2.2 Boundary conditions
2.2.3 Projection method
2.2.4 Poisson equation
2.3 Steady state solutions
2.3.1 Selective Frequency Damping (SFD)
2.3.2 Recursive Projection Method (RPM)
2.3.3 Pseudo-arclength RPM continuation [97, 50]
2.4 Validation
2.4.1 Steady states
2.4.2 Temporal Dynamics
2.4.3 Influence of the Reynolds number Re
3 Viscous Effect
3.1 The Bifurcation Structure of Viscous Steady Axisymmetric Vortex Breakdown with Open Lateral Boundaries
3.1.1 Introduction
3.1.2 Mathematical Model
3.1.3 Numerical Simulations
3.1.4 The critical state of inviscid swirling flow
3.1.5 Asymptotic expansion of near-critical swirling flows in the large Reynolds number limit
3.1.6 Grabowski profile
3.1.7 Discussion and Conclusions
3.2 Convergence of the various branches in the bifurcation diagram
3.2.1 Saddle node
3.2.2 Solution with recirculation region
4 Stability of Swirling Jets
4.1 The effect of an axial pressure gradient
4.1.1 General aspects
4.1.2 Flow configuration and theoretical prediction
4.1.3 Numerical procedure
4.2 Hopf bifurcation
4.3 Spiral vortex breakdown
4.3.1 Previous studies
4.3.2 Governing equations
4.3.3 Arnoldi method used by ARPACK
4.3.4 Results
5 Conclusions and outlook
5.1 Conclusion
5.2 Outlook
6 Appendix: The pseudo-arclength continuation algorithm based on the Recursive Projection Method (RPM)
6.1 Bratu-Gelfand equation
6.2 Matlab script
6.2.1 Time step subroutine
6.2.2 Jacobian matrix-vector product and F subroutines
6.2.3 Basis increase subroutine
6.2.4 Basis decrease subroutine
6.2.5 RPM subroutine
6.2.6 Main routine
Bibliography
