Application of the Turbulent Potential Model to Unsteady Flows and Three-Dimensional Boundary Layers

The turbulent potential model is a Reynolds-averaged (RANS) turbulence model that is theoretically capable of capturing nonequilibrium turbulent flows at a computational cost and complexity comparable to two-equation models. The ability of the turbulent potential model to predict nonequilibrium turbulent flows accurately is evaluated in this work. The flow in a spanwise-driven channel flow and over a swept bump are used to evaluate the turbulent potential model’s ability to predict complex three-dimensional boundary layers. Results of turbulent vortex shedding behind a triangular and a square cylinder are also presented in order to evaluate the model’s ability to predict unsteady flows. Early indications suggest that models of this type may be capable of significantly enhancing current numerical predictions of turbomachinery components at little extra computational cost or additional code complexity.


INTRODUCTION
Three-dimensional boundary layers and unsteady vortex shedding represent situations that are common in many engineering flows.These flow situations are examples in which the turbulence does not have time to reach equilibrium with the mean flow.Flows on aircraft wings, inside curved ducts or bends, on rotating disks, propellers, junctions of stationary walls of fluid machinery are good examples of non-equilibrium flows.Previous research has indicated that two equation models, which assume that the Reynolds-stress tensor is aligned with the velocity-gradient tensor, are not accurate for non-equilibrium flows.When such models were applied to three-dimensional shear-driven or pressure driven flows (non-equilibrium flows), they yielded poor results (Fannelop et. al. 1975, Bradshaw et. al. 1996).Ölçmen & Simpson (1993) reviewed the performance of more complicated algebraic (or nonlinear) eddy-viscosity models and concluded that none of the models could perform well in all the cases studied.The models that accounted for the anisotropy of the eddy viscosity in general performed better.However, the anisotropic constants had to be changed for different flows.
Reynolds stress transport models can predict nonequilibrium flows for the fact that they do not hypothesize constitutive relation between the Reynolds stress and mean flow gradients.However the Reynolds stress transport models tend to be more expensive and numerically stiff.Recently, Durbin (1993a) proposed an elliptic Reynolds stress model that could reproduce some of the features observed by Moin et. al. (1990) in a three-dimensional channel flow.For the unsteady vortex shedding behind bluff body, Franke and Rodi (1991) compared the ability of different models to predict turbulent vortex shedding from a rectangular cylinder.Their conclusion is that some k-ε models do not predict the right shedding frequency but more expensive Reynolds stress transport models can produce results in good agreement with the experiments.
The turbulent potential model should be able to capture nonequilibrium situations at a cost comparable to two equation models.The model does not hypothesize an explicit relationship between the turbulence and the mean flow and it is the least expensive RANS model that is theoretically capable of capturing non-equilibrium turbulent flows.Non-equilibrium occurs when the mean flow changes rapidly either in time or spatially along the flow direction more quickly than the turbulence can respond.Turbulent non-equilibrium is a persistent characteristic of large scale unsteady flows such as bluff body vortex shedding, and accelerated flows such as three-dimensional boundary layers.Detailed information on the development and formulation of the turbulent potential model can be found in Perot (1999) and Perot and Wang (1999).The particular emphasis in this work is evaluating and enhancing the model's performance in threedimensional boundary layers and unsteady flows.

TURBULENT POTENTIAL MODEL
The transport equations that constitute the turbulent potential model are summarized below.
where φ and ψ are the scalar and vector potentials respectively.

NUMERICAL METHOD
The numerical method uses an unstructured staggered mesh scheme which can conserve mass, momentum, and kinetic energy to machine precision.The turbulence quantities are advected using an upwinding scheme to guarantee positivity constraints.The model integrates up to the wall, so wall functions are not used, but the first grid point should be in the laminar sub-layer to obtain accurate predictions.The details of this numerical method, including accuracy analysis and conservation properties are discussed in Perot (2000) and Perot & Zhang (1999).

