12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

ON THE NEW UNSTABLE MODE IN THE BOUNDARY LAYER FLOW OF
SUPERCRITICAL FLUIDS

Benjamin Bugeat
Process and Energy Department
Delft University of Technology
Leeghwaterstraat 39,
2628 CB Delft, The Netherlandsb.bugeat@tudelft.nl

Pietro Carlo Boldini
Process and Energy Department
Delft University of Technology
Leeghwaterstraat 39,
2628 CB Delft, The Netherlands
p.c.boldini@tudelft.nl

Rene Pecnik
Process and Energy Department
Delft University of Technology
Leeghwaterstraat 39,
2628 CB Delft, The Netherlands
r.pecnik@tudelft.nl

ABSTRACTwave. This mode was furthermore different from the acoustic unstable Mack’s modes observed in compressible boundary layers. In their study, the authors considered temperature gradients generated by viscous heating. The new mode
was then detected when sufficient viscous heating, controlled
by the Eckert number, made the temperature profile cross the
Widom line. This line, in the supercritical region of the thermodynamic space, divides the gas-like from the liquid-like behaviour of a fluid, and is then also called pseudo-boiling line
(Banuti, 2015). It can be seen as an extension of the coexistence line between liquid and gas but greatly differs fromit as thermodynamic properties continuously vary across it.
Still, steep variations of the supercritical fluid properties areobserved around the Widom line, in particular the density andthe viscosity.
The aim of the present paper is to shed light on the origin
and the nature of the new unstable mode by understanding therole played by the Widom line. A linear stability analysis using the Navier-Stokes equation is first carried out on a Van der
Waals fluids, allowing us to recover the qualitative behaviourof Ren et al. (2019). Inviscid stability calculations are thenperformed in order to show the inflexional nature of this instability. Subsequently, we theoretically show that a generalised
inflexion point (GIP) can appear when the temperature profile
crosses the Widom line is crossed, and that it is fundamentally
linked to the presence of a minimum of kinematic viscosity. A
discussion on the universality of this instability, beyond supercritical fluids, is then opened.

Ren et al. (2019) recently studied the stability of theboundary layer flow over a flat plate for supercritical CO2.
While only one unstable mode usually exists for boundary
layer flows, the authors found an additional unstable mode,
whose origin has so far not been identified. In the present
work, we carry out a stability analysis in the general case of
a fluid following the Van der Waals equation of state and flowing over a heated flat plate in the limit of zero Eckert number.
In this framework, the second unstable mode is also recovered, ruling out an acoustic origin. From the Rayleigh equation derived in the presence of density gradients, a generalised
inflection point (GIP) criterion of instability exists, similar tothat of fully compressible flows. Inviscid stability calculations
confirm the existence of an unstable mode in the presence ofa GIP, which is linked to the additional second mode found at
finite Reynolds numbers. A theoretical analysis is then carriedout by approximating the momentum equation for a base flow
exhibiting strong gradients of dynamic viscosity. It is shownthat the origin of the GIP, and hence the additional unstablemode, is associated with a minimum of kinematic viscosity at
the Widom line. The universality of this result beyond supercritical fluids is eventually discussed.

CONTEXT AND OBJECTIVES
Above certain conditions of temperature and pressure, fluids reach a supercritical state in which the distinction between
liquid and gas states no longer exists. These fluids have important applications in industrial processes, in particular in power
cycles (Iverson et al. (2013)) for which accurate models of heattransfers are required. Heat transfers depend dramatically on
the flow regime; yet, little is known about the transition toturbulence in supercritical fluids. Ren et al. (2019) recently
studied the stability of the boundary layer over a flat platefor supercritical CO2. The authors discovered a new linearly
unstable mode in addition of the usual Tollmien-Schlichting

APPROACH
The scope of this study is limited to a supercritical fluid
flowing over a flat plate with an Eckert number set to zero(equivalently, the Mach number is also zero). This allows us
to ignore acoustic effects, which were taken into account by
Ren et al. (2019). The physical interpretation of the results isthen simplified while the equations are less cumbersome. To

