Modeling Transition…

Turbulence prediction has an enormous impact airliners fuel budget, weather information, automotive industry, turbomachinery, biomedical industry and could be easily linked to many other aspects of our daily lives. A lot of focus has been given to the prediction of turbulence through numerical modeling to level of physical fidelity depending on the goal of the practitioner (ranging from calculating lift and drag coefficients of an Aircraft wing or identifying turbulence coherent structures in a boundary layer) and the computational resources at hand. Commercial CFD packages are now offering a diversity of turbulence models for a range of applications and complexity. Yet, a lot less attention has been given for satisfactory simulation of the transition to turbulence onset and process. Moreover, even theoretically speaking, a full understanding of transition onset mechanisms is not available even today.

Transition Mechanisms

Transition from an organized laminar flow to a chaotic, seemingly random turbulent state is strongly dependent on a specific non-dimensional number reflecting on how well momentum is diffused relative to the flow velocity (in the cross-stream direction) and on the thickness of a boundary layer relative to the body – The Reynolds Number. Although the laminar flow solution will remain a solution to the Navier-Stokes equations, above a critical Reynolds number it shall become unstable to small disturbances. Then, a series of events shall occur, some of which are linear and some non-linear that act as instabilities in the transitional process and will eventually lead to a fully turbulent state. The modes of transition onset may be predicted (to some extent) by the linearized stability equations. These equations, derived separately in the beginning of the 20th century by Orr and Sommerfeld (hence termed Orr-Sommerfeld equations) are still investigated by researchers to these day. Orr-Sommerfeld+and+Squire+equations

Orr-Sommerfeld (equation for normal velocity) and Squire equation (for the normal vorticity) for the stability of a parallel undisturbed laminar base-flow

Modal analysis performed on the equations may achieve modal solutions characterized by initial exponential growth, such as Klebanoff (K-type, classical),  Novosibisrsk or Herbert (N-type or H-type, sub-harmonic) all typical to low intensity of incoming turbulence.

Modal transitionTypes of modal transition in a boundary layer (left to right): K-type and H-type (from Berlin et al.)

The initial breakthrough in the field of transition onset and process was the description both theoretically and experimentally of what is termed Tollmien-Schlichting wave instabilities for a low incoming turbulence intensity boundary layer.

TS wavesTollmien-Schlichting wave instability subsequently followed by a series of secondary instability and nonlinear breakdown

A different growth mechanism characterizing a type of transition onset typical to turbomachinery, where high intensity incoming turbulence exist, is the bypass transition (as to bypassing the Tollmien-Schlichting waves instability – my favorite transition mechanism 🙂 ).

Turbulent Boundary Layer (P. Schlatter and D. Henningson of KTH)

A short description of the transition mechanism (my favorite 😉 ) was given in my post: A Forest of Hairpins – on the quest for turbulence coherent structures . Bypass transition.png

Bypass transition mechanism description

Simulating transition

A method based on linear stability theory is the e^n method, originally devised by J.L. van Ingen (TU Delft).  The method is shown to be applicable to boundary layers with pressure gradient, suction and separation. e n methodwhere x0 is the station where the disturbance with frequency ω and amplitude a0 first becomes unstable, n is the “amplification factor” while -αi  is the “amplification rate”. So the “amplification factor” is calculated as a function of x for a range of frequencies giving a set of n-curves. The envelope of these curves gives the maximum amplification factor N which occurs at any x. The method could be used while the stability analysis is based on velocity profiles extracted from a highly resolved boundary layer code. This means that one should construct the infrastructure for such calculation, which is off course not readily available in commercial CFD packages, moreover,  the “amplification factor” needs to be calibrated according to wind tunnel or free-stream flows. e n transition

e^n transition prediction methodology

As transition onset prediction has a direct impact on flow control Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES) are obvious candidates for the job. But the resourcefulness of the plea to a direct numerical description of the equations is a mixed blessing due to the fact that the computational effort in DNS of the Navier-Stokes equations rises as Reynolds number in the power of 9/4 which renders such calculations as prohibitive for most engineering applications of practical interest and it shall remain so for the foreseeable future, its use confined to simple geometries and a limited range of Reynolds numbers  in the aim of supplying significant insight into transition physics that can not be attained in the laboratory.