THREE-DIMENSIONAL BOUNDARY LAYERS
The turbulent body force potential model's ability to predict three-dimensional boundary layers is studied using two different three-dimensional boundary layers.The first was a spanwise driven channel flow.In this flow, a large spanwise pressure gradient is suddenly applied to a fully developed channel flow.The pressure gradient produces a horizontally homogeneous flow (Figure 1) that develops in time.The mean velocity vector skews with height as the flow turns most quickly near the wall.Spanwise driven channel flow at bulk velocity Reynolds number of 3300 was simulated by Moin et al (1990) using direct numerical simulation (DNS).Durbin (1993a) also modeled this flow with a Reynolds stress transport model.A two-equation RANS model would not be expected to produce good results for this flow, since the turbulence effects would appear instantly in the spanwise direction rather than taking some time to develop.The turbulent potential model predictions of the turbulence quantities are improved by adding an extra production term to the dissipation evolution equation.We add a term ω ψ × where ω is the mean vorticity and ψ is the turbulent potential.This term is zero in two-dimensional problems and has some similarities with the classic production term, ω ψ ⋅ .This extra term captures the breakup of turbulent structures due to threedimensional acceleration, and the resulting drop in turbulence intensities.
Figure 1.Schematic of two-dimensional channel flow subjected to spanwise pressure gradient The DNS data is for a channel flow where the suddenly applied spanwise pressure gradient is ten times the streamwise pressure gradient.The data is given at the time of the spanwise pressure gradient application (t=0) and at time increments of 0.3 thereafter up to a nondimensional time of 0.9.All quantities are non-dimensionalized by the channel half-width and the initial shear velocity.The dimensionless viscosity is 0.00556 = 1/180.The spanwise velocity increases rapidly during this time and reaches the same order of magnitude as the streamwise velocity.The turbulence cannot be considered to be in equilibrium at any time during this simulation because the eddy turn over time is of the order of one.Curiously in this flow, the turbulent intensities decrease initially.This is counterintuitive, as one would expect an increase in the turbulent intensity because of the increase in magnitude of the shear due to the additional pressure gradient.The mean streamwise velocity is shown in figure 2. The symbols are the DNS data and the solid lines are the model predictions.This convention remains the same in all future graphs involving the spanwise driven channel flow.It is seen clearly that, in both the DNS data and the model predictions, there is not any significant change in the streamwise velocity during the initial development of the flow.
The mean spanwise velocity is shown in figure 3. The mean spanwise velocity is roughly half the mean streamwise velocity at the final measurement time.The agreement with the DNS data is good.Visual comparison with the fully turbulent streamwise velocity profile makes it clear that the spanwise velocity boundary layer is essentially laminar at these early times.The turbulence may just have begun to effect the spanwise velocity at the final time.The standard ε and other models, which assume equilibrium, would have applied the full turbulent eddy viscosity to the spanwise velocity and thereby caused the spanwise boundary layer to grow much more rapidly as a turbulent boundary layer.Since the turbulence potential model is non-equilibrium in nature, the spanwise velocity predictions of the model are in good agreement with the DNS data.In flows with a single inhomogeneous direction, such as spanwise driven channel flow, there is a direct correspondence between the turbulent potentials and some of the Reynolds stresses allowing a direct comparison to take place.In the more complex flows shown later comparison with the turbulence quantities is theoretically possible but highly impractical.The v u ′ ′ and w v ′ ′ shear stresses are shown in figure 4. The upper set of curves are the v u ′ ′ shear stress, and it is this stress that directly influences the evolution of the streamwise velocity.The curves actually proceed from top to bottom as t= 0.0, 0.3, 0.9, 0.6, showing that this stress initially drops very slowly, between t=0.3 and t=0.6 it drops much more rapidly, and then after that it begins to increase towards its initial value.The model predictions decrease monotonically but display similar qualitative behavior.The initial drop is very small, and speeds up at later times.If the calculation allowed to proceed past t=0.9 the stress begins to increase as with the DNS data.While the predictions do not match exactly, the discrepancy appears to be a simple time lag in the predictions since the model predictions at t=1.2 (not shown) are close to DNS data at t=0.9.The mean flow predictions for the streamwise velocity show that, at least at these early times, that the time-lag defect is not fundamentally important in predicting the mean flow.
The lower group of curves in figure 4 are for the shear stress.The DNS data shows that this stress starts at zero and then increases monotonically as time proceeds.After a long time, this stress can be expected to be ten times larger than the v u ′ ′ shear stress.However, at these early times it is relatively small.The model predictions closely match the DNS data at these early times.The small magnitude of this stress at early times leads to the essentially laminar velocity profile for spanwise mean velocity.While the turbulent potential model predicts the spanwise shear stress well at these early times and is expected to perform equally well at later times, it is also true that the spanwise shear stress is relatively small at these early times.The development of the spanwise velocity profile is therefore essentially laminarlike at these early times and any non-equilibrium model (that starts with small spanwise shear stresses) would also be able to predict the spanwise development at these early times quite well.
Hence, a more stringent test of the turbulent potential was performed by predicting the flow over an infinite swept bump.If sidewall effects and the boundary layer growth on the top wall are neglected, this flow can be computed as a two-dimensional problem (it is still a three-dimensional boundary layer, however).The flow over a swept bump is an interesting test case because it test the model's ability to predict a turbulent boundary layer that is subjected to both streamwise pressure gradients and changes in curvature.One more reason for interest in this flow is because of its similarity to the boundary layer over a wing, turbine blades, etc. High quality experimental data is available for this flow as well as information about the inlet conditions (Webster et. al, 1996).In flows with a single direction of inhomogeneity (such as a flat plate boundary layer) 22 R φ = and 3 1 2 R ψ = so the turbulent potentials can be prescribed accurately at the inlet.Wu & Squires (1998) were able to predict three-dimensional boundary layer over a swept bump using the RANS model developed by Durbin (1993b) and using large eddy simulations.The bump is defined by three circular arcs (Figure 5 and 6), and allows the examination of combined effects of the surface curvature and the streamwise pressure gradients on the mean flow.The boundary layer experiences several changes in the streamwise pressure gradient: firstly mildly adverse, then strongly favorable, strongly adverse and finally mildly favorable.
The flow is nearly homogeneous in the direction parallel to the bump because a suction slot is used.If the top boundary layer and the channel sidewall effects are neglected then a twodimensional domain (Figure 7) can be used to calculate the flow.This reference frame will be referred to as the computational frame of reference.It is at an angle of 45 to the original (experimental) frame of reference.An unstructured grid, shown in (Figure 8), with high near wall resolution is used to calculate the flow.    .Symbols are the model predictions and the lines are the experimental data.(a) Streamwise velocity (b) Spanwise velocity.Downstream of the bump the model seems to recover too quickly compared to the experiment.It is unclear how much of this recovery difficulty is due to the unknown upper boundary layer and sidewall effects, and how much is due to the model.This is the opposite of the recovery problem after a backward facing step where the model solution does not recover quickly enough.
Similar results were observed when a zero pressure gradient two-dimensional boundary layer with 3260 Re = θ was used as the upstream condition.The streamwise and spanwise velocity profiles scaled using the initial boundary layer edge velocity U 0 for an initial 3260 Re = θ , are shown as function of the distance from the channel floor y in figure 10.