112th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
VISCOUS STABILITY RESULTS
The results of the linear stability analysis are shown infigure 2, in which neutral curves are plotted in the Reynolds
number-frequency space (Re, ω). When the wall temperature
is below the pseudo-boiling temperature, one and only one unstable mode, termed mode A, is found. It is the TollmienSchlichting wave reminiscent of the incompressible case in theabsence of temperature gradient. When Tw is increased andthe Widom line is crossed, a new mode, termed mode B, appears in addition to mode A. It can be seen that it results fromthe destabilisation of a distinct eigenvalue from mode A since
the two modes can be unstable at the same values of Re andω. Mode B is further destabilised as Tw is increased, with thecritical Reynolds number of instability dramatically decreasing and the magnitude of the growth rate increasing.
These results are qualitatively similar to those of Ren
et al. (2019). This shows that mode B is not associated withacoustic effects since Ec = 0 in the present calculation. Moreover, while the authors studied supercritical CO2 with nonanalytical thermodynamic functions obtained from look-up tables, we here see that the simple Van der Waals equation ofstates associated with analytical diffusion laws produces thesame observations. This indicates some degree of robustnessof this new mode and shows that this simpler and more universal model contains the physical ingredients responsible forits growth, offering a framework that may be more suited for a
fundamental analysis of this instability.

Figure 1: Base flow profiles for different wall temperatures.
Left: temperature. Profiles that crosses the pseudo-boilingtemperature, indicated by the vertical black line, are called
transcritical. Right: dynamic viscosity.

INVISCID STABILITY ANALYSIS
A linear inviscid stability analysis is carried out to determine the nature of mode B. When the wall temperature is
below the pseudo-boiling temperature (Tw /Tpb < 1), no unstable mode is found. However, as soon as Tw /Tpb > 1, positive
growth rate are detected (figure 3). In addition to the evidenceof the inflexional nature of mode B, this result represents a
new strong indication that this mode is intrinsically linked tothe crossing the Widom line by the base flow profile. Wall temperature very close to the threshold of the transcritical regime
are indeed studied here, with Tw /Tpb = 1.04 being equivalent
to the Tw /T∞ = 1.15 case presented in figure 2.
The inviscid growth rate increases as the wall temperature increases in the transcritical regime (figure 3). The cut-off
frequency of the instability, at which the perturbation becomesneutral, also increases with Tw , extending the range of unstablefrequencies. Note that negative growth rate are found becauseof the complex mapping mentioned before. Real mappings
would only produce neutral modes if no unstable mode exists
given the existence of pairs of eigenvalues and their complex
conjugate in this case.
It has been seen that the growth rate is small when Tw /Tpb
is small. In viscous stability calculation, due to viscous damping, this means that this mode could be absent at Reynolds
number of the order of those usually studied in boundary layer.
In figure 2, a very small portion of the parameter space wasunstable at Tw /T∞ = 1.15, while the neutral curve invaded alarge region at Tw /T∞ = 1.20. This behaviour is expected for
an inflexional instability: the larger the inviscid growth rate is,the more viscous damping is needed to stabilise it, hence thesmaller the critical Reynolds number must be. But another viscous effect is actually at play in addition to viscous damping.
In figure 4, the growth rate of mode B is plotted for different
Reynolds numbers and is compared to the inviscid case. Larger
growth rates are obtained for smaller Reynolds numbers presented, indicating a destabilising effect of the viscosity (fortrigger non-ideal effects in the absence of viscous heating generated by finite Eckert numbers, the flat plate is maintained at
a prescribed temperature Tw larger than the freestream temperature T∞ . The Van der Waals equation of state is used alongwith the dynamic viscosity and thermal conductivity laws forsupercritical fluids proposed by Jossi et al. (1962) and Stiel& Thodos (1964), respectively. The reduced pressure being
fixed at p/pc = 1.08, the pseudo-boiling temperature is foundto be equal to Tpb /T∞ = 1.11. Two wall temperatures below
and two above Tpb are considered (figure 1). In the two latter cases, the temperature profile then crosses the Widom line:
there is a wall-normal location at which the base flow verifies
T̄ = Tpb . These cases are referred to as transcritical.
The Navier-Stokes equations are solved numerically inthe limit of the boundary layer approximation in order to obtain steady, self-similar base flows for different wall temperatures (as just shown in figure 1). A local stability analysis
is then carried out by decomposing small perturbations into
normal modes. This involves linearising the Navier-Stokes
equations written in the limit of zero-Eckert number. In thisframework, density may only vary due to the effect of temperature, while the pressure is purely of hydrodynamic nature.
An eigenvalue problem is eventually solved using a pseudospectral method (Orszag, 1971), from which the spatial growthrate, noted −αi , of each mode is obtained.
An inviscid linear stability solver is also used in this workin order to study the stability of the flow in the limit of infinitely large Reynolds numbers. It is based on the Rayleighequation extended to a varying-density base flow. Following
Boyd (1985), a parabolic complex-mapping is used to avoid
the well-known critical layer singularity occurring for neutralmodes. This allows us to accurately calculate small values of
the growth rate.

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 4: Effect of the Reynolds number on the growth rate ofmode B at Tw /T∞ = 1.2. The Re → +∞ line is obtained frominviscid analysis.

Figure 5: Quantity associated with the generalised inflexion
point (GIP) criteria, calculated for different wall temperatures.
A GIP is detected when this quantity is equal to zero.

small enough Reynolds numbers, viscous damping eventuallytakes over and stabilises the flow, as already seen in figure 2).
This behaviour will be discussed at the end of the paper as it
opens interesting fundamental questions.

