NASA Technical Reports Server (NTRS) 20080014294: Arbitrary Steady-State Solutions with the K-epsilon Model

1.3 MB·English

by  NASA Technical Reports Server (NTRS)

Arbitrary Steady-State Solutions
with the K-Epsilon Model
C. L. Rumsey
NASA LangleyResearch Center, Hampton,Virginia, 23681-2199,USA
B. A. Pettersson Reif
Norwegian Defence Research Establishment, P.O. Box 25, NO-2027,Kjeller, Norway
T. B. Gatski
NASA LangleyResearch Center, Hampton,Virginia, 23681-2199,USA
SubmissiontoAIAAJournal
(cid:0)
(cid:1)
Widely-usedformsofthe - turbulencemodelareshowntoyieldarbitrarysteady-
stateconvergedsolutionsthatarehighlydependent onnumericalconsiderationssuch
asinitialconditionsand solutionprocedure. These solutionscontain pseudo-laminar
regions of varyingsize. By applyinga nullclineanalysisto the equation set, it ispos-
sibletoclearlydemonstratethereasonsfortheanomalousbehavior. Insummary,the
degenerate solutionacts asa stable (cid:2)xed point under certain conditions, causing the
numerical method to converge there. The analysis also suggests a methodology for
preventingtheanomalousbehaviorinsteady-statecomputations.
I. Introduction
The predictionof turbulent (cid:3)ow (cid:2)elds forengineering purposes continues to be dominated
byemployingtheReynolds-averagedNavierStokes(RANS)approach. Althoughtherearesev-
eral classes ofclosure methodologiesavailable, the class of two-equationlinear eddy viscosity
(cid:2) (cid:3)
models may be the most popular. Included in this class is the - model, which has many
variants. Such models are most effectively utilized when the focus is on mean (cid:2)eld dynamics
1of24
ratherthandetailedbehaviorofthestatisticalmoments. Thetwo-equationformulationyieldsan
(cid:4)
eddy viscosity that directlyin(cid:3)uences the mean (cid:3)owbehavior. The equation can be directly
derived from the full Reynolds stress transport equation (by taking the trace). Closure of this
(cid:4)
equation then requires a constitutive equation for the Reynolds stress tensor as well as the
transportandpressure-diffusionterms. FortheReynolds stresstensor, alinearrelationwiththe
mean strain rate tensor is assumed with the proportionalitycoef(cid:2)cent being the turbulenteddy
viscosity. Forthetransportandpressure-diffusiontermsgradienttransportmodelsaregenerally
(cid:5)
assumed. The equationis based ontransportprocesses in thedissipative range, butit is often
viewedas beingan empiricalmodel.1,2
Naturally, implicit in the use of all RANS models is the assumption that for a given set of
initial and boundaryconditions, a unique solution will be obtained. However, it will be shown
(cid:4) (cid:5)
here that for certain formulations of the - model, portions of the (cid:3)ow (cid:2)eld can converge to
a degenerate solution that is (cid:147)pseudo-laminar(cid:148) in nature. The terminology (cid:147)pseudo-laminar(cid:148)
alludes to the fact that the model does not predict the correct laminar limit in terms of the tur-
(cid:6)(cid:7)(cid:4)(cid:9)(cid:8)(cid:10)(cid:5)
bulent tomean (cid:3)owtime scale ( ), butthe resultingeddy viscosity is still nearzero so the
mean (cid:3)ow behaves as a laminar (cid:3)ow. Unfortunatley, this disturbing behavior is further exas-
cerbated by the fact that the location and spatial extent of these pseudo-laminar regions in the
converged solutions can depend on initial conditions and method of solution! This paper both
demonstrates the problem and performs an analysis that explains why the system of equations
behavesinthismanner. Itshouldbenotedthatthisanomalousbehavioroccursonlyindevelop-
ing (cid:3)ows, i.e., boundary layers. Flowcon(cid:2)gurations that utilize streamwise periodic boundary
conditions(such as fullydevelopedchannel (cid:3)ow)donotin generalencounterthisproblem.
Such issues have not been addressed previously, and the purpose ofthis analysis is to criti-
(cid:4) (cid:5)
cally examine the characteristics of the - model froma dynamical systems standpoint. Mo-
(cid:4) (cid:5)
hammadi and Pironneau3 reported extensive mathematical analysis on the - model, but did
not consider the anomalies described here. Dynamical system analyses have previously been
utilized4(cid:150)6 in order to identify the (cid:2)xed points of two-equation and second-moment closures
in homogeneous shear (cid:3)ow, and to calibrate closure models for equilibrium (cid:3)ows. The issues
addressed in the presentstudy, however, requirethat the evolutiontowardthe equilibriumstate
is understood, particularly the solution trajectories from a given initial condition, and not the
(cid:2)xed points per se. Although the theoretical approach is inherently based on the simplifying
assumptionofhomogeneousturbulence,anattemptisalsomadetoaccountforinhomogeneous
effectsso thatthetheoreticalresultscan bemoreeasilyrelatedtothe fullnumericalsolutionof
the modelinrealisticwall-bounded(cid:3)ows.
2of24
As already alluded to, the dynamical systems analysis undertaken here is not only con-
cernedwith(cid:2)xedpointidenti(cid:2)cationbutalsoaboutdependentvariablephase-planetrajectories
and dependence on initial conditions. The analysis is shown to isolate de(cid:2)ciencies in some
(cid:11) (cid:12)
formulations of the - model that lead to initial condition and solution-method dependent
converged solutions.
II. Illustration of the Problem
(cid:11) (cid:12)
To illustrate the anomalies arising in numerical solutions of the - model, two cases are
considered: (1) a (cid:3)at plate boundary layer (cid:3)ow, and (2) (cid:3)ow over an airfoil. The computer
code CFL3D7 was employed.a As will be shown, the capacity to reach arbitrary steady-state
(cid:11) (cid:12)
solutionsisapropertyofthe - equationsthemselves,soanynumericalmethodwillencounter
the problem.
(cid:11) (cid:12)
A basic formofthe - modelcan bewrittena(cid:31)"s!
(cid:13) (cid:23) (cid:28)(cid:30)(cid:29) (cid:31) (cid:23)
(cid:11) (cid:11)
(cid:13)(cid:15)(cid:14) (cid:16) (cid:17)(cid:19)(cid:18)(cid:20)(cid:12)(cid:22)(cid:21) (cid:23)(cid:25)(cid:24)(cid:27)(cid:26) (cid:21) (cid:23)(cid:25)(cid:24)(cid:27)(cid:26)"’
#%$(cid:15)& (cid:31) ! (1a)
(cid:13) (cid:23) (cid:28)(cid:30)(cid:29) (cid:31) (cid:23)
(cid:12) (cid:12) (cid:12)
(cid:13)((cid:14) (cid:16) (cid:17)(cid:19)(cid:18)1032 24(cid:12)657(cid:21) (cid:23)(cid:25)(cid:24) (cid:26) (cid:21) (cid:23) (cid:24) (cid:26)"’98
(cid:11)*),+(cid:7)-/. +(cid:7)- # & (1b)
-
(cid:13);:3(cid:13)((cid:14) (cid:23) : (cid:23) (cid:14) (cid:26) (cid:23) : (cid:23) (cid:24) (cid:26) # $ # 2 :KJML
(cid:16) (cid:21)=< (cid:16)?>(cid:30)@BA(cid:30)A 2C(cid:16)?>3@ED3F (cid:16)?>3@EG (cid:16)IH 2O(cid:18)
where , +(cid:7)-,. , +(cid:7)- , , - +ONP)/+(cid:7)-
5RQ HS(cid:16)TGU@VA%> (cid:16)XGU@EG3Y
+(cid:7)-/. , , and +WN . The set of values used here for the closure constants are
representative of values typically chosen. However, as will be seen, the anomalous behavior
(cid:17)
exists regardless of their exact numerical values. The rate of turbulent energy production is
(cid:26)^(cid:23) : (cid:23) (cid:24)U(cid:26)
(cid:18) Z\[]Z <_[
de(cid:2)ned as , and the componentsof theReynolds stress tensor is modeledusing
the Boussinesq assumption,
(cid:31)"!(cid:7)(cid:29) (cid:23) (cid:23) (cid:26)
(cid:26) (cid:26) Z%[ Z
Z [ Z (cid:16)?‘ (cid:11)ba [ (cid:18) (cid:23) (cid:24) (cid:26) (cid:21) (cid:23) (cid:24) 8
F [ & (2)
withthe eddyviscositygiven by (cid:31) ! 2
(cid:11)
(cid:16) @
+ N (cid:12) (3)
0P2
In this type of formulation, the function is introduced into Eq. (1b) so that as a solid
boundaryis approachedthe destruction-of-dissipationterm does not become ill-behaved.9 It is
aAnothercode,ISAAC,8hasalsobeenusedandyieldssimilarresults.
3of24
also used to calibrate the model in the log-layerregion. For example, two of the most widely-
c3d
used formsfor are:
d
c dWe(cid:19)fhgjilk/m(cid:10)nUo;pqgCirdtsvu w(cid:25)x(cid:7)y
(4)
w dt{}|R~U(cid:127)6(cid:128)
sCu eSz
where , and
| (cid:128)
c d(cid:129)e(cid:19)f(cid:22)g1iq(cid:130)qm(cid:131)nKo gCiq(cid:132)(cid:131)sCur(cid:133) y
(5)
{"~
svu(cid:134)(cid:133)Se(cid:19)(cid:135) zb(cid:136) (cid:136) i(cid:30)k iq(cid:130)
where and is the distance fromthe wall. The constants and are always
i d i (cid:132) fvg(cid:137)i kv(cid:138)
positive and less than or equal to 1.0, and and are positive constants. Note that
c3d(cid:139)(cid:138) f f^g(cid:140)i (cid:130)(cid:139)(cid:138)9c3d(cid:22)(cid:138) f
forEq. (4)and forEq.(5). Forthepresentdemonstrations,onlyEq.(5)
{
i (cid:130) e f i (cid:132) e(cid:142)(cid:141) (cid:141)3(cid:143)
was used, with the commonly used values , , but Eq. (4) exhibits similar
behavior.
In order to rule out the possibility that the observed arbitrary behavior occurred only for
(cid:127)
z
a particular set of boundary conditions, various freestream boundary conditions for and
were employed. The anomalies occurred regardless of these variations, as exempli(cid:2)ed below.
Becausethesolutionswerehighlydependentonnumericalparameters,thegridsizeitselfcould
also have an in(cid:3)uence on the (cid:2)nal solution. Reasonably (cid:2)ne grids were used for both cases,
but a formal grid independence study was not conducted because it has no meaning when the
equations themselves (andnotthenumerics)canyieldarbitraryresults.
A. FlatPlateBoundaryLayer
(cid:144)K(cid:145)(cid:141)P(cid:146)
In the (cid:3)at plate boundary layer (cid:3)ow considered (cid:2)rst, a freestream turbulent intensity of
{l(cid:150) (cid:150)W(cid:152)v(cid:153)
eI(cid:147) (cid:141)Pz(cid:149)(cid:148) f(cid:10)(cid:151) (cid:143)
( )wasimposedeverywhereinthecomputationaldomainusingagridof at
(cid:153)(cid:154)(cid:152)
sCuCe f (cid:144)}(cid:155) z(cid:157)(cid:156) z
aReynoldsnumber(basedonplatelength)of . Thus,theinitialcondition on
{"£⁄d (cid:153) (cid:152)
z(cid:159)(cid:156)(cid:158) e9z(cid:160)(cid:148)(cid:158) e¡z¢(cid:148) (cid:148) e (cid:145)E(cid:144) f (cid:144)6¥6(cid:155)
waseverywherethesameastheboundarycondition( ). The
dissipation rate boundary condition at in(cid:3)ow, which determines the freestream eddy viscosity
(cid:127) {"£ (cid:130) {l£ d (cid:127)
(cid:148)Cƒ (cid:148) e¤§ (cid:145) f"z¢(cid:148) (cid:148) (cid:158)(cid:148) e9§ (cid:145) f'z (cid:148)(cid:158)
and turbulencefreestreamdecayrate, was set at , or . This
{
“‹« “ (cid:148)›e (cid:144)U(cid:145)Vfi6(cid:144)
yielded a freestreameddy viscosity forthis case of at the boundaries. A sample
£ {3fl
(cid:148)
solution showing the (cid:2)nal -velocity contours is plotted in Fig. 1, for the case with initial
(cid:127) {
(cid:158)(cid:156) e9§ (cid:145) f'z(cid:160)(cid:156)(cid:158) (cid:176) ƒ†– (cid:144)K(cid:145) f
condition . Thesolutionexhibitedapseudo-laminarsolutionupstreamof ,
then a turbulentsolution downstream. Next, differentcases were run with fourdifferentinitial
(cid:127)
(cid:158)(cid:156) (cid:144)K(cid:145)‡(cid:144)(cid:30)(cid:144) §Uf"z (cid:156)(cid:158) §Uf'z (cid:156)(cid:158)
conditions in the(cid:2)eldvaryingfrom to .
As Fig. 2 shows, each converged solution yielded a differentapparent (cid:147)transition(cid:148) location
that was located further downstream with increasing initial dissipation rate values. At suf(cid:2)-
4of24
·3(cid:181)¶(cid:160)•*‚U„"”(cid:160)¶(cid:181)
ciently high levels ( ), the (cid:3)ow remained laminar throughout. Initially, this might
seem like consistent behavior: an increased dissipation should reduce the level of turbulent
”
kinetic energy and therefore should shift the position of (cid:147)transition(cid:148) further downstream.
However, this rationaleis (cid:3)awed fortwo reasons; one more serious than the other froma CFD
practitioner’spointofview.
»(cid:137)•I…7‰(cid:22)(cid:190)"…(cid:25)¿
First, in the presence of mean shear (e.g. ), a laminar solution is characterised
»(cid:7)”(cid:9)(cid:190)(cid:131)·(cid:140)(cid:192) „
by very large values of the time-scale ratio, i.e., . Physically this implies that the
turbulenttimescale is muchgreaterthan the mean(cid:3)ow timescale and as such turbulencedoes
»(cid:7)”(cid:9)(cid:190)(cid:131)·
not persist. But inspection of the (cid:147)laminar(cid:148) regions in the current solution shows that
` „ ” ·
is in fact . In other words, the - model is not really predicting laminar (cid:3)ow at all, but
rather a pseudo-laminar behavior that will be shown in section III to be a stable (cid:2)xed point of
the equations.
Second, andevenmoreseriousisthatthe(cid:2)nalconvergedsolutionshouldnotdependonthe
initial condition in a steady-state computation at all! The converged solution should depend
on the boundary conditions only, and in this case the boundary conditions were the same in
all computations. It should be stressed that the solutions presented in Fig. 2 were very well
´(cid:7)ˆ
converged solutions. The -norm of the density residual dropped by more than 6 orders of
„(cid:154)˜(cid:149)„(cid:131)¯U˘%˙(cid:201)¨
magnitude,toapproximately ,asshowninFig.3. Thesolutionsafter25,000multigrid
cycles showed noperceptibledifferencesfromthosesolutions obtainedafter2500cycles.
B. RAE2822Airfoil
As a second example of anomalous behavior, the (cid:3)ow over an RAE 2822 airfoil at freestream
˚(cid:204)¸˝•¤¯U˛—ˇl(cid:209) (cid:210)(cid:211)•(cid:213)(cid:212)U˛(cid:214)ˇ"(cid:212) (cid:215)C(cid:216)C•9(cid:217)K˛‡(cid:212)(cid:218)˜C„(cid:10)¯3(cid:219)
Machnumber ,angle-of-attack ,and wascomputed. Aplot
showingtheairfoilshapeandresultingpressurecontoursfortheseconditionsisgiveninFig.4.
There is a strong shock wave present on the airfoilupper surfacenear 65% chord, whereas the
(cid:3)owonthelowersurfaceremainssubsonic. Inthisexample,theinitialandboundaryconditions
were kept (cid:2)xed but the numerical solution method was changed. The initial conditions and
” ¸(cid:181) •S(cid:212)K˛‡(cid:209)(cid:220)˜(cid:221)„(cid:10)¯ ˘6(cid:222) · (cid:181)¸ •9(cid:223)U˛(cid:214)ˇ"(cid:209)(cid:224)˜(cid:221)„(cid:10)¯ ˘6(cid:222)
far(cid:2)eldboundaryconditionsin thiscase were set to and .
Thiscorrespondedtoaverylowfreestreamturbulentintensityof0.013%andafreestreameddy
Æ (cid:226)R(cid:190)(cid:131)Æ ¸ •9¯U˛E¯3¯(cid:30)ª
viscosity of . These aretypicalvalues used inCFL3D.7
Figure5showstheskin-frictiondistributiononthelowersurfaceoftheairfoilobtainedusing
two different numerical solution strategies for obtaining converged solutions. The converged
results were completely different, with each suggesting a (cid:147)transition(cid:148) location in a different
5of24
place. The (cid:2)rst case was run using multigrid and 3 levels of mesh sequencing, with 2500
iterations on the coarse grid, followedby 2500 iterationson the medium grid, and (cid:2)nally3000
(cid:228)(cid:30)(cid:229)P(cid:230);(cid:231)bŁP(cid:230)
iterationson the(cid:2)nest grid( ). The second case wasrunwith multigridand2 levels of
mesh sequencing, with 5000 iterations on the medium grid followed by 3000 iterations on the
ØWŒ
(cid:2)nest grid. Although not shown, both cases converged very well, with the -norm of density
residualreduced morethan 3ordersofmagnitude.
º (cid:236)
Both these examples show that caution needs to be exercised when using the - model.
A numerically converged solution does not necessarily constitute the intended solution to the
set of governing equations; it may depend on numerical parameters such as initial conditions
and solution procedure. It is also importantto mention here that many CFD practitioners have
º (cid:236)
noticed that the - equations often fail to go fully turbulent, although the cause has never
been identi(cid:2)edbefore. In fact, it is customaryto buildin ad-hoc (cid:2)xes to attempt to ensure that
º (cid:236)
turbulence always develops. Some of these (cid:2)xes include: (1) restarting - solutions from
another turbulentsolution, (2)setting initialconditions to have turbulent-likelevels ratherthan
freestreamlevels,and(3)imposingatemporarysourcetermintheboundarylayer. Allofthese
(cid:2)xes, in general, are workable ways to avoid the problem; but they do not shed any light on
the reasons behind the problemand were not developedbased on any (cid:2)rmrational foundation.
As a consequence, their generality cannot be assured. In the following section an analysis is
conductedthat makesthereasons clear.
III. Dynamical Systems Analysis
A dynamical systems analysis can be used to determine the temporal dynamics associated
with the numerical solution of systems of equations.10 A so-called nullcline analysis will also
º (cid:236)
be used toidentifysome parametricrestrictionson the - equationsEq. (1),inordertoavoid
arbitrarypseudo-laminarconvergedsolutions.
It is possible togain criticalinsightinto thesolution behaviorforinhomogeneousturbulent
(cid:3)owsthroughananalysisofthehomogeneousformofEq.(1). Initshomogeneouslimit,Eq.(1)
can bewrittenas anonlinear,autonomousequationsystem in theform
(cid:237) Œ
ºb(cid:238) ºb(cid:238)
(cid:237)3(cid:239) (cid:236)
(cid:238) (cid:240) æO(cid:242) (cid:236) (cid:243) (6a)
(cid:237) Œ
(cid:236) (cid:236)
(cid:237)3(cid:239) º (cid:238) Œ Œ
(cid:238) (cid:240) æ(cid:7)(cid:244),ı(cid:246)æ (cid:242) (cid:243)(cid:20)(cid:247) æ(cid:7)(cid:244) º (cid:238)(cid:15)ł (6b)
6of24
(cid:1)(cid:3)(cid:2)
(cid:3)(cid:4) (cid:5) (cid:6)(cid:5)
øbœ(cid:149)ß(cid:253)(cid:252)(cid:7)ø (cid:252)=(cid:254) (cid:255) (cid:255) œ ß(cid:253)(cid:252)
where ( is (cid:2)xed in the analysis), and . There can be either one
(cid:8)(cid:7) (cid:10)(cid:9)
ø(cid:221)œ ß (cid:9)ß
or two critical points in this system: (i) the null vector wi(cid:2) th elements , and (ii), if
(cid:11) (cid:11)(cid:12)(cid:5) (cid:13)(cid:9)
Uø œ œ ß ø œ
it exists, is the intersection of the set of points with that lie on the -nullcline
described by
(cid:7) (cid:15)(cid:14)(cid:17)(cid:16) (cid:18)(cid:20)(cid:19) (cid:22)(cid:21)
(cid:15)ß 3ø œ
(7a)
(cid:2)
(cid:11)(cid:12)(cid:7) (cid:11)(cid:12)(cid:5) (cid:23)(cid:9) (cid:7)
œ ß
and theset ofpointswith thatlieonthe -nullclinedescribedby
(cid:18) (cid:19) (cid:18)(cid:26)(cid:25)(cid:28)(cid:27)
(cid:7) (cid:24)(cid:14) (cid:29)(cid:12)(cid:30) (cid:30) (cid:31)
ß (cid:18) (cid:25) ø œ
(7b)
(cid:29)!(cid:30) (cid:30)
Realizability considerations dictate that only the positiv(cid:2) e roots need to be considered. At the
(cid:18) (cid:25)(cid:28)(cid:27) (cid:18) (cid:25)
(cid:29)(cid:12)(cid:30) (cid:30) ß
intersection point of Eq. (7a) and Eq. ((cid:2) 7b), . Thus, this second critical point
(cid:18) (cid:25)"(cid:27) (cid:18) (cid:25)
exists onlyif canachieveavalue (whichforthecurrentsetofclosurecoef(cid:2)cientsis
0# .7869). From(cid:2) Eq(cid:30) . (4#))( or Eq. (5), t(cid:2)his (cid:30)means that the second critical po(cid:29)*i(cid:30)nt (ii)can exist only if
(cid:27)(cid:22)$(cid:24)%’& (cid:18) (cid:25)"(cid:27) (cid:18) (cid:25) $(cid:24)%’& (cid:18) (cid:25)"(cid:27) (cid:18) (cid:25)
or forthese particularchoices of .
If stable, the critical points represent the possible steady-state solutions to the system,
Eq.(6). Notethatneitherofthesecriticalpointsistheso-called(cid:147)turbulent(cid:148)solution. Intheanal-
,+ -
ø(cid:204)œ
ysis of homogenous turbulent (cid:3)ow, the (cid:147)turbulent(cid:148) solution grows without bound ( ).
The reason the analysis yields an unbounded growth in this case is that the mean (cid:3)ow (cid:2)eld is
(cid:2)xed and unaffected by the turbulence, and this provides an in(cid:2)nite source of energy for the
turbulence. In practical computations, this behavior is not seen because there is a two-way
coupling between the turbulent and mean (cid:3)ow (cid:2)elds. The coupling allows for a steady-state
to be reached. Diffusion, to be introduced later in the analysis, also plays an important role in
practicalcomputations.
The stabilitypropertiesofthetwocriticalpointscanbeexaminedbylinearizingabouteach
criticalpoint. The coef(cid:2)cientmatrixofthis linearsystemis theJacobian matrix
/0 1 (cid:30) FG
. (cid:18)(cid:20)(cid:19)3254 (cid:30)(cid:25)86 7 1 & (cid:18)(cid:20)(cid:19)’254 (cid:25)86 7 &9%
ß (cid:18)(cid:20)(cid:19)(cid:12)(cid:18) (cid:25)"(cid:27)(cid:6): (cid:29);(cid:30) (cid:18) (cid:25) (cid:30) 2 (cid:25) 7 & (cid:18) (cid:25) (cid:30)(cid:22)< (cid:25)>= = & (cid:29)(cid:12)(cid:30) (cid:18) (cid:25) (cid:30) 2 (cid:25) 7 & (cid:18) (cid:25) (cid:30)(cid:22)< (cid:25)D= (cid:25)= (cid:31) (8)
4 6 4 6@?(cid:15)AC4B 6 4 6 4 6@?EA@B
A A
Todeterminethenatureofthecriticalpoints(e.g.,iftheyarestableorunstable),theeigenvalues
ofthismatrixarefoundatthecriticalpoints. Table1liststhepossibletypes. Acenterindicates
that trajectoriesorbitaround thecritical point. A stable critical pointmeans that, whensolving
the equations, the particular point can be reached; an unstable critical point or a saddle point
7of24
means that a numerical scheme will not converge to it (a saddle point is not stable in prac-
tice because orbits approach the critical point along one eigenvector, but then recede along the
eigenvectorassociated withtheunstablesolution).
H
Table1. Typeofcriticalpointasdeterminedbyeigenvaluesof .
typeofeigenvalues criticalpointtype
complexwithzero realpart center
complexwithnegative realpart stable focus
complexwithpositive realpart unstable focus
realandbothnegative stablenode
realandbothpositive unstablenode
realwithone positive,one negative saddlepoint
It can be shown that, when it exists, the second critical point (ii) of the system Eq. (6) is
I!J
alwaysa saddlepoint. WhenEq. (4)is usedfor ,the eigenvectorsat thesaddle pointare:
KL8MON KL V
MQPSRT9U RTpo
W X(cid:20)Y[Z]\ ^_X(cid:26)‘"a)bc^_X(cid:26)‘ W X(cid:20)Y(cid:12)dfehgiZ"jkb(cid:3)l9m n
J (9)
n X Y*Z]\(cid:17)qrX(cid:26)‘(cid:28)a)b ^(cid:13)X(cid:20)Y(cid:12)X(cid:26)‘ d ehg Z"jkb’^(cid:10)s!X(cid:20)Y*X(cid:26)‘ dtZ"X(cid:26)‘(cid:28)a8qvu*b*ehgcZ(cid:28)jkb j Z"X(cid:26)‘ q
U J J J U
J J J
X(cid:26)‘(cid:28)a)b]w*Zyxza)X(cid:26)‘ b d \{q(cid:8)X|‘(cid:28)a)w(cid:12)X(cid:26)‘
where U ,
J J I}J
, and . When Eq. (5) is used for , the eigenvectors at the
saddle pointare:
KL M N KL V
MSP RT U RT o
W X(cid:20)Y*Z]\(cid:22)^~X(cid:26)‘(cid:28)a)bcl W X(cid:20)Y[Z]\(cid:22)q~X(cid:26)‘(cid:28)a)b q(cid:127)X(cid:20)Y*X(cid:26)‘ dfe(cid:128)g(cid:3)Z(cid:28)(cid:129){b
J
J (10)
(cid:129) Z(cid:28)X ‘ q(cid:24)X ‘"a b(cid:130)w*Zyx(cid:132)(cid:131)5X ‘ b
U
J J
where . The effect of the saddle point on the solution trajectories
willbe shownbelowwhen phase plotsaredU(cid:13)raw(cid:135)EnU(cid:13). B(cid:136)ut themoreinterestingcriticalpointinthis
(cid:133)p(cid:134)
analysis is the (cid:2)rst (i) degenerate one ( ), which corresponds to a pseudo-laminar
I!J
solution. WhenEq. (4)orEq. (5)isused for , and ifthefollowingconditionis true,
(cid:133)E(cid:134) \ (cid:138)(cid:12)(cid:139) (cid:139)Q(cid:138)(cid:143)(cid:142)(cid:144)(cid:142)(cid:144)(cid:145)](cid:146)}(cid:138)(cid:12)(cid:147))(cid:148)(cid:144)(cid:149)(cid:143)(cid:150)
U9(cid:135){U(cid:140)(cid:136)~(cid:141) o
(cid:133) (cid:134)
(cid:135) (cid:137)
(11)
then this degenerate point is either a center or a stable point. This stability of the degenerate
8of24
criticalpointturnsouttobethecauseoftheapparentlyarbitrarysolutionbehaviordemonstrated
insection II.
InsectionC,theconditionsforwhichEq.(11)occurwillbeexploredingreaterdetail. Itwill
also be shown that stability of the degenerate point requires the second critical (saddle) point
to exist. For now, however, the assertion is made that in practice Eq. (11) is quite often true
(cid:151)3(cid:152)(cid:143)(cid:153)f(cid:154) (cid:155) (cid:151)(cid:156)(cid:152)C(cid:157)(cid:13)(cid:154) (cid:155)
for the system of equations given by Eq. (1) or Eq¡. (6). As a result, and ¡
(cid:158);(cid:159)p(cid:160) (cid:155) ¢*£ (cid:155){⁄¥(cid:158)5(cid:159)
near the cri¡tical point as well. Therefor¡e, for E£]fi;q£(cid:132).fl (4) is approximately¡ ,
(cid:158)z(cid:159) (cid:155) ¢}£ (cid:158)C£)ƒE§'¤5“}«"‹(cid:6)¤@› (cid:158)(cid:132)(cid:176)_(cid:160) (cid:155)
and for it¡ is approximately ¡ . Simi¡larly, fo(cid:159)yr–(cid:28)£]† (cid:159)y–(cid:28)E£Cflq. (5) is
¢}£ (cid:155)3⁄f(cid:158))(cid:176) (cid:158)(cid:132)(cid:176) (cid:155) ¢*£ (cid:158) ƒE§ “}«‡›·‹
¤
approximately , and for it is approximately . Using
these exprefissions, alongwithexpressions forthederivativesofEq. (4)and Eq. (5)withrespect
ƒO§
to and , the eigenvalues of the matrix in Eq. (8) near the degenerate critical point can be
¢[£
determined. These are givenin Table 2 forthetwo expressions and variouscombinations of
theircoef(cid:2)cients. Inallcases butone, thedegeneratecriticalpointis stable. (The(cid:147)center(cid:148)type
ofcriticalpointis also notdesirable inthis context, because solutionsca¡ n become lockedinan
(cid:158)(cid:12)(cid:159) (cid:155)
orbit and fail to converge. However, no known models actually use in Eq. (4), so from
nowon theanalysis focuses onlyonthemorecommoncases where thecriticalpointisstable.)
(cid:181)8¶(cid:22)•„‚(cid:22)•»” (cid:181),¶]…)‚(cid:22)‰¿(cid:190)
Table2. Eigenvaluesnear the criticalpoint,with . Thecomputationsinsection II
employedamodelthatcorrespondstothethirdalternative.
¢(cid:12)£ (cid:158)C£ (cid:158)
¤3(cid:192)_`
function( , £ fl ) type ofeigenvalues criticalpointtype
(cid:155)´⁄f(cid:158)5(cid:159)‡ˆ¯˜˙˘c«(cid:130)⁄3(cid:158)@£(cid:132)(cid:151)3(cid:152) (cid:153)
¡ , complex,near-zerorealpart center
(cid:158)5(cid:159) (cid:155)
£ fl fl
(cid:155)´⁄f(cid:158)5(cid:159)‡ˆ¯˜˙˘c«(cid:130)⁄3(cid:158)@£(cid:132)(cid:151)3(cid:152) (cid:153) (cid:158)¨(cid:159)˚(cid:201)(cid:24)«"¸(cid:26)(cid:204)y£|⁄_(cid:155) “(cid:12)¸(cid:26)(cid:204)y£
, real,both negativefor fl stablenode
(cid:201)9(cid:158)5(cid:159)(cid:20)(cid:201)(cid:24)(cid:155) (cid:158)¨(cid:159) «(cid:28)¸(cid:26)(cid:204)‡£(cid:26)⁄~(cid:155) “;¸(cid:26)(cid:204)‡£
` (cid:192)
fl complex, negativerealpartfor stable focus
(cid:155)(cid:22)⁄f(cid:158))(cid:176))ˆ5˜˙˘˝«]⁄(cid:156)(cid:158) (cid:151)(cid:156)(cid:152)C(cid:157)
¤
¡ , complex, negativerealpart stable focus
(cid:158))(cid:176) (cid:155)
fl fl
(cid:155)(cid:22)⁄f(cid:158))(cid:176))ˆ5˜˙˘˝«]⁄(cid:156)(cid:158) (cid:151)(cid:156)(cid:152)C(cid:157) (cid:158)C(cid:176)(cid:156)(cid:201)(cid:24)«"¸(cid:26)(cid:204)y£|⁄_(cid:155) “(cid:12)¸(cid:26)(cid:204)y£
¤
, real,both negativefor fl stablenode
(cid:201)9(cid:158))(cid:176)’(cid:201)(cid:24)(cid:155) (cid:158)C(cid:176) «(cid:28)¸(cid:26)(cid:204)‡£(cid:26)⁄~(cid:155) “;¸(cid:26)(cid:204)‡£
` (cid:192)
complex, negativerealpartfor stable focus
Equations (7a) and (7b) can be sketched for a case when the two nullclines intersect. In-
spection ofEq. (6a)andEq. (6b)thenshows thefollowingtoapply:
† † †
ƒE§ ¡ fl ƒE§ fl ƒE§ fl
†(cid:12)˛ «(cid:28)(cid:215)(cid:143)(cid:216) †(cid:12)˛ (cid:201) ˆ@(cid:219)D(cid:220) (cid:219) «(cid:28)(cid:215)(cid:143)(cid:216) †(cid:12)˛ (cid:220) (cid:219) «(cid:28)(cid:215)(cid:143)(cid:216)
`[ˇ(cid:209)—(cid:12)(cid:210)(cid:8)(cid:211)(cid:26)(cid:212)(cid:214)(cid:213) `[ˇ(cid:218)(cid:217) — (cid:211)(cid:26)(cid:212)(cid:214)(cid:213) (cid:192)9`[ˇ(cid:209)(cid:221)](cid:222)(cid:128)(cid:223)(cid:12)(cid:224) — (cid:211)(cid:26)(cid:212)(cid:214)(cid:213)
§ § §
(12a)
9of24
Æ(cid:12)(cid:226) Æ(cid:12)(cid:226) Æ(cid:12)(cid:226)
Æ(cid:12)ª](cid:228){(cid:229)¥(cid:230)[(cid:231)(cid:209)Ł(cid:12)Ø(cid:8)Œ(cid:26)º(cid:214)(cid:236)˙(cid:237)(cid:28)(cid:238)5(cid:239)(cid:3)(cid:240) Æ(cid:12)ª(cid:130)(cid:228)(cid:242)æ9(cid:230)[(cid:231)(cid:209)(cid:243)¨(cid:239)(cid:214)Ł(cid:143)(cid:244)*ıfŒ(cid:26)º(cid:214)(cid:236)˙(cid:237)(cid:28)(cid:238)5(cid:239)(cid:3)(cid:240) Æ(cid:12)ª](cid:228)(cid:17)(cid:246)9(cid:230)˙(cid:231)(cid:6)(cid:239)(cid:214)ı5(cid:247)hŁ(cid:143)łøŒ|º(cid:214)(cid:236)[(cid:237)"(cid:238)5(cid:239)˝(cid:240))(cid:236)
(12b)
Combining this information with the knowledge of the behavior at the critical points, a clear
mapping of the phase-plane trajectories can be drawn, as shown in Fig. 6. The eigenvectors
œOß œQ(cid:252)
( , ) at the saddle point are shown along with dividing curves (dashed lines) that separate
regions of different trajectory behavior. Clearly, any initial condition above dividing curve 1
willeventuallyend upatthe degeneratecriticalpoint.
(cid:228) (cid:226)
An actual phase-plane portrait can be constructed by computing the right hand sides of
(cid:253)
(cid:254) (cid:228)C(cid:255) (cid:254) ª (cid:228) (cid:254) (cid:226) (cid:255) (cid:254) ª (cid:228) (cid:228) (cid:226)
Eq. (6) for a large number of and values. These computed values then correspond to
(cid:253) (cid:253)
(cid:228) (cid:226)
and at each particular point in nondimensional - phase space, and the
(cid:253)
(cid:255)
trajectory(how and change withtimeorwithiteration)canbe computedas well.
(cid:229) (cid:229) (cid:229)
Æ An example is shown in Fig. 7, using Eq. (5) for (cid:0)(cid:2)(cid:1) with (cid:3)(cid:5)(cid:4) (cid:7)(cid:6) , (cid:3)(cid:5)(cid:8) (cid:10)(cid:9) (cid:9)(cid:12)(cid:11) , (cid:13)(cid:15)(cid:14) (cid:7)(cid:6)(cid:16)(cid:11)(cid:18)(cid:17)(cid:12)(cid:19) ,
(cid:229) [(cid:236) 5(cid:230)
(cid:8)
(cid:14) (cid:20)(cid:9) (cid:21)(cid:19)(cid:23)(cid:22)(cid:24)(cid:11)(cid:26)(cid:25)(cid:27)(cid:6) (cid:24)(cid:28) (typical results using Eq. (4) yield a similar behavior). The nullclines are
(cid:228) (cid:226)
shown along with the phase space trajectories. It is clear that this (cid:2)gure matches the sketch
(cid:253) (cid:229) (cid:229) (cid:230)
(cid:228) (cid:226)
in Fig. 6 in character: is a stable attractor, and the other critical point (near
(cid:253) (cid:229)(cid:13)(cid:230)[(cid:236) (cid:229) (cid:230)˙(cid:236)(cid:230) !(cid:238)
(cid:228)
(cid:29)(cid:6)(cid:30)(cid:19) , (cid:12)(cid:11) in this particular case) is a saddle point. The true turbulent solution is
(cid:253)
(cid:228) (cid:226)
obtained when grows (exponentially); in other words, the expected turbulent solution only
(cid:253)
occurs when the and values follow the trajectories in the lower right or far upper right
(cid:228) (cid:226)
parts oftheplot. This (cid:2)guredemonstrates thattherearemany regionsin themapforwhichthe
(cid:253) (cid:229) (cid:229)(cid:140)(cid:230)
solutionconverges towardthedegenerate criticalpoint .
(cid:229)
Itisalsointerestingtolookatthephasespace trajectoriesofthecase with (cid:0) (cid:1) (cid:31)(cid:6) , shownin
Fig. 8. Here, the twonullclinesdo notcross, so thereis no saddle point(i.e., no second critical
point). The trajectories now behave like the ones to the right of the saddle point in Fig. 7;
regardless of the initial condition, the solution always goes to the turbulent solution and never
tothe degenerateone.
A. AccountingforDiffusion
Inhomogeneous effects that necessarily occur in any practical calculation can be represented
(cid:255) (cid:255)
by the viscous diffusion term and the gradient diffusion models for the turbulent transport,
#"%$ #"(’
represented by (cid:18)! and &! . Since the consideration here is for high Reynolds number
(cid:3)ows,viscousdiffusioneffectscanbeneglectedintheturbulentregion,anditisonlynecessary
toadequatelyrepresenttheeffectsoftheturbulenttranportinboththekineticenergyandenergy
dissipation rate equations. It suf(cid:2)ces for the purposes here to simply assume that the transport
10of24