COMSOL 6.4 - Spanwise Rotating Turbulent Channel Flow

Spanwise Rotating Turbulent Channel Flow

In this model, the Single-Phase Flow, Elliptic Blending R-ε interface is employed to investigate fully developed turbulent flow in a spanwise rotating planar channel, as well as the subsequent flow development in a suddenly expanded section. The dynamic characteristics of the resulting flows are analyzed in detail and visualized.

Model Definition

Correctly accounting for the impact of system rotation on turbulent flows is a challenge. Eddy-viscosity models are built so as to properly react to shear and strain but, commonly, not to vorticity. Activating mends this deficiency to a certain extent. Nevertheless, the Reynolds stress tensor remains strictly proportional to the strain-rate tensor, which is too limiting in the general case. Differential Reynolds stress models of turbulence (RANS-RSM) are based on a deeper foundation, and appropriately respond both to the in-frame flow vorticity and to the frame rotation. However, even in the simplest case of a spanwise rotating turbulent channel, Wilcox R-ω and SSG-LRR underestimate the extent of the region with zero absolute vorticity, and predict laminarization at a too low rotation rate. Fortunately, the Elliptic Blending R-ε model from Ref. 1 properly blends turbulence behavior in the bulk and near-wall regions, which in turn leads to improved reaction to the system rotation.

In this example, spanwise rotating turbulent channel flow is extensively analyzed. The channel Reynolds number,

is Re = 14,000. Above, Ub is the streamwise bulk velocity and H is the height of the channel. The range of rotation numbers investigated,

is between Ro = 0 and Ro = 3 (note that Rossby number has the same denotation, Ro, although they are inverses of each other). Ω is the angular rotation rate. A Periodic Flow Condition in a geometry that is short in the streamwise direction is employed to compute the fully developed channel flow (periodic channel study). Then, the channel suddenly and symmetrically expands with a 1.5 height ratio. The deflection of the expanding jet and the asymmetry between the separation bubbles are the most apparent outcomes of the system rotation. The case Ro = 0.1 is presented below (expanded channel study).

Figure 1: The model geometry of the rotating periodic channel (top) and expanded rotating channel (bottom). The converged fully developed state of the periodic channel study is used as an inlet condition for the expanded channel study.

Assuming positive (counterclockwise, cyclonic) frame rotation and rightward flow direction, the Coriolis force acts downward. The in-frame flow vorticity is cyclonic near the top wall of the channel and anticyclonic near the bottom wall. At the anticyclonic side of the channel, turbulence is enhanced, while it is strongly suppressed at the cyclonic side. Correspondingly, the anticyclonic side is referred to as the “unstable side” and the cyclonic side as the “stable side”. In the central region of the turbulent channel flow (far away from the walls), the flow velocity possesses a constant cross-stream slope with the value of the anticyclonic in-frame vorticity, −2Ω, which exactly compensates the cyclonic vorticity, +2Ω, inherent to the system rotation.

Figure 2: Meshes used for the studies of the rotating periodic channel (top) and expanded channel (bottom). Very thin cells at the walls of the rotating channel are used to ensure accurate evaluation of traction. The initial developing region of the expanding jet is meshed fine, while the remaining near-outlet region serves mostly to reduce the effect of the outlet boundary condition.

Implementation in COMSOL Multiphysics

Figure 1 contains geometries along with boundary conditions for both studies, while Figure 2 presents the corresponding meshes. A feature is employed to impose the frame rotation. It is crucial for the efficient computation of the periodic channel, since the feature is not inherently compatible with a moving wall-based approach. An with continuation is used for quicker computation through the bunch of rotation numbers. The expanded channel flow is computed subsequently; values of variables at its inlet are taken from the converged state of the periodic channel flow. Thus, the inlet represents a fully developed turbulent channel flow in the rotating frame. The near-wall cells of the channel should be meshed very well, and the regions with separation bubbles and the ensuing wake should have very fine or at least fine meshing. The remaining part mainly serves to relieve the influence of the outlet boundary condition and might have quite coarse mesh. The pressure at the outlet is taken as if a Coriolis force of a constant velocity flow would be balanced by the cross-stream pressure gradient. Note that is activated to filter out centrifugal pressure.