Figure 2: Stability diagram in the Reynolds-frequency spacefor different wall temperatures. Colours indicate the value ofthe growth rate of mode B when it is unstable. The red dashedline shows the neutral curve of mode A, which is unstable ineach case. Note that the centre and right plots corresponds totranscritical cases.

ORIGIN OF THE ADDITIONAL UNSTABLE
MODE
Discarding viscous effects, the inviscid linearised NavierStokes equation can be recast as the Rayleigh equation, solvedin the previous section. In the presence of density gradients,
it can be shown that a necessary condition of instability is theexistence of a GIP in the base flow profile, defined as the location where the quantity (ρ̄ ū0 )0 is zero (the superscript 0 stands
for wall-normal derivative). This quantity is plotted in figure
5. A GIP is found only in the transcritical cases, indicating thatsome features associated with the Widom line are responsiblefor it.
To investigate this further, the base flow is examined in
the vicinity of the pseudo-boiling point where strong, local
gradients of dynamic viscosity are expected. In the boundary
layer flow over a flat plate, viscous stresses are usually balanced with advection effects, the latter decreasing and eventually vanishing at the wall. However, if the viscosity is allowed
to vary, a term associated with viscosity gradients appears andcan, in some region, even be dominant compared to advection.
This is more likely to occur near the wall where advection is
weak and, of course, for fluids exhibiting large viscosity variations. Under the assumption that viscosity effects dominate the

Figure 3: Inviscid linear stability analysis: growth rate for different wall temperatures close to the pseudo-boiling temperature.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
balance of streamwise momentum of the boundary layer equations (the validity of this assumption will be discussed later),
we can then writeµ̄ 0 ū0 + µ̄ ū00 ' 0.

(1)

Furthermore, one can simply rewrite (ρ̄ ū0 )0 as(ρ̄ ū0 )0 = ρ̄ ū00 + ρ̄ 0 ū0 .

(2)

Combining equations (1) and (2) and making use of µ̄ 0 =T̄ 0 (∂ µ̄/∂ T ) and ρ̄ 0 = T̄ 0 (∂ ρ̄/∂ T ) at constant pressure resultin the following approximation:(ρ̄ ū0 )0 ' ρ̄ ū0 T̄ 0!
1 ∂ ρ̄
1 ∂ µ̄−.
ρ̄ ∂ T p µ̄ ∂ T p

Figure 6: Kinematic viscosity of different supercritical fluids
around the pseudo-boiling point (p/pc = 1.1).

(3)

Under this approximation, it can be seen that a GIP will appear as soon as the quantity in brackets in the right-hand side iszero. Importantly, this quantity is only related to the fluid properties rather than the flow properties. Noticing that the derivative of the kinematic viscosity ν can be introduced, equation(3) eventually reduces to(ρ̄ ū0 )0 ' −ρ̄ ū0 T̄ 0

∂ ν̄
∂T p

Figure 7: Comparison between the quantity plotted in figure 5
at Tw /T∞ = 1.20 and the approximation proposed in equation(3). The horizontal black dotted-line shows the location of the
Widom line.

(4)

This implies that the existence of an extremum of ν is, under
the working assumptions, a sufficient condition for a GIP toexist (which is, in turn, a necessary condition for an inviscid
instability to exist). The profile of kinematic viscosity is plotted against temperature for real, supercritical fluids in figure 6.
All of the presented fluids - O2 , CO2 and toluene - admit a minimum of ν at the pseudo-boiling line. To verify that this minimum is indeed the key ingredient at the origin of the GIP insupercritical fluids, the approximation obtained in equation (4)is tested in figure 7. It is verified near the wall where viscosity
effects are always dominant since advection effects go to zero
given the no-slip condition. An excellent agreement is also observed in the immediate vicinity of the pseudo-boiling point,where the above derivation was indeed intended to hold because of the strong gradients of viscosity in supercritical fluids
in this thermodynamic region. Crucially, it is at this locationthat the GIP is found. Therefore, equation (4) indeed provides
the interpretation that the GIP originates from the extremumof kinematic viscosity observed at the Widom line. More precisely, it can be shown that a minimum of kinematic viscosity
is in fact required in order to satisfy the Fjørtoft’s theorem extended to varying density fluids. Ultimately, these results suggest that any supercritical fluids, whose flow can be describedby the boundary layer approximation and whose temperature
profile crosses the Widom line, will exhibit an inviscid instability associated with a GIP located around the Widom line.