UNSTEADY VORTEX SHEDDING
To test this model's ability to predict unsteady nonequilibrium turbulent flow, the problem of vortex shedding behind a 2D triangular cylinder was chosen.The flow is inherently unsteady.Complicated phenomena such as separation and large-scale vortices coexist with the turbulence.This geometry is slightly easier to simulate than the circular cylinder, since the separation points are fixed.The total mass flow was 1 m =0.6 kgs -1 in their experiment, and the inlet velocity is evaluated based on that value.These values are also used as the initial value for the whole domain.
is the height of the duct (which is 3 times of the height of the triangle).A zero gradient boundary condition is used for all the variables at the outlet.Slip-wall boundary conditions are used for the duct wall.The Reynolds number of this simulation is in Re=U d/ν=45,000 , where d is the side length of the triangle.
Unsteady behavior is due to vortices alternately shedding from the upper and lower edges of the cylinder, forming a Von Karmann vortex street behind the triangle.No special triggering measure is taken to start the vortex shedding, the unsteadiness in the computational results evolve naturally.It was triggered by the machine error and asymmetry of the mesh.
To illustrate the periodicity of the flow, the stream function of a point about one triangle height behind the triangle near the centerline is shown in Figure 12.It can be seen that an almost The corresponding Strouhal number defined by, Sr=fd/U in is 0.257, which should be compared with experimental data of 0.25 and the computed value of 0.27 in Johnnasson (1993).Figure 13 shows an instantaneous velocity vector plot, we can see that the center of a vortex is rolled up at the lower edge and a new vortex is beginning to roll up at the upper edge.The vortex street can also be seen in the instantaneous vorticity contours plot shown in Figure 14.Although the instantaneous flow is asymmetric, the timeaveraged fields are always symmetric or anti-symmetric.Figure 15 shows the mean stream-wise velocity at the centerline.The length of recirculation zone is accurately predicted, while the location of the maximum negative velocity is slightly upstream compared with the experiments.The magnitude of the maximum negative velocity is also a little lower than the experiment data.
Figure 16 shows the stream-wise velocity at different cross sections behind the triangle.The calculated velocity profiles are in reasonable agreement with the experimental data.However, it is hypothesized that the boundary layer on the triangle is not fully resolved due to mesh size restrictions.The computed boundary layer is thicker than the real one, thus close to the back of the triangle, the fluid is slowed down and driven backwards more than it should be.This would explain the mean velocity profile close to the centerline at x=15mm where the velocity is under-predicted.In the another simulation, we use half of the domain mentioned above and imposed a symmetric boundary condition along the centerline.Figure 17 shows the comparison of the mean stream-wise velocity contours.The top simulation is the time averaged unsteady streamwise velocity.The bottom simulation is the steady solution imposed by forcing a symmetry condition at the centerline.The contour levels in each plot are the same.The steady solution has a recirculation zone that is much longer than the time-averaged unsteady solution.The reason for this is that the unsteady flow increases the momentum exchange between the wake and it's surrounding, thus reducing the recirculation zone.In the steady simulation the turbulence model does a poor job representing the momentum exchange due to the large shedding vortices.In this work, a flow past a square cylinder at Re = 21,400 (which is based on upstream velocity and cylinder side height d) is also simulated.The geometry and mesh of this simulation is shown in figure 18.On the left boundary, the x-component velocity is set to a constant while the y-component is set to zero.On the right boundary, dynamic pressure is a constant and the gradient of y component velocity is zero.The top and bottom boundaries are slip walls.The boundary on the cylinder is a solid wall.The turbulence fluctuation velocity is about 1% of the mean velocity.The turbulence quantities evaluated based on this is taken to be the initial condition and is also fixed on the left boundary (taken as a boundary condition).Figure 19 shows the mean stream-wise velocity at the centerline.There are two experimental data sets available.Lyn's (1995) andDurao's (1988) data are quite different in the downstream region (x/d>4).When x/d < 4, our simulation matches both experiments well.When x/d >4, our simulation is closer to Durao's data.Similar LES results have been reported by Bouris and Bergeles(1996), Murakami and Mochida(1995).