Osborne Reynolds pipe flow: DNS visualization, from laminar through gradual transition to fully developed turbulence (P.Moin et al., Stanford, 2015)

In LES on the other hand, the large energetic scales are resolved while the effect of the small unresolved scales is modeled using a subgrid-scale (SGS) model and tuned for the generally universal character of these scales. LES has severe limitations in the near wall regions, as the computational effort required to reliably model the innermost portion of the boundary layer (sometimes constituting more than 90% of the mesh) where turbulence length scale becomes very small is far from the resources available to the industry. and although it remains a viable and suitable research tool for transition onset and process description, anecdotally, best estimates speculate that a full LES simulation for a complete airborne vehicle at a reasonably high Reynolds number will not be possible until approximately 2050…

Wall-modeled compressible LES of the Canonical Shock-Turbulent-Boundary-Layer Interaction (P. Moin et al. – Stanford university)

Most of nowadays CFD simulations are conducted with the Reynolds Averaging approach. Reynolds Averaged Navier-Stokes  (RANS) Simulation is based on the Reynolds decomposition according to which a flow variable is decomposed into mean and fluctuating quantities. When the decomposition is applied to Navier-Stokes equation an extra term known as the Reynolds Stress Tensor arises and a modelling methodology is needed to close the equations. The “closure problem” is apparent as higher and higher moments of the set of equations may be taken, more unknown terms arise and the number of equations never suffices. This is of course an obvious consequence to the fact that taking these higher moments is simply a mathematical endeavor and has no physical contribution what so ever.

Reynolds StressReynolds-stress tensor

Levels of modeling are related to the number of differential equations added to Reynolds Averaged Navier-Stokes equations in order to “close” them. 0-equation (algebraic) models are the simplest form of turbulence models, a turbulence length scale is specified in advance through experimenting. 0-equations models are very limited in applications as they fail to take into account history effects, assuming turbulence is dissipated where it’s generated, a direct consequence of their algebraic nature. 1-equation and 2-equations models, incorporate a differential transport equation for the turbulent velocity scale (or the related the turbulent kinetic energy) and in the case of 2-equation models another transport equation for the length scale (or time scale), subsequently invoking the “Boussinesq Hypothesis” relating an eddy-viscosity analog to its kinetic gasses theory derived counterpart (albeit flow dependent and not a flow property) and relating it to the Reynolds stress through the mean strain. In this sense 2-equation models can be viewed as “closed” (as explained in the following post: Understanding The k-ω SST Model) because unlike 0-equation and 1-equation models (with exception of 1-equations transport for the eddy viscosity itself as described in a post: Understanding The Spalart-Allmaras Turbulence Model) these models possess sufficient equations for constructing the eddy viscosity with no direct use for experimental results.

2-equations models do however contain many assumptions along the way for achieving the final form of the transport equations and as such are calibrated to work well only according to well-known features of the applications they are designed to solve. Nonetheless although their inherent limitations, today industry need for rapid answers dictates CFD simulations to be mainly conducted by 2-equations models whose strength has proven itself for wall bounded attached flows at high Reynolds number (thin boundary layers) due to calibration according to the law-of-the-wall.

IMG_0646 The turbulent boundary-layer and the “law of the wall”

y+_Calculation Near wall cell size calculation

There is an inherent problem with RANS as far as transition flows are concerned. The averaging conducted in RANS eliminates the effect of the linear mechanism of disturbance growth. While there are low-Reynolds formulations for 2-eq closures, they should not be confused with transitional Reynolds modeling, as many practitioners rely on low Reynolds models to achieve a measure of transition from laminar to turbulence prediction, these low Reynolds corrections are actually devised to handle the near wall viscous/laminar sublayer. There is no well supported reason to believe one specific calibration of the damping functions in standard low-Reynolds model shall apply to both the viscous/laminar sublayer and transition, especially noting that there are many mechanisms for transition onset. A review of ANSYS Fluent different approaches to turbulence modeling could be found in the presentation below (on going work containing specific subject expanding links and workshops).

