Transition to Turbulence

Created by Ian Nesbitt and Jukes Liu on 2019-02-19

<Insert fun image of model here>

Overview

In this module, we will examine the transition to turbulence and learn how viscosity affects the dynamics of a simple fluid system with a single cylindrical impediment placed in the flow path. Our tools in this endeavor are COMSOL Multiphysics, SPH, and ParaView. At the end of the module, you should be able to describe the initial conditions of a system, the transition from laminar to turbulent flow, and the effect viscosity has on the inertia, complexity, and predictability of a system!

We will move from looking at models with higher viscosity regimes and move to lower viscosity. Viscosity is the measure of a fluid's resistance to deformation.

First, let's start by learning how to talk about the various components of flow.

Characterizing the fluid dynamics

When we talk about fluid flow, we need to first describe the fluid. Some of the fundamental characteristics of fluid are density, viscosity, compressibility, heat capacity, and thermal conductivity. In this exercise we will ignore the thermal properties entirely, define a constant density and compressibility, and focus most of our attention on how viscosity affects the kinematics of a fluid system.

Compressibility and density are fairly easy to define. Compressibility is the tendency of a fluid to change its volume due to changes in pressure. Since we are not dealing with immense pressure, we assume that compressibility is zero across all models. Density is the quantity of mass per unit volume. For the purposes of this model, we will hold density constant at 1 kg m^-3 to try and isolate the effects of viscosity.

These parameters will be fed into the Navier-Stokes fluid momentum equation (below) for each location in the model at each time step.

Navier-Stokes momentum equation

To characterize the motion of particles in the fluid, we must account for the changes in momentum. This Navier-Stokes equation is the foundation of our models. It is derived based on the conservation of momentum in a system. This equation includes the influence from inertial forces, body forces, viscous forces, and the pressure gradient throughout the system.

Let's take a look at the structure of the Navier-Stokes equation and what each term means before diving into the math.

Acceleration = Body forces + Shear forces - Pressure gradient

These terms are represented by some greek symbols which may or may not mean anything to you right now, but we'll break them down shortly.

∂v / ∂t = ρg + µ∇²U - ∇P

where
∂v	= change in velocity
∂t	= change in time
ρ	= density
g	= acceleration due to gravity
µ	= viscosity
∇2U	= divergence of velocity (gradient)
∇P	= divergence of pressure (gradient)

The inertial forces (acceleration) is represented by ∂v / ∂t. On the right hand side of the equation, we have the body forces (ρg), viscous forces (µ∇²U), and the pressure gradient (∇P). Solving this equation for each grid point or particle at each time step will define how our fluid will behave in a computational sense. Once our fluid starts to move, we need ways to describe how it is flowing. One of those ways is whether or not the flow is turbulent. The parameter that is used to describe this characteristic is called the Reynolds number.

Reynolds number

The flow regime of a fluid is defined by the dimensionless Reynolds number, Re. Reynolds number is a way of quantifying whether a flow is laminar, transitional, or fully turbulent. Re can also be thought of as inversely related to how well flow conditions can be described mathematically across a given dimension (x, y, z, or time). Reynolds number is equal to the inertial forces (velocity times characteristic length) over the resisting forces (viscosity):

Re = ρVD/µ = VD/v

where
V	= fluid velocity (m s-1)
D	= linear dimension (m)
µ	= dynamic viscosity (Pa s)
v	= kinematic viscosity (Pa s)
ρ	= density (kg/m^3)

The Reynolds number is a function of viscosity and velocity: f(µ, v). Velocity and viscosity have inverse influences on the Reynolds number (i.e. as velocity increases, the Reynolds number increases while increasing viscosity lowers the Reynolds number). We are choosing to focus on examining the effects of lowering viscosity by examining fluid flow fluids with different viscosities (honey, water, air) while holding the velocity constant. This would have a similar effect as increasing velocity and holding viscosity constant. Using our models, we examine the effects of changing viscosity, effectively changing the Reynolds number, on fluid flow.

The table below lists the qualitative descriptions for each quantitative Re range.

Quantifying qualitative descriptions of Reynolds number
Laminar	Transitional	Turbulent
Re < 100	100 < Re < 10^5	10^5 < Re

Flows with Reynolds numbers of less than 100 are considered "laminar," while those from 100 to 10^5 are "transitional," and above 10^5 are "turbulent." At low Reynolds numbers (Re < 100), viscous forces dominate and flow is laminar. As the Reynolds number increases, due to lower velocity or higher viscosity, inertial forces begin to play a greater role and flow becomes increasingly turbulent. In this transition between laminar and turbulent flow, we begin to observe some unique fluid phenomena. Let's describe our model and then take a look at these phenomena!