which supercritical fluids satisfy. These conditions are thepresence of strong local gradients of dynamic viscosity coupled with a minimum of kinematic viscosity; no features thatare specifically linked to supercritical fluids are eventually invoked. The universality of this result could motivate the investigation of other flows exhibiting these properties. Somefluid mixtures, such as methanol and toluene, are for example known to exhibit a minimum of kinematic viscosity forintermediate mole fractions (McAllister, 1960). The stabilityof shear flows involving layers of fluids of different viscosities was studied by Yih (1967) who showed that viscosity gradients could cause instability at low Reynolds numbers. Ern
et al. (2003) investigated this case in the presence of a finite
interface between the fluids. This work could be revisited for
fluid mixtures featuring a minimum of kinematic viscosity inwhich an inviscid instability may enrich the dynamics of thesystem. In this case, it can be noted that the viscosity profile is driven by the concentration of a fluid species, which isa scalar following a transport-diffusion equation analogous to
the temperature studied in the present paper. The two cases are
therefore very similar from a theoretical point of view.
Different shear flows could also be considered. Plane
Couette flow is, for example, expected to lead to similar results
as those presented in this work since it strictly satisfies equation (1). This could also constitute a more canonical framework than the flat-plate boundary layer to study mode B. Conversely, plane Poiseuille flow would tell a different story as thestreamwise pressure gradient driving the base flow would preclude the use of equation (1) and, subsequently, the validity of
the criterion of minimum kinematic viscosity.

DISCUSSION
The additional unstable mode found in the flat-plateboundary layer flow of supercritical fluids has been examined.
It has been shown to result from an inviscid instability associated with a GIP that is generated under certain conditions

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Ern, Patricia, Charru, François & Luchini, Paolo 2003 Stability
analysis of a shear flow with strongly stratified viscosity.
Journal of Fluid Mechanics 496, 295–312.
Iverson, Brian D, Conboy, Thomas M, Pasch, James J &Kruizenga, Alan M 2013 Supercritical co2 brayton cycles
for solar-thermal energy. Applied Energy 111, 957–970.
Jossi, John A, Stiel, Leonard I & Thodos, George 1962 The
viscosity of pure substances in the dense gaseous and liquid
phases. AIChE Journal 8 (1), 59–63.
Mack, Leslie M 1984 Boundary-layer linear stability theory.
Tech. Rep.. California Inst of Tech Pasadena Jet Propulsion
Lab.
McAllister, RA 1960 The viscosity of liquid mixtures. AIChE
Journal 6 (3), 427–431.
Orszag, Steven A 1971 Accurate solution of the orr–sommerfeld stability equation. Journal of Fluid Mechanics
50 (4), 689–703.
Ren, Jie, Marxen, Olaf & Pecnik, Rene 2019 Boundary-layer
stability of supercritical fluids in the vicinity of the widom
line. Journal of Fluid Mechanics 871, 831–864.
Stiel, Leonard I & Thodos, George 1964 The thermal conductivity of nonpolar substances in the dense gaseous and liquid
regions. AIChE Journal 10 (1), 26–30.
Yih, Chia-Shun 1967 Instability due to viscosity stratification.
Journal of Fluid Mechanics 27 (2), 337–352.

Finally, investigating the presence of larger growth at finite Reynolds number than at infinitely large one (observed
earlier in figure 4) would be an interesting follow-up study.
This phenomenon is not unusual as it occurs, for example,
with the first mode of a supersonic boundary layer which, in
addition of being inviscid unstable, bears the viscous destabilisation mechanism of the incompressible TS waves (see forexample Mack (1984)). The TS-viscous mechanism is alreadyassociated with mode A in supercritical fluids. This opens fundamental questions regarding the destabilisation mechanismof mode B at finite Reynolds number, which would merit a
dedicated investigation. Given the aforementioned similarity
with miscible fluids, a potential approach could be to study if
the low Reynolds number instability existing in miscible fluids
could, to some extent that should be understood, be invoked to
explain this behaviour.

REFERENCES
Banuti, DT 2015 Crossing the widom-line–supercriticalpseudo-boiling. The Journal of Supercritical Fluids 98, 12–16.
Boyd, John P 1985 Complex coordinate methods for hydrodynamic instabilities and sturm-liouville eigenproblems withan interior singularity. Journal of Computational Physics
57 (3), 454–471.

5