turbulence-the-gist-3 Turbulence Modeling – by Tomer Avraham 

The Local Correlation-Based Transition (LCTM) Methodology

F. Menter, R. Langtry (Ansys) and S. Volker (GE) have devised a local correlation-based transition model, mainly based on the following favorable features:
  1. Calibrated prediction of transition onset and length.
  2. Simulation of a diversity of transition mechanisms whether they are a consequence of modal or transient growth.
  3. Formulated locally, as quantities such as the integral thickness of the boundary layer are non-local making calculations problematic to carry out in commercial codes which simply do not provide an infrastructure for such calculations.
  4. Do not affect the turbulence model while in a fully turbulent regime.
As former correlation-based models correlated the momentum thickness to local free-stream conditions, in order to achieve a local correlation-based model, a local property which may be calculated for every grid point as the strain-rate Reynolds number takes part in the development of the local correlation-based trnasition formulation (S is the strain rate): LCTM Re strain rate The strain-rate Reynolds number is proportional to the momentum thickness as follows: LCTM Re strain rate This scaling of the strain-rate Reynolds number means that: LCTM Re strain rate Now there is a framework to develop a local correlation-based transition model.

The γ-Rθt Transition Model

Menter et al. LCTM is based on two additional transport equations for the intermittency-γ and the momentum thickness (or the transition Reynolds number)-Reθt. The intermittency is used to trigger transition locally. In the original formulation the intermittency function is coupled with the turbulence kinetic energy equation in the k-ω SST formulation (see my post: Understanding The k-ω SST Model) to impact the production of turbulence kinetic energy downstream of the transition point in the boundary layer. It could be coupled with a different 2-equation turbulence model, but that means that a re-calibration of the model constants shall be mandatory. In addition to the intermittency transport equation, a transport equation for the transition Reynolds number (momentum thickness) is added which captures the non-local effects of changes in turbulence intensity and free-stream velocity outside the boundary layer. This equation also relates empirical correlations to the transition onset in the intermittency equation. The relation allows for the model incorporation in commercial CFD package as “ready to use”, without the need for interaction from the user due to diverse geometry setups. The intermittency equation is as follows (usual transport equation form): lctm-re-strain-rate3.jpg Now enters the most important (and fun?… 🙂 ) part in each and every construction of closure transport equations, the construction by physical reasoning of the terms in the transport equation. For the intermittency equation (as the left hand side of the above equation is the advection of intermittency ) they are identified as follows: LCTM Re strain rate where: LCTM Re strain rate LCTM Re strain rate F-onset-1 includes the critical Reynolds number,  representing the location where turbulence starts to grow, while the transition Reynolds number represents a location upstream to that, where the velocity profile deviates from the undisturbed (nearly) parallel laminar base-flow profile. The connection between them is obtained by empirical correlation, meaning: LCTM Re strain rate 2 The same notion applies for the F-length in the production term of the intermittency transport equation, which controls the length of transition region, i.e.: LCTM Re strain rate 4 Both correlations are determined through a series of numerical experiments conducted on a flat-plate along with free-stream turbulence intensity, where for the first the critical Reynolds numbers varies with turbulence intensity and the transition Reynolds number is measured based on the most upstream location that the skin friction starts to decrease, and for transition length, experiments were reproduced to extract a validated curve fitting relating the transition length to the transition Reynolds number. These values were then used to develop a correlation to match transition length for a full range of transition Reynolds numbers. The destruction term in the intermittency transport equation reads: LCTM Re strain rate 4 The term ensures for intermittency stays low in the laminar region and enables the representation of re-laminarization conditions, through a return of the intermittency to zero when F-onset is no longer achieved. To avoid the destruction term for fully turbulent regime: lctm-re-strain-rate-5.jpg Subsequently to performing the identification of the different terms in the transport equation, it should be remembered that we are still left out with some added constants to be calibrated. In turbulence modeling calibration of the model is at least as important as the derivation of the model itself. Calibration is achieved with the help of experimental and numerical results of the type of flow that should be modeled. The calibration process is also the first step in which the range of validity of the model would be revealed to close inspection and not just postulated from physical reasoning. For the intermittency equation the calibrated closure constants are: LCTM Re strain rate 5
The transitional Reynolds number shall be extracted from the second transport equation: LCTM Re strain rate 3 This equation, as explained in the above paragraph, takes the turbulence intensity-Tu and the streamwise pressure gradient as empirically correlated with the transitional Reynolds number: LCTM Re strain rate 6 The only unknown in this equation is the transitional Reynolds number. This means that essentially what the transport equation does is taking a non-local correlation of the effect of the change in turbulence intensity and free-stream velocity outside the boundary layer and transforms it to a local transported scalar used to calculate transition length and the critical Reynolds number as explained above. lctm-re-strain-rate-31.jpg The source term in this equation shall force the transported transition Reynolds number in the equation (with the overbar) to match the local value of the transitional Reynolds number (without overbar). The time scale appears for dimensional reasons and has only to scale with both the convective and diffusive terms in the transport equation. F-θt is a blending function (which is used by Menter quite often…) that shuts down the source term to allow the transported scalar to be diffused from the free-stream and effect the boundary layer, meaning it is zero in the free-stream and one in the boundary layer and reads as follows: LCTM Re strain rate 7 F-wake ensures the blending function is inactive in wake regions and the boundary condition for the transported scalar is zero flux at a solid boundary. The constants for the transition Reynolds number transport equation are: LCTM Re strain rate 7 As far as the mesh goes, the first grid must be located at y+=1 (wall units). A coarse mesh shall predict transition onset further upstream with increasing y+ (highly undesirable).
Finally, it should be mentioned that the main achievement of the model is the locality feature which is crucial as far as commercial CFD packages are concerned (one of the reasons why the WALE LES model is highly incorporated in commercial CFD packages while dynamic and mixed models are left out). This not to say non-local operation may not be more efficiently incorporated in an inhouse code assuming you have a code developer on your side… lctm-re-strain-rate-73.jpg