Model setup

Geometry

Below is the geometry of our simple model in COMSOL.

This is a two-dimensional COMSOL model, which assumes that the channel is of infinite depth and thus does not produce drag from the bottom of the flow in any way. Both the horizontal (long) walls and the walls of the cylinder have a "no slip" condition. This means that as the flow gets closer to the wall, velocity goes to zero. This condition sets up a velocity gradient between the walls at zero and the center of the flow. The wall on the far left side has an inflow condition of 1 meter per second, and the wall on the right side has a zero-pressure outlet condition, which lets flow pass out of the model unhindered.

When we run the model, we are executing a variation of the Navier-Stokes momentum equation for each mesh element of the model at each time step. This equation, shown above, is fundamental to describing the behavior of fluids in continuum mechanics.

Experimental variable: viscosity

Let's look at the dynamic viscosities of fluids we'll be working with.

Fluid	µ (viscosity at 20 ºC)
1. Honey	1 * 101.5 Pa s
2. Water	0.89 * 10-3 Pa s
3. Air	1.81 * 10-5 Pa s

As mentioned above, we use an initial fluid velocity of 1 m/s in all of our models.

Questions to consider

From the diagram of the model geometry, what is the linear dimension (D)? In this model, it is equivalent to the diameter of the cylinder.

Before we go through the model results, try calculating the initial Reynolds numbers for each fluid viscosity using the linear dimension and the fluid velocity. Based on the Reynolds number calculated for the honey, do you think the fluid flow will be more laminar or turbulent? What about for water? Air?

Models

Honey (μ = 10^1.5 Pa s)

In the first example, we have a fluid with the viscosity of about that of honey, μ = 10^1.5 Pascal seconds (Pa s). You might imagine honey to be a pretty "stiff" liquid, which has something to do with its viscosity. What do you think will happen? Let's look at a few parameters of flow in the honey model.

Velocity 1 (x component)

Things are pretty orderly and logical here. Not a whole lot changes over time. You might describe this flow system as "laminar," since you can describe what will happen at any given point in the model space at any given time fairly easily. Note that velocity is highest near the cylinder where the honey must "squeeze by" the cylinder in order to keep up with the 1 m s^-1 inflow rate. This sets up a steep velocity gradient near the no-slip boundary of the cylinder, which we can see in the shear rate:

Shear rate 1

Shear rate is by definition the derivative of the velocity field, so it is no surprise that the highest shear values are in the areas next to the cylinder where the steep velocity gradient exists. If we break shear down into its component parts, we get vorticity. Vorticity is the curl of the velocity field, or in other words, the rate of rotation of the angle between two particles as they slip past each other.

Vorticity 1

Again, unsurprisingly, the curl is highest where the gradient of velocity is high. This makes sense because as two particles slip past each other at different velocities, their orientation changes (rotates) with regards to one another. Vorticity will be highest near the no-slip boundaries in a simple laminar system like this.

Pressure 1

The pressure field should also make sense, because our one obstacle in this model sets up the overall gradient from inflow to outflow pretty nicely. Notice that the pressure values are very slightly higher on the stoss side and slightly lower on the lee side of the cylinder. Note also that pressure gradients define the shortening and stretching of a system. Flow is shortened as it is forced towards the front of the cylinder and stretched as it goes past the cylinder.

Reynolds number 1

Reynolds number describes how chaotic a system is. Since our system has almost no chaos, this plot is pretty boring at this viscosity.

Explain difference between cell Reynolds number and regular Reynolds number.

Discussion 1

"The high viscosity of honey results in perfectly laminar flow when poured from a bucket, while the low surface tension allows it to remain sheet-like even after reaching the fluid below. ... when the flow meets resistance it slows and [shortens by folding] upon itself."[1] (used under Creative Commons 2.0)

Okay, so what's happening here? Things look pretty "boring" so far. With the exception of the very first couple of frames, in which the fluid has not achieved its full velocity yet, all of the frames are visually nearly the same. As you can tell by the color scale on the velocity plot, the fluid is moving as fast as 2 m s^-1 past the cylinder, but there seems to be no instability here. This system is the result of kinematic forces at work, but the system seems to be stable, i.e. there is very little dynamic element to this system. Another way to say this is that there is not much variance in any given cell through time.