The details of the implementation of the Single-Phase Flow, Elliptic Blending R-ε interface can be found in the section Theory for the Turbulent Flow Interfaces of the CFD Module User’s Guide, while specifics of the Rotating Frame feature are listed in Rotating Frame description of Theory for the Single-Phase Flow Interfaces.

Results and Discussion

Figure 3 shows velocity profiles in rotating turbulent channel flow. With increasing Ro, the profiles first become anticyclonically skewed with boosted slopes near the anticyclonic (unstable) wall and reduced slopes near the cyclonic (stable) wall. Consequently, the velocity maximum moves closer to the stable wall; for Ro in the interval (0.3, 1.1) the position of the maximum is at y > H/4, which means that it is near the stable side within one-quarter of H. At higher rotation rate, say at Ro larger than 1.1, the profiles show a clear tendency to laminarization, taking on progressively more parabolic form. At Ro = 3, laminarization is almost complete. Note that at small Ro the increase in the velocity slope near the unstable side is much less visible than the slope’s decrease near the stable wall.

Figure 3: Rotating channel flow at various rotation numbers Ro: velocity profiles. The channel is situated within [−H/2,H/2].

Indeed, Figure 4 presents the dependence of the friction velocity uτ = τw/ρ on Ro (normalized by the nonrotating value). Clearly, while uτ at the unstable side grows only slightly with Ro (by 10% at Ro = 0.1), uτ at the stable side falls abruptly (by 42% at Ro = 0.2). The picture is consistent with Ref. 1. The stable-side plunge in uτ has correct magnitude, but happens too quickly (Ro = 0.5 is an approximate rotation number when uτ reduces by 42% according to Large-Eddy Simulation).

Figure 5 illustrates flow pattern in a suddenly expanded channel section at Ro = 0.1 (rotation number is based on the parameters of the periodic channel). The jet becomes quickly deflected to the unstable side by the Coriolis force, and then slowly approaches strictly straight direction in the expanded section. Apparently, two major separation bubbles are formed. The one at the unstable wall is quite compact, shorter than in nonrotating case, with easily distinguishable boundary. Meantime, that on the stable wall is very extended in the streamwise direction and has long fuzzy streaks, which make it hard to determine reconnection point visually. Indeed, according to Ref. 1, the Elliptic Blending R-ε model accurately reproduces behavior on the unstable side, but is less reliable at the stable side of the expanded channel setup. Allegedly, imperfect sensitivity of the model to the laminarization tendency is responsible for this. It might be suggested that true longitudinal length of the stable-side bubble is approximately half the length predicted by the current version of the Elliptic Blending R-ε. The major bubbles are accompanied by small counterrotating secondary “corner” bubbles (pushed to the corners of the expanded channel).

Figure 4: Rotating channel flow at various rotation numbers Ro: uτ(Ro)/uτ(0) is presented; y is within [−H/2,H/2].

Figure 6 shows friction coefficients along the stable and unstable walls of the expanded section,

with streamwise length normalized by the step height hstep. Notice that even on unstable side it takes quite a distance before friction coefficient magnitude settles to a constant value. On the stable side Cf is negative for very long stretch, which indicates that the turbulence model overestimates rotation-induced tendency to laminarization. Thus, flow recovery to the new fully developed state is significantly delayed. The friction coefficient variations inside the secondary bubbles are also revealed.

Figure 5: Suddenly expanded rotating channel. The jet is deflected down by the Coriolis force, the negative maximum of the y-velocity is twice higher than its positive maximum, and two different separation regions are formed. The separation regions are in color, and each consists of a major bubble and a counterrotating secondary“corner bubble”. Re = 14,000 and Ro = 0.1(parameters for the periodic channel).

Figure 6: Friction coefficients along the walls of the expanded channel at Re = 14,000 and Ro = 0.1(parameters for the narrower section). The length is normalized by the step height hstep = (Hexp − H)/2.