Contours of skin friction of a fully turbulent (left) and transitional (right) of a Eurocopter cabin (from F. R. Menter · R. Langtry · S. Völker – 2006)

Join An exclusive CFD, Thermal analysis, design, physics, and much more…

7 thoughts on “Modeling Transition…

  1. Very nice post as always…
    Just a quick remark about placing your first grid point : velocity and pressure are computed in the CENTER of the cells. So first height should be 2 times the calculated y (if prism) 😉

  2. Pingback: Understanding the Boussinesq Hypothesis and the Eddy-Viscosity Concept – TENZOR (ANSYS Channel Partner) Blog – Mechanical Analysis to the Level of ART

  3. Pingback: Some Fundamental Thoughts About Turbulence Modeling… – Tomer's Blog – All About CFD…

  4. sebastian soto

    Greetings to all.
    I comment on the following concern and problem that I currently have. I am working with a pipe of length 20 [m] and radius 1 [m], where a Bingham fluid flows, within the parameters I have is that its inlet velocity is 0.1 [m / s], Re = 20, Pressure output 0 [Pa] and its respective rheology of the Bingham fluid. As an additional fact, for each variation of the length meter, the pressure must vary 2 [Pa]. (Important data), that is to say that in this case the inlet pressure should be 40 [Pa]. But at the time of solving it in CFX the inlet pressure is much higher even having it placed in the setup with 40 [Pa]. How can I solve or change the setup to achieve the desired inlet pressure?

  5. Pingback: All About CFD… – Index - All About CFD...

  6. Pingback: All About Turbulence Modeling - The “All You Can Eat” Menu - All About CFD...

  7. Pingback: Some Fundamental Thoughts About Emergent Theories... - All About CFD...

Leave a Reply