Why is this? Well, we said before that honey is "sticky." In more technical terminology, higher viscosity means higher internal friction and resistance to deformation. This has the effect of dampening the inertial term in Navier-Stokes, and thus causes the flow to gain stability and to behave in a "boring," predictable manner. Take another look at the Reynolds number plot. All areas on the plot are Re < 1. This indicates a very stable and predictable flow, in which the kinematics of any one particle anywhere in the continuum can be described quite well mathematically. Let's take another look at the table we saw above:

Quantifying qualitative descriptions of Reynolds number
Laminar	Transitional	Turbulent
Re < 100	100 < Re < 10^5	10^5 < Re

Flows with Reynolds numbers of less than 100 are considered "laminar," while those from 100 to 10^5 are "transitional," and above 10^5 "turbulent."

Water (μ = 10^-3 Pa s )

Now let's look at how flow changes when viscosity is lowered more than four orders of magnitude to that of water, μ = 10^-3 Pa s (remember, we're holding density constant). What do you think will be different? What do you think the system will look like?

Velocity 2 (x component)

Whoa! That's quite different. Notice that things initially begin laminar but start to decay into a regular pattern of "wobbling." This is called a von Kármán vortex street. Notice also how violently the velocity vectors swing around as the vortices begin forming and shedding and the flow field oscillates from side to side. Flow is actually going backwards in some places! Why?

There's no question this flow is more dynamic than the first one. Technically, as you'll see just below, this is still a laminar flow. Why? Can you still predict what's happening in certain places or at certain times? Why or why not? Keep these questions in mind as you scroll through the next few plots.

Shear rate 2

Here we see the shear plot, which still tracks the derivative of the velocity field, but in this case it's much more dynamic as well. We can see the formation of vortices as the dominant shear jumps back and forth between the two sides of the cylinder. Why do you think this is happening?

Vorticity magnitude 2

What a gorgeous vorticity plot. After you've admired it for a bit, can you describe what's happening here? Look back at the shear plot. When the newest vortex is about to detach from the lee side of the cylinder, where is the shear rate highest?

Pressure 2

Look at the similarities between this pressure plot and the vorticity plot, and you'll notice that pressure is inversely related to vorticity. Why? Can you think of real-life scenarios (other than flow past a cylinder in a flume) where you can observe this phenomenon? Recall the flow vectors going backwards. Given this pressure plot, can you give an explanation for this phenomenon?

Reynolds number 2

As discussed above, the Reynolds number doesn't quite make it above the critical number of 10 in order to be considered a "transitional" flow. Why do you think this is? What is happening in this flow that keeps it relatively "predictable?"

Note: The stippling is an artifact of the meshing process in COMSOL.

Discussion 2

Quasi-laminar flow past a cylinder (a GPS survey staff) on Third River in Montclair, NJ. (Photo: Ian Nesbitt, licensed under Creative Commons 4.0)

This specific type of oscillating flow is called a von Karman vortex street, after Theodore von Karman, the Hungarian-American physicist who first described it in detail. The process that leads to this type of flow is called vortex shedding.

The vortex shedding pattern shown in these examples develops after a couple of brief moments in which the flow is fairly laminar and there are few dynamics acting on the system. Then the system begins to develop an instability. In real life this instability is the result of tiny perturbations on the surface of the cylinder that cause differences in the boundary layer thickness as flow passes, causing one side to briefly slip faster than the other. The momentum of the vortex continues over to the other side of the cylinder, which increases the fluid pressure on that side and thus slows flow. As pressure always wants to equilibrate, this process reverses and a competing vortex is sent spinning to the opposite side. If flow conditions do not change, this process will continue indefinitely.

In mathematical terms, fluid particles will accelerate over time—∂v / ∂t—towards areas of low pressure—- ∇P (or decelerate if they are already traveling away from an area of low pressure). In order to figure out whether a particular particle will accelerate (or decelerate) towards a certain area of pressure, you'd have to examine the direction magnitude component partial differentials of the velocity and term, ∇²U. In other words, the grad (∇) symbol is a shorthand for the fact that you're looking at something that changes in more than one linear dimension (in the case of this simple 2D model, ∇ denotes that the variable will change in both x and y).

In a 3D model, things progress in much the same way. Let's look at a particle model of the same process. Particles in the first plot here are given colors based on the x component of their velocity.

In the second plot, particles are given colors based on their vorticity.

As a fluid becomes less viscous, the momentum term will begin to dominate the Navier-Stokes solution, meaning that the internal shear strength of the fluid is less able to keep it from deforming. As we said before, in this flow past a cylinder model, that means vortex shedding will become more and more viable at either lower viscosities or higher velocities.