Figure 7 illustrates contours of reduced pressure, which are almost perpendicular to the Coriolis force and are strongly distorted by the minima associated with pressure losses in the separation regions. The trough on the stable side is especially strong. Velocity magnitude is included too. The recirculation streamline line at the edge of the bottom bubble is well distinguishable, while the recirculation streamline at the edge of the top bubble is sticking to the wall almost in parallel fashion.

Figure 8 demonstrates Reynolds stresses. All the components possess much higher maxima near the unstable side. In particular, the region on the stable side where vv is suppressed extends downstream a lot. This might be one of the reasons for the observed behavior of Elliptic Blending R-ε on the stable side, since turbulent diffusion is proportional to vv in the Daly–Harlow turbulent diffusion model. The shear layers of uu and ww are more elongated at the stable side than at the unstable side, while the opposite applies to the shear layers of uv and vv. The boundary layer of uu (with peak values) is well pronounced near the unstable wall, while it does not develop near the stable wall at the distances shown.

Figure 7: Levels of the reduced pressure with pronounced troughs (especially in the top bubble), and velocity magnitude in the expanded channel. Recirculation lines have purple color. Notice that the unstable side recirculation line terminates abruptly, while the stable side recirculation line leans to the wall smoothly.

Summary and Outlook

The Elliptic Blending R-ε Reynolds stress model implemented in COMSOL Multiphysics can predict most features of rotating turbulent channel flow by properly accounting for the influence of rotation and near-wall effects. Both periodic rotating channels and suddenly expanded rotating channels can be analyzed with the approach.

In rotating periodic channel flow, velocity profiles are captured correctly at various rotation numbers. Friction velocities at the anticyclonic (unstable) side are quite reliable. Due to oversensitivity to rotation, friction velocities at the cyclonic (stable) side abruptly fall with slightly nonzero rotation number.

In rotating expanded channel flow, the major separation bubble and the friction coefficient are predicted well at the unstable side. At the stable side, laminarization tendencies are too strong, which results in an overpredicted bubble size with very high recirculation length.

Figure 8: Reynolds stresses in the expanded channel section. Top to bottom: shear stress uv, diagonal components (common color legend)- streamwise uu, cross-stream vv, spanwise ww. There are clear bands of increased stresses in the shear layers, which are stronger at the unstable side. uu has very pronounced region of high values in the boundary layer along the unstable wall, but along the stable wall the boundary layer is not visible at the distances shown.

Subsequently, the friction coefficient recovers, approaching its fully developed value slowly, too.

To summarize, Elliptic Blending R-ε is superior to other RSMs in prediction of near-wall effects, such as friction and Reynolds stress profiles. For example, it is able to correctly reproduce near-wall peaks of the diagonal streamwise component of Rij. The model is a good tool for capturing the effect of system rotation on turbulence development, although its predictions at the stable side (lying opposite to the direction of the Coriolis force) should be taken with caution.

Reference

1. R. Manceau, “Recent progress in the development of the Elliptic Blending Reynolds-stress model,” Int. J. Heat Fluid Flow, vol. 51, pp. 195–220, 2015; doi.org/10.1016/j.ijheatfluidflow.2014.09.002.

CFD_Module/Single-Phase_Flow/rotating_turbulent_channel

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > > ε (spf).

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file rotating_turbulent_channel_parameters_1.txt.

Parameters 2

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file rotating_turbulent_channel_parameters_2.txt.

Geometry 1

Periodic channel

In the toolbar, click

In the window for , type Periodic channel in the text field.

Locate the section. In the text field, type 2*l_delta.

In the text field, type H.

Locate the section. From the list, choose .

In the text field, type -L_i-l_delta.

Locate the section. Find the subsection. Click .

In the dialog, type Channel in the text field.

Click .

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type -L_i-l_delta.

In the text field, type -H/2.

Locate the section. From the list, choose .

In the text field, type -L_i-l_delta.

In the text field, type H/2.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type L_i.

In the text field, type H.

Locate the section. From the list, choose .

In the text field, type -L_i/2.

Rectangle 3 (r3)

Right-click and choose .

In the window for , locate the section.

In the text field, type H/2.

In the text field, type H_exp.

