Week 8: Literature review RANS derivation and analysis
Derivation of Reynold's stresses from the Reynold's Averaged Navier-Stokes equations
The incompressible Navier-Stokes equations are given as:
Continuity equation
`(partialu)/(partialx) + (partialv)/(partialy) + (partialw)/(partialz) = 0`
Momentum equation in the x-direction
`(partialu)/(partialt) + u(partialu)/(partialx) + v(partialu)/(partialy) + w(partialu)/(partialz) = -1/rho(partialP)/(partialx) + nu((partial^2u)/(partialx^2)+(partial^2u)/(partialy^2)+(partial^2u)/(partialz^2))`
where,
u, v, and w are the velocity components,
x, y, and z are the spatial directions,
`rho` is the density of the fluid,
P is the static pressure,
`nu` is the kinematic viscosity
Reynolds decomposition
Flow variable (`phi`) can be decomposed into a time-averaged component (`bar phi`) and a fluctuating component (`phi'`).
`phi (x,t) = bar phi (x) + phi' (x,t)`
where x is the spatial position and t is the time.
Time-averaged variable `bar phi ` is defined as
`bar phi = lim_(T->infty) (1/T) int_0^T phi dT`
Properties of a time-averaged variable
- Time-averaging of a fluctuating component is zero.
`bar (phi ') = 0` - Time-averaging does not affect an already time-averaged variable.
`bar(bar phi) = bar phi` - Differential operators commute with time averaging.
`(bar(partial phi))/(partial x) = (partial bar phi)/(partial x)` - Averaging a constant `lambda` times a variable `phi`:
`bar(lambda phi) = lambda bar phi`
Therefore, if a and b are field variables,
`bar (bar a b) = bar a* bar b`
However,
`bar (a*b) ne bar a *bar b`
Using these principles, Reynold's decomposition of velocity variables u, v, and w will be:
`u = bar u + u'`
`v = bar v + v'`
`w = bar w + w'`
Plugging the Reynold's decomposed velocity variables in the continuity equation gives:
`(partialbaru)/(partialx) +(partialu')/(partialx) + (partialbarv)/(partialy) + (partialv')/(partialy) + (partialbarw)/(partialz) + (partialw')/(partialz)= 0` `rightarrow ` Eq 1
Applying time-averaging on the above equation gives:
`(partialbaru)/(partialx) + (partialbarv)/(partialy) + (partialbarw)/(partialz)+ (partialbar(u'))/(partialx) + (partialbar(v'))/(partialy) + (partialbar(w'))/(partialz)= 0`
We know that, `bar(u') = 0`
this implies `implies` `(partial bar(u'))/(partialx) = 0 `
likewise, `(partial bar(v'))/(partialy) = 0 ` and `(partial bar(w'))/(partialz) = 0 `
The above equation then reduces to:
`(partialbaru)/(partialx) + (partialbarv)/(partialy) + (partialbarw)/(partialz) = 0` `rightarrow ` Eq 2
Another important result is Eq2 - Eq1:
`(partialu')/(partialx) + (partialv')/(partialy) + (partialw')/(partialz)= 0` `rightarrow ` Eq 3
Plugging the Reynold's decomposed velocity variables in the momentum equation in the X-direction gives:
`(partialbaru)/(partialt) + (partialu')/(partialt) + [(baru*(partialbaru)/(partialx) + baru*(partialu')/(partialx) + u'*(partialbaru)/(partialx) + u'*(partialu')/(partialx)),( + barv*(partialbaru)/(partialy) + barv*(partialu')/(partialy) + v'*(partialbaru)/(partialy) + v'*(partialu')/(partialy)), (+ barw*(partialbaru)/(partialz) + barw*(partialu')/(partialz) + w'*(partialbaru)/(partialz) + w'*(partialu')/(partialz) ) ] = -1/rho ((partialbarP)/(partialx) + (partialP')/(partialx)) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2))`
`\ \+ nu((partial^2u')/(partialx^2)+(partial^2u')/(partialy^2)+(partial^2u')/(partialz^2))`
Applying time-averaging on the above equation:
`(partialbaru)/(partialt) + (partialbar(u'))/(partialt) + [(baru*(partialbaru)/(partialx) + baru*(partialbar(u'))/(partialx) + bar(u')*(partialbaru)/(partialx) + bar(u'*(partialu')/(partialx))),( + barv*(partialbaru)/(partialy) + barv*(partialbar(u'))/(partialy) + bar(v')*(partialbaru)/(partialy) + bar(v'*(partialu')/(partialy))), (+ barw*(partialbaru)/(partialz) + barw*(partialbar(u'))/(partialz) + bar(w')*(partialbaru)/(partialz) + bar(w'*(partialu')/(partialz)) ) ] = -1/rho ((partialbarP)/(partialx) + (partialbar(P'))/(partialx)) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2))`
`\ \+ nu((partial^2bar(u'))/(partialx^2)+(partial^2bar(u'))/(partialy^2)+(partial^2bar(u'))/(partialz^2))`
Since `bar(phi') = 0`, multiple terms in the above equation can be eliminated.
And `(partial baru)/(partial t) = 0 `, as `bar u` is constant with respect to time.
The above equation now becomes:
`cancel((partialbaru)/(partialt)) + cancel((partialbar(u'))/(partialt)) + [(baru*(partialbaru)/(partialx) + cancel(baru*(partialbar(u'))/(partialx)) + cancel(bar(u')*(partialbaru)/(partialx)) + bar(u'*(partialu')/(partialx))),( + barv*(partialbaru)/(partialy) + cancel(barv*(partialbar(u'))/(partialy)) + cancel(bar(v')*(partialbaru)/(partialy)) + bar(v'*(partialu')/(partialy))), (+ barw*(partialbaru)/(partialz) + cancel(barw*(partialbar(u'))/(partialz)) + cancel(bar(w')*(partialbaru)/(partialz)) + bar(w'*(partialu')/(partialz)) ) ] = -1/rho ((partialbarP)/(partialx) + cancel((partialbar(P'))/(partialx))) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2))`
`\ \+ nu(cancel((partial^2bar(u'))/(partialx^2)) + cancel((partial^2bar(u'))/(partialy^2)) + cancel((partial^2bar(u'))/(partialz^2)))`
After elimination:
`[(baru*(partialbaru)/(partialx)+ barv*(partialbaru)/(partialy) +barw*(partialbaru)/(partialz)), (+ bar(u'*(partialu')/(partialx)) + bar(v'*(partialu')/(partialy))+ bar(w'*(partialu')/(partialz)))] =-1/rho (partialbarP)/(partialx) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2)) ` `rightarrow ` Eq 4
Considering the fluctuating components in the above equation:
`implies` ` bar(u'*(partialu')/(partialx)) + bar(v'*(partialu')/(partialy))+ bar(w'*(partialu')/(partialz))`
Using the product chain rule of differentiation to the above, we get:
`bar(u'*(partialu')/(partialx)) + bar(v'*(partialu')/(partialy))+ bar(w'*(partialu')/(partialz)) = (bar(partialu'u'))/(partialx) - (bar(u'partialu'))/(partialx) + (bar(partialv'u'))/(partialy) - (bar(u'partialv'))/(partialy) + (bar(partialw'u'))/(partialz) - (bar(u'partialw'))/(partialz)`
`\ \ = (bar(partialu'u'))/(partialx) +(bar(partialv'u'))/(partialy) + (bar(partialw'u'))/(partialz) - bar((u'partialu')/(partialx) - (u'partialv')/(partialy) - (u'partialw')/(partialz))`
Since, `bar(phi_1 + phi_2) = bar phi_1 + bar phi_2`
`implies``bar(u'*(partialu')/(partialx)) + bar(v'*(partialu')/(partialy))+ bar(w'*(partialu')/(partialz)) = (bar(partialu'u'))/(partialx) +(bar(partialv'u'))/(partialy) + (bar(partialw'u'))/(partialz) - bar(u'((partialu')/(partialx)+(partialv')/(partialy)+(partialw')/(partialz)))`
From Equation 3 we know that,
`(partialu')/(partialx) + (partialv')/(partialy) + (partialw')/(partialz)= 0`
`implies` `bar(u'*(partialu')/(partialx)) + bar(v'*(partialu')/(partialy))+ bar(w'*(partialu')/(partialz)) = (bar(partialu'u'))/(partialx) +(bar(partialv'u'))/(partialy) + (bar(partialw'u'))/(partialz)`
Plugging above into Equation 4 gives:
`[(baru*(partialbaru)/(partialx)+ barv*(partialbaru)/(partialy) +barw*(partialbaru)/(partialz)), (+ (bar(partialu'u'))/(partialx) +(bar(partialv'u'))/(partialy) + (bar(partialw'u'))/(partialz))] =-1/rho (partialbarP)/(partialx) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2)) `
The X-component of the Momentum equation can now be written as
`baru*(partialbaru)/(partialx)+ barv*(partialbaru)/(partialy) +barw*(partialbaru)/(partialz) =-1/rho (partialbarP)/(partialx) -((bar(partialu'u'))/(partialx) +(bar(partialv'u'))/(partialy) + (bar(partialw'u'))/(partialz)) + nu((partial^2baru)/(partialx^2)+(partial^2baru)/(partialy^2)+(partial^2baru)/(partialz^2)) `
The fluctuating components in the above equation are:
`(-(bar(partialu'u'))/(partialx) -(bar(partialv'u'))/(partialy) - (bar(partialw'u'))/(partialz))`
Multiplying and dividing the respective terms by density `rho` :
`implies` `1/rho*(-(partial bar(rhou'u'))/(partialx) - (partial bar(rho v'u'))/(partialy) - (partial bar(rho w'u'))/(partialz))`
The terms `-bar(rhou'u')`, `-bar(rho v'u')`, `-bar(rho w'u')` are the X-component of the Reynold's stresses. These terms are not inherently related to stress but the units of these terms resemble the units of stress.
In three dimensions, the Reynolds stress are:
`[(-bar(rhou'u'), -bar(rhov'u'), -bar(rhow'u')), (-bar(rhou'v'), -bar(rho v'v'), -bar(rhow'v')), (-bar(rhou'w'), -bar(rho v'w'), -bar(rhow'w'))]`
The above matrix with the 3D Reynold's stress is called the Reynold's stress tensor.
The fluctuating Reynolds stress terms are undesirable for numerical analysis and hence need to be modeled as a function of the mean flow. This closing of the fluctuating terms is called the closure problem. One such method of closing the problem is by using Eddy viscosity models.
The turbulent eddy viscosity is a constant of proportionality that relates the Reynolds stress terms to their respective velocity gradients using the Boussinesq hypothesis:
`-bar(rho(u_i'u_j')) = nu_t ((partialu_i)/(partialx_j)+(partialu_j)/(partialx_i)) -2/3 k delta_(ij)`
where,
i,j corresponds to x-direction if =1, corresponds to y direction if =2, corresponds to z direction if =3,
`nu_t ` is the turbulent eddy viscosity,
k is the turbulent kinetic energy = `1/2 ( bar(u^2') + bar(v^2')+ bar(w^2'))`,
`delta_(ij)` is the kroenecker delta
Molecular viscosity is different from turbulent eddy viscosity such that molecular viscosity is the measure of a fluid's internal resistance to flow under the force of gravity.