CONCLUSIONS
The turbulent potential model demonstrates the ability to accurately predict the mean flow behavior of the suddenly spanwise driven channel flow.The effects observed in this flow cannot be reproduced by most RANS models.The turbulent potential model was able to reproduce the observed effects in the turbulence quantities by a small addition to the ε -equation.The change in the ε -equation is similar to the suggestion proposed by Durbin (1993a).The turbulence quantities show a reasonable agreement with the DNS data and the correct qualitative trends both spatially and temporally.The predictions of the turbulent potential model for the three-dimensional boundary layer over a swept bump are in good agreement with the experimental data for the mean streamwise and spanwise velocity.
Given the uncertainties of the experiment (i.e. the top boundary layer and sidewall effects), we can not expect the model predictions to compare any more favorably.
In predicting unsteady vortex shedding, the turbulent potential model shows good overall agreement with the experiments.The shedding frequency is with 5% and the mean flow cross sections and centerline predictions are within the experimental error.Our investigations showed that the model naturally models the turbulence and not the large-scale vortices.The shedding occurs spontaneously in the simulations.
This work has demonstrated the efficacy of the turbulent potential modeling approach for complex non-equilibrium flows.The model accurately predicted both three-dimensional boundary layers and unsteady vortex shedding.Its capability of capturing non-equilibrium turbulent flows and inexpensive computing cost make the turbulent potential model an attractive approach for engineering applications.

Figure 3 .
Figure 3. Spanwise mean velocity at times of 0.0, 0.3, 0.6 and 0.9.Symbols are DNS data of Moin et.al., Solid lines are the turbulent potential model predictions.

Figure 4 .
Figure 4. Turbulent shear stress profiles at 0.0,0.3,0.6 and 0.9.Symbols are DNS data of Moin et.al., Solid lines are the turbulent potential model predictions.The upper group of curves are v u ′ ′ , and the lower group of curves are w v ′ ′ .

Figure 5 .
Figure 5. Top view of the swept bump.

Figure 8 .
Figure 8. Mesh for the swept bump computational domain The streamwise velocity profiles for an initial Re θ =1400 are shown as function of the distance from the channel floor y in Figure 9a.U 0 is the initial boundary layer free-stream velocity and distances are given in meters .The results are shown at various positions downstream of leading edge of the bump.The model predictions are the symbols and the experimental results are the lines.The model predictions are in good agreement with the experimental results.The streamwise component of velocity increases and reaches a maximum value at the apex of the bump.On the trailing edge, the flow is very close to separation because of the adverse pressure gradient caused by the flow expansion.The flow relaxes to a two-dimensional boundary layer as it moves downstream of the bump.

Figure 9 .
Figure 9. Velocity profiles in the experimental frame of reference for inlet Re 1400 θ = .Symbols are the model predictions and the lines are the experimental data.(a) Streamwise velocity, (b) Spanwise velocity.

Figure 11 .
Figure 11.Computational domain and mesh for flow past triangular cylinder.The computational domain is shown in figure 11.It is the same as Sjunnesson's (1991) experiment.Approximately 25,000 triangles are used in this simulation.The inlet mean stream-wise velocity is set to a constant and the vertical velocity is set to zero.For turbulent kinetic energy and dissipation rate, we use the same boundary conditions described Johnasson (1993).

Figure 12 .
Figure 12.The stream-function of one point about one cylinder height behind the triangle near the centerline.

Figure 17 .
Figure 17.Mean stream-wise velocity contours of the timeaveraged unsteady and steady solution.

Figure 18 .
Figure 18.Computational domain and mesh for flow past square cylinder.

Figure 19 .
Figure 19.The mean stream-wise velocity at the centerline for a flow over square cylinder.
Mean stream-wise velocity behind the triangle: calculations; experiments.(a) x/d =1 (b) x/d=3 (c) x/d=6 downstream of the center of the square.

Figure 21
Figure 21 is an instantaneous streamwise velocity contours plot.As shown in this figure, separation occurs on the top and bottom wall of the square cylinder.The flow pattern of this geometry is more complicated than that of the triangular cylinder where separation occurs on the back corner.