Locate the section. In the text field, type H/4.

Rectangle 4 (r4)

Right-click and choose .

In the window for , locate the section.

In the text field, type L_1.

Locate the section. In the text field, type L_1/2.

Rectangle 5 (r5)

Right-click and choose .

In the window for , locate the section.

In the text field, type L_2.

Locate the section. In the text field, type L_1+L_2/2.

Definitions

Expanded channel

In the toolbar, click

In the window for , type Expanded channel in the text field.

Locate the section. Under , click

In the dialog, select in the list.

Click .

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Turbulent Flow, Elliptic Blending R-ε (spf)

In the window for ε, locate the section.

From the list, choose .

Locate the section. Select the checkbox.

Select the checkbox.

Fluid Properties 1

In the window, under > ε (spf) click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type rho_i.

From the μ list, choose . In the associated text field, type mu_i.

Initial Values 1

In the window, click .

In the window for , locate the section.

Specify the u vector as


U_i	x

Periodic Flow Condition 1

In the toolbar, click

and choose .

Select Boundaries 1 and 7 only.

In the window for , locate the section.

From the list, choose .

In the

text field, type rho_i*U_i*H*1[m].

Pressure Point Constraint 1

In the toolbar, click

and choose .

Select Point 3 only.

Rotating Frame 1

In the window, click .

In the window for , locate the section.

In the Ωf text field, type Ro.

Definitions

Domain Point Probe 1

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to -L_i-l_delta.

In row , set to H/2.

Point Probe Expression 1 (ppb1)

In the window, expand the node, then click .

In the window for , type utau_c in the text field.

Locate the section. In the text field, type sqrt(-spf.nu*ppr(uy)).

Domain Point Probe 2

In the window, under > right-click and choose .

In the window for , locate the section.

In row , set to -H/2.

Point Probe Expression 1 (ppb2)

In the window, expand the node, then click .

In the window for , type utau_ac in the text field.

Locate the section. In the text field, type sqrt(spf.nu*ppr(uy)).

Mesh 1

Mapped 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Distribution 1

Right-click and choose .

Select Boundaries 2, 3, 5, and 6 only.

In the window for , locate the section.

In the text field, type 1.

Distribution 2

In the window, right-click and choose .

Select Boundaries 1, 4, and 7 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 200.

In the text field, type 15.

Select the checkbox.

Click

Study 1

Step 1: Wall Distance Initialization

In the window, under click .

In the window for , locate the section.

In the column of the table, under , clear the checkbox for ε.

Step 2: Stationary

In the window, click .

In the window for , locate the section.

In the column of the table, under , clear the checkbox for ε.

Click to expand the section. Select the checkbox.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
Ro (Rotation number)	range(0,0.025,0.2) range(0.3,0.1,1.4) range(1.5,0.25,3)

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Under , click

In the dialog, select in the list.

Click .

In the toolbar, click

Turbulent Flow, Elliptic Blending R-ε 2 (spf2)

In the window, under click ε.

In the window for ε, locate the section.

From the list, choose .

Locate the section. Select the checkbox.

Select the checkbox.

Fluid Properties 1

In the window, under > ε 2 (spf2) click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type rho_i.

From the μ list, choose . In the associated text field, type mu_i.

Initial Values 1

In the window, click .

In the window for , locate the section.

Specify the u vector as


U_i*H/H_exp	x

Rotating Frame 1

In the window, click .

In the window for , locate the section.

In the Ωf text field, type Ro_2.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 7 only.

In the window for , locate the section.

Click the button.

Specify the u0 vector as


withsol('sol1',u,setval(Ro,Ro_2))	x
withsol('sol1',v,setval(Ro,Ro_2))	y

Locate the section. Click the button.

Specify the R0 matrix as


withsol('sol1',uu,setval(Ro,Ro_2))	withsol('sol1',uv,setval(Ro,Ro_2))	0
withsol(’sol1’,uv,setval(Ro,Ro_2))	withsol('sol1',vv,setval(Ro,Ro_2))	0
0	0	withsol('sol1',ww,setval(Ro,Ro_2))