Let's take a look at this behavior in a fluid with two orders of magnitude less viscosity than water: air.

Air (μ = 10^-5 Pa s)

Air has a viscosity of μ = 10^-5 Pa s (still holding density and velocity constant). How do you think this will affect the deformation of the fluid?

Velocity (x component) 3

Shear rate 3

Vorticity magnitude 3

Discuss dissipation of turbulence.

Pressure 3

Reynolds number 3

The stippling is an artifact of the meshing process in COMSOL.

Complex flow

Although our models don't show it here, we think it would be good to show some examples of more complex flows, which make up the vast majority of flow behavior you'll encounter in the real world. It's important to note that the models above are oversimplifications meant to highlight the physical forces acting on fluids and making them behave in certain ways. They aren't wrong per se, but they are simplistic.

Fully turbulent flows

A hydraulic jump causing a transition from laminar to turbulent flow at a US Army Corps of Engineers flood control project on the Hoosac River in North Adams, Massachusetts, four years to the day after peak flow recorded during Hurricane Irene (see video, left). Photo is taken near USGS gaging station 01332500 (green shed with solar panel). (Photo: Ian Nesbitt, August 28, 2015, licensed under Creative Commons 4.0)

Fully turbulent flow, where 2000 < Re and predictability over a scale larger than a couple of particles or beyond a couple of time steps approaches zero, is very difficult to model. In some cases, flows of this nature are easier to compute and the solutions are better represented by Smoothed Particle Hydrodynamic models, in which individual particle smoothing kernels, rather than mesh elements, define the behavior of the fluid. However, flows in high Reynolds cases are are difficult to describe mathematically, and thus solutions that approach full turbulence in continuum or smoothed particle models can take orders of magnitude longer to compute.

The gif above was filmed on August 28, 2011, near peak flow during Hurricane Irene, the current flood of record for USGS gaging station 01332500 near Williamstown, MA. The video shows flow over a US Army Corps of Engineers (USACE) flood control structure (also pictured at baseflow on right). The concrete-bedded flood control flume has low bed roughness and thus high velocity, and comparatively lower Re. Flows like this can be very erosive when they encounter a jump in bed roughness due to a transition from a flood-control flume with a concrete bed to natural (rougher) bed conditions. The purpose of induced turbulence—from USACE's perspective—is to limit scour caused by a sudden increase in bed roughness. By reducing friction at the boundary layer, the flow will not be quite as erosive in its transition between bed conditions. Here, however, it is near its design capacity to limit the destructive nature of the flow, and the flow will probably still erode smaller particles from the natural bed.

"Multiphase" flows

Using a particle physics solution, one can model fluid particle behavior without the need for computationally remeshing, which can have substantial processing power overhead in a continuum model. In this case we use the word "multiphase" to describe the use of two fluids, each with different properties, to initialize the model. The following GIF shows a DualSPHysics two-phase example model of a dam break, where a parcel of water flows over and deforms denser substrate (sediment) material. This type of model can be used to examine physical processes of erosion at a small scale, and potentially could be used to reproduce and expound on results of physical modeling studies of these phenomena.

As you can see here, water flowing on the surface of the sediment creates shear stress at the bed, which both deforms the sediment to certain depths, and entrains grains from the surface. You'll notice that as the front of the flow takes on sediment it behaves more like a debris flow than a water wave. This has the effect of buoying dense sediment up to the surface, which then remains entrained as the flow hits the wall and begins to curl around. From the DualSPHysics documentation:

These multi-phase sediment scouring phenomena are induced by rapid flows creating shear forces at the surface of the sediment which causes the surface to yield and produce a shear layer of sediment suspended particles at the interface and finally sediment suspension in the fluid. Applications include scouring in industrial tanks, port hydrodynamics, wave breaking in coastal applications and scour around structures in civil and environmental engineering flows among others.

Some of the more relevant physical properties of this model are shown in the following table.

Fluid properties of multiphase model
Fluid	Density (kg m-3)	Viscosity (Pa s)	Speed of sound (m s-1)	Cohesion coeff.	Angle of internal friction (deg)
Water	1000	0.001	80	0	0
Saturated sediment	1500	0.002	81	1	35

A physical multiphase model typically involves a flume, a pump, and a specific particle size or range of sizes to test the effects that flow-induced boundary stress will have on particle saltation, sliding, and entrainment. A simple sediment flume demonstration conducted by Dr. Ronadh Cox (Williams College) is shown below. Apart from the geometry, one major difference between this physical model and the numerical one is that the grain density is that of quartz (2650 kg m^-3) in the physical model.