Week 8: Literature review - RANS derivation and analysis
AIM: To demonstrate an expression for Reynolds stress by applying Reynolds decomposition to the Navier Stokes equation.
INTRODUCTION: In this project, I will demonstrate a detailed explanation of substituting Reynolds decomposition equations in the Navier Stokes equation.
THEORY:
Navier Stokes Equation: The Navier Stokes equation govern the motion of fluids and can be seen as Newton's second law of motion for fluids. In the case of a compressible Newtonian fluid, this yields
`rho((delu)/(delt)+u.gradu) = -gradp+grad.(mu(gradu+gradu)^T)-2/3mu(grad.u)I)+F`
where u is the fluid velocity, p is the fluid pressure, `rho` is the fluid density, and `mu` is the fluid dynamic viscosity. The different terms correspond to the inertial forces (1), pressure forces (2), viscous forces (3), and the external forces applied to the fluid (4). The Navier-Stokes equations were derived by Navier Poisson, Saint-Venant, and Stokes between 1827 and 1845.
These equations are always solved together with the continuity equations. The Navier-Stokes equations represent the conservation of momentum, while the continuity equation represents the conservation of mass.
How Do They Apply to Simulation and Modeling?
These equations are at the heart of fluid flow modeling. Solving them, for a particular set of boundary conditions (such as inlets, outlets, and walls), predicts the fluid velocity and its pressure in a given geometry. Because of their complexity, these equations only admit a limited number of analytical solutions. It is relatively easy, for instance, to solve these equations for a flow between two parallel plates or for the flow in a circular pipe. For more complex geometries, however, the equations need to be solved.
Governing Equations For Incompressible Flow:
Continuity Equation - Conservation of Mass:
`(delu)/(delx)+(delv)/(dely)=0`
Momentum Equation - Conservation of Momentum (Newton's Second Law):
`rho[(delu)/(delt)+u(delu)/(delx)+v(delu)/(dely)]= -(delP)/(delx)+mu((del^2u)/(dely^2))`
Energy Equation - Conservation of Energy (First Lae of Thermodynamics):
`rhoC_p[(delphi)/(delt)+u(delphi)/(delx)+v(delphi)/(dely)]=K[(del^2phi)/(delx^2)+(del^2phi)/(dely^2)]+mu((delu)/(dely))^2`
Reynolds-Averaged Navier-Stokes Equations: The RANS equations are time-averaged equations of motion for fluid flow. The idea behind the equations is Reynolds Decompositions, whereby an instantaneous quantity is decomposed into its time-averaged and fluctuating quantities, an idea first proposed by Osborne Reynolds. The RANS equations are primarily used to describe turbulent flows. These equations can be used with approximate time-averaged solutions to the Navier-Stokes equations.
Need of RANS: Many engineering simulations need to be computationally effective (less computational effort). To capture the minute effects of turbulence DNS cannot be applied effectively(high computation effort required) and it requires the mesh to be of order (1e-09 or less) to be in compliance with CFL number and small-time steps, which is going to be computationally inviable. Hence, to be computationally viable, engineers depend upon the time-averaged solutions, where the turbulence of engineering needs is captured effectively. Here the time step used is large compared to the relevant period of fluctuations of components. Hence a coarse grid than the one used for DNS can be used which will reduce the computational effort.
Reynolds Decomposition:
Through the Reynolds equation, we believe that the fluid is in a random unsteady turbulent state and work with the time-averaged or mean equations of motion. Any variable 'U' is resolved into a mean value of `barU` plus a fluctuating value `U'`, defined as.
`barU = 1/T_(i)int_0^(t_0+T_i)Udt`
Instantaneous quantity at any point in space and time =Time-averaged quantity +Fluctuating component.
Where T is large compared to the relevant period of fluctuations. The time gradient is fast when compared to the mean value of `barU`.
Let us consider the incompressible turbulent flow with constant transport properties
`u(x,y,t)=baru(x,y)+u'(x,y,t)`
`v(x,y,t)=barv(x,y)+v'(x,y,t)`
`w(x,y,t)=barw(x,y)+w'(x,y,t)`
`phi(x,y,t)=barphi(x,y)+phi'(x,y,t)`
`barU = 1/T_(i)int_0^(t_0+T_i)Udt`
mean values are not the function of 't'.
Moreover, the averaging of fluctuating terms must be equal to zero.
`=>u(x,y,z,t)= baru(x,y)+u'(x,y,z,t)`
`=>1/Tint _0^Tu(x,y,z,t)dt=1/Tint_0^Tbaru(x,y)=1/Tintu'(x,y,z,t)`
`=>cancel(baru(x,y,z))^0=cancel(baru(x,y,z))^0=1/Tint_0^Tu'(x,y,z,t)`
`1/Tint_0^Tu'(x,y,z,t)=0`
Substituting the time-averaged equation into the continuity equation.
continuity Equation:
`=> (delu)/(delx)=(delv)/(dely)=0`
`=> (del(baru+u'))/(delx)+(del(barv+v'))/(dely)=o`
`=> (delbaru)/(delx)+(delu')/(delx)+(delbarv)/(dely)+(delv')/(dely)=0`
averaging with respect to time 'T'
`=>1/Tint_0^T(delbaru)/(delx)dt+1/Tint_0^T(delu')/(delx)dt+1/Tint_0^T(delbarv)/(dely)dt+1/Tint_0^T(delv')/(dely)dt`
The time integral of the fluctuation component is zero and integral limit not a function of space
`=>1/Tint_0^T(delbaru)/(delx)dt+cancel(1/Tint_0^T(delu')/(delx)dt)^0+1/Tint_0^T(delbarv)/(dely)dt+cancel(1/Tint_0^T(delv')/(dely)dt)^0 =0`
by further simplifying we get
`=>1/Tint_0^T(delbaru)/(delx)dt+1/Tint_0^T(delbarv)/(delx)dt=0`
Since the `baru` and `barv` is the function of space only, the above equation can be re-written as:
`=>(delbaru)/(delx)1/Tint_0^Tdt+(delbarv)/(delx)1/Tint_0^Tdt =0`
`=>1/Tint_0^Tdt=(T-0)/T=1`
The continuity equation after the substitution of averaged quantities depends on only the mean values of decomposed quantities as shown below.
`=>(delbaru)/(delx)+(delbarv)/(dely)=0`
Momentum Equation:
`=> rho[(delu)/(delt)+u(delu)/(delx)+v(delu)/(dely)]= - (delP)/(delx)+mu((del^2u)/(dely^2))`
`=> mu/rho =v`
`=> (delu)/(delt)+u(delu)/(delx)+v(delu)/(dely)=-1/rho(delP)/(delx)+v(del^2u)/(dely^2)`
Multiply the above equation with x-component of velocity with continuity equation and add it to the LHS and its results will be zero.
`=> (delu)/(delt)+u(delu)/(delx)+v(delu)/(dely)+u((delu)/(delx)+ (delv)/(dely))=-1/rho(delP)/(delx)+v(del^2u)/(dely^2)`
`(delu)/(delt)+2u(delu)/(delx)+v(delu)/(dely)+u(delv)/(dely)= -1/rho (delP)/(delx)+v(del^2u)/(dely^2)` --------(1)
From product rule.
`=> (del(uv))/(dely)=u(delu)/(dely)+v(delu)/(dely)` -----------------------(2)
and
`=> 2u (delu)/(delx)=(delu^2)/(delx)` ---------------------------(3)
Now substitute the equation 2 and 3 in 1
`=> (delu)/(delt)+(delu^2)/(delx)+(del(uv))/(dely)=-1/rho (delP)/(delx) +v(del^2u)/(dely^2)`
Now substitute the averaged quantities
`=> (del(baru+u'))/(delt)+(del(baru+u')^2)/(delt)+(del(baru+u')(barv+v'))/(delt)= -1/rho (del(barP+P'))/(delx)+v(del^2(baru+u'))/(dely^2)`
`=> (delbaru)/(delt)+(delu')/(delt)+(delbaru^2)/(delx)+(del(u')^2)/(delx)+(del(2baru u'))/(delx)+(del(baru barv))/(dely)+(del(baru v'))/(dely)+(del(u' barv))/(dely)+(del(u' v'))/(dely)= -1/rho (delbarP)/(delx)-1/rho (delP')/(delx)+v(del^2abaru)/(dely^2)+v(del^2u')/(dely^2)`
averaging using the time step
`=> 1/Tint_0^T (delbaru)/(delt)dt+1/Tint_0^T (delu')/(delt)dt+1/Tint_0^T (delbaru^2)/(delx)dt+1/Tint_0^T (del(u')^2)/(delx)dt+1/Tint_0^T (del2(baru u'))/(delx)dt+1/Tint_0^T (delbaru barv)/(delx)dt+1/Tint_0^T (delbaru v')/(delx)dt+1/Tint_0^T (delu' barv)/(delx)dt+1/Tint_0^T (delu'v')/(delx)dt=-1/Tint_0^T 1/rho (delbarP)/(delx)dt-1/Tint_0^T 1/rho (delP')/(delx)dt+1/Tint_0^T v(del^2baru)/(dely^2)dt+1/Tint_0^T v(del^2u')/(dely^2)dt`
`=> 1/Tint_0^T (delbaru)/(delt)dt+cancel(1/Tint_0^T (delu')/(delt)dt)^0+1/Tint_0^T (delbaru^2)/(delx)dt+1/Tint_0^T (del(u')^2)/(delx)dt+cancel(1/Tint_0^T (del2(baru u'))/(delx)dt)^0+1/Tint_0^T (delbaru barv)/(delx)dt+cancel(1/Tint_0^T (delbaru v')/(delx)dt)^0+cancel(1/Tint_0^T (delu' barv)/(delx)dt)^0+1/Tint_0^T (delu'v')/(delx)dt=-1/Tint_0^T 1/rho (delbarP)/(delx)dt-cancel(1/Tint_0^T 1/rho (delP')/(delx)dt)^0+1/Tint_0^T v(del^2baru)/(dely^2)dt+cancel(1/Tint_0^T v(del^2u')/(dely^2)d)^0`
`=> 1/Tint_0^T (delbaru)/(delt)dt+1/Tint_0^T (delbaru^2)/(delx)dt+1/Tint_0^T (del(u')^2)/(delx)dt+1/Tint_0^T (delbaru barv)/(delx)dt+1/Tint_0^T (delu'v')/(delx)dt=-1/Tint_0^T 1/rho (delbarP)/(delx)dt+1/Tint_0^T v(del^2baru)/(dely^2)dt`
Upon integrating the equation we get
`=> (delbaru)/(delt)+(delbaru^2)/(delx)+1/Tint_0^T(del(u')^2)/(delx)dt+(del(baru barv))/(dely)+1/Tint_0^T(del(u'v'))/(dely)dt=-1/rho (delbarP)/(delx)+v(del^2baru)/(dely^2)`
`=> (delbaru)/(delt)+(delbaru^2)/(delx)+(del(baru barv))/(delx)=-1/rho (delbarP)/(delx)+v(del^2baru)/(dely^2)-1/Tint_0^T(del(u')^2)/(delx)dt-1/Tint_0^T(del(u'v'))/(dely)dt`
Through the above equations, it is clear that the momentum equations obtained are mean values of components is similar to momentum equations of actual transport properties except for the integration part.
Let us consider that u' value is fluctuating between -1 and 1 then by averaging the value of u' will be zero but if we consider the same case for the `(u')^2` then the square of -1 and 1 will be between 0 and 1 which when further integrated will be non zero.
`1/T int _0^T(del(u')^2)/(delx)dt !=0`
Similarly, in the term, `1/Tint_0^T (del(u'v'))/(dy) it's the fluctuations along the x-axis u' and the fluctuations along the y-axis v' are correlated. That is when a low momentum y-component is flowing in the positive y-direction then the high momentum (i.e high velocity) x-component gets slowed down which is caused by v'. So, the positive increase in the v' is going to cause a negative result in u', and the product of two will come out to be negative value only. Hence, this term will also not equal to zero.
`=> (delbaru)/(delt)+(delbaru^2)/(delx)+(del(baru barv))/(delx)=-1/rho (delbarP)/(delx)+v(del^2baru)/(dely^2)-1/Tint_0^T(del(u'v'))/(dely)dt-1/Tint_0^T(del(u')^2)/(delx)dt`
`(delbaru)/(delt)` = change with respect to time
``(delbaru^2)/(delx)+(del(baru barv))/(dely)` = convective term`
`-1/rho (delbarP)/(delx)` = pressure gradient
`v(del^2baru)/(dely^2)` = shear stress
``1/Tint_0^T(del(u')^2)/(dely)dt-1/Tint_0^T(del(u'v'))/(dely)dt` = Reynolds stres`
Applying average terms in integral values.
`=> (del baru)/(delt)+(del baru^2)/(delx)+(del(baru barv))/(dely)=-1/rho (del barP)/(delx)+v(del^2baru)/(dely^2)-(del(bar(u'v')))/(dely)-(del(bar(u')^2))/(delx)`
`=> (del baru)/(delt)+(del baru^2)/(delx)+(del(baru barv))/(dely)=-1/rho (del barP)/(delx)+v(del^2baru)/(dely^2)+1/rho[(del(-rho(bar(u'v'))))/(dely)-(del(-rho(bar(u'))^2))/(delx)]`
Stress Term = `mu nabla^2u_i - rho(del(bar(u'_i u'_j)))/(delx_j)`
`mu nabla^2u_i` = Shear Stress gradient
`rho(del(bar(u'_i u'_j)))/(delx_j)` = Reynolds Stress Gradient
`(del)/(delx_j)[mu ((delu_i)/(delx_j))-rho(bar(u'_iu'_j))]`
Total Shear Stress = `mu ((delu_i)/(delx_j))- rho(bar(u'_iu'_j))`
`=> mu ((delu_i)/(delx_j))` = Laminar Stress
`=> rho(bar(u'_iu'_j))` = Turbulent Stress
Reynolds Stress: It is a component of the stress tensor in the fluid which is obtained by averaging operation over the Navier Stokes equation used to account for turbulent fluctuation in the fluid.
Turbulent Viscosity: It is the function of fluctuations in velocity and is influenced by the randomness of flow. It is non-linear and difficult to estimate. It depends upon flow regimes and is not a property. It is dominant in turbulent flow regimes with high Reynolds numbers.
Molecular Viscosity: It is the resistance to the movement of one layer of fluid over the adjacent layer. It is linear and easy to measure. It is the fluid property and independent of geometry.
The ratio of turbulent to molecular viscosity gives an indication of how strong the Reynolds stress is when compared to the molecular stresses. Typically, an eddy viscosity ratio of more than 100 to 1000 indicates a turbulent flow. Although, the evolution of turbulence models has revolved around the definition of eddy viscosity. Hence, different models may have different values of it.
Transport Equation:
`rhoC_p[(delphi)/(delt)+u(delphi)/(delx)+v(delphi)/(dely)]=K[(del^2phi)/(delx^2)+(del^2phi)/(dely^2)]+mu((delu)/(dely))^2`
Applying decomposed equations to the quantities
`=> rhoC_p[(del(barphi+phi'))/(delt)+(baru+u')(del(barphi+phi'))/(delx)+(barv+v')(del(barphi+phi'))/(dely)]=K[(del^2(barphi+phi'))/(delx^2)+(del^2(barphi+phi'))/(dely^2)]+mu((del(baru+u'))/(dely))^2`
`=> rhoC_p[(delbarphi)/(delt)+(delphi')/(delt)+baru(delbarphi)/(delx)+baru(delphi')/(delx)+u'(delbarphi)/(delx)+u'(delphi')/(delx)+barv(delbarphi)/(dely)+barv(delphi')/(dely)+v'(delbarphi)/(dely)+v'(delphi')/(dely)] = K [(del^2barphi)/(delx^2)+(del^2phi')/(delx^2)+(del^2barphi)/(dely^2)+(del^2phi')/(dely^2)]+mu[((delbaru)/(dely))^2+((delu')/(dely))^2+2(delbaru)/(dely)(delu')/(dely)]`
applying time step T
`=> rhoC_p[1/Tint_0^T(delbarphi)/(delt)dt+1/Tint_0^T(delphi')/(delt)dt+1/Tint_0^Tbaru(delbarphi)/(delx)dt+1/Tint_0^Tbaru(delphi')/(delx)dt+1/Tint_0^Tu'(delbarphi)/(delx)dt+1/Tint_0^Tu'(delphi')/(delx)dt+1/Tint_0^Tbarv(delbarphi)/(dely)dt+1/Tint_0^Tbarv(delphi')/(dely)dt+1/Tint_0^Tv'(delbarphi)/(dely)dt+1/Tint_0^Tv'(delphi')/(dely)dt] = K [1/Tint_0^T(del^2barphi)/(delx^2)dt+1/Tint_0^T(del^2phi')/(delx^2)dt+1/Tint_0^T(del^2barphi)/(dely^2)dt+1/Tint_0^T(del^2phi')/(dely^2)dt]+mu[1/Tint_0^T((delbaru)/(dely))^2dt+1/Tint_0^T((delu')/(dely))^2dt+1/Tint_0^T2(delbaru)/(dely)(delu')/(dely)dt]`
`=> rhoC_p[1/Tint_0^T(delbarphi)/(delt)dt+cancel(1/Tint_0^T(delphi')/(delt)dt)^0+1/Tint_0^Tbaru(delbarphi)/(delx)dt+cancel(1/Tint_0^Tbaru(delphi')/(delx)dt)^0+cancel(1/Tint_0^Tu'(delbarphi)/(delx)dt)^0+1/Tint_0^Tu'(delphi')/(delx)dt+1/Tint_0^Tbarv(delbarphi)/(dely)dt+cancel(1/Tint_0^Tbarv(delphi')/(dely)dt)^0+cancel(1/Tint_0^Tv'(delbarphi)/(dely)dt)^0+1/Tint_0^Tv'(delphi')/(dely)dt] = K [1/Tint_0^T(del^2barphi)/(delx^2)dt+cancel(1/Tint_0^T(del^2phi')/(delx^2)dt)^0+1/Tint_0^T(del^2barphi)/(dely^2)dt+cancel(1/Tint_0^T(del^2phi')/(dely^2)dt)^0]+mu[1/Tint_0^T((delbaru)/(dely))^2dt+1/Tint_0^T((delu')/(dely))^2dt+cancel(1/Tint_0^T2(delbaru)/(dely)(delu')/(dely)dt)^0`
`=> rhoC_p[1/Tint_0^T(delbarphi)/(delt)dt+1/Tint_0^Tbaru(delbarphi)/(delx)dt+1/Tint_0^Tu'(delphi')/(delx)dt+1/Tint_0^Tbarv(delbarphi)/(dely)dt+1/Tint_0^Tv'(delphi')/(dely)dt] = K [1/Tint_0^T(del^2barphi)/(delx^2)dt+1/Tint_0^T(del^2barphi)/(dely^2)dt]+mu[1/Tint_0^T((delbaru)/(dely))^2dt+1/Tint_0^T((delu')/(dely))^2dt]`
`rhoC_p[(delbarphi)/(delt)+baru(delbarphi)/(delx)+1/Tint_0^Tu'(delphi')/(delx)dt+barv(delbarphi)/(dely)+1/Tint_0^Tv'(delphi')/(dely)dt] = K [(del^2barphi)/(delx^2)+(del^2barphi)/(dely^2)]+mu[((delbaru)/(dely))^2+1/Tint_0^T((delu')/(dely))^2dt]`
From the above equation, it is clear that the transport equation obtained with mean values of components si similar to the transport equation of actual transport properties except with the red terms.
By replacing the integral terms by averaging values we get
`=> rhoC_p[(delbarphi)/(delt)+baru(delbarphi)/(delx)+bar(u'(delphi')/(delx))+barv(delbarphi)/(dely)+v'bar((delphi')/(dely))] = K [(del^2barphi)/(delx^2)+(del^2barphi)/(dely^2)]+mu[((delbaru)/(dely))^2+((delbaru')/(dely))^2]`
Adding `bar(phi(del(u'))/((delx)))+bar(phi(del(v'))/((delx)))=0` on the LHS
`=> rhoC_p[(delbarphi)/(delt)+baru(delbarphi)/(delx)+bar(u'(delphi')/(delx))+barv(delbarphi)/(dely)+v'bar((delphi')/(dely))+[bar(phi(del(u'))/((delx)))+bar(phi(del(v'))/((dely)))]] = K [(del^2barphi)/(delx^2)+(del^2barphi)/(dely^2)]+mu[((delbaru)/(dely))^2+((delbaru')/(dely))^2]`
`=> rhoC_p[(delbarphi)/(delt)+baru(delbarphi)/(delx)+barv(delbarphi)/(dely)+bar(phi(del(u'phi'))/((delx)))+bar(phi(del(v'phi'))/((dely)))] = K [(del^2barphi)/(delx^2)+(del^2barphi)/(dely^2)]+mu[((delbaru)/(dely))^2+((delbaru')/(dely))^2]`
`bar(phi(del(u'phi'))/((delx)))+bar(phi(del(v'phi'))/((dely)))` = Turbulent Flu`
`(del^2barphi)/(delx^2)+(del^2barphi)/(dely^2)` = Molecular Flu`
`mu((delbaru)/(dely))^2` = Viscous Stress Dissipatio`
`mu((delbaru')/(dely))^2` = Turbulent Stress Dissipatio`
CONCLUSION:
- I would like to conclude that, the above RANS modeled Navier Stokes equation has some unknown quantities like time-mean correlations of fluctuating velocities and temperature.
- The solution cannot be attempted without additional relations or empirical modeling ideas.
- These forms are the basis for the modeling of any available time averaging turbulence models by using: K-epsilon of K-Omega.