 In the previous lecture, we appreciated that in a turbulent layer the inner layer can be nearly universal at least when pressure gradients are moderate, but the outer layer universality is very difficult to establish. It is here that we principally need the turbulence models and that is the topic of discussion today. So, I will explain what the main task is of modeling. I will take up the so called ad viscosity models and in that there are three variants. One is called the general mixing length model. The second one is higher and also number one equation model and higher and also number two equation model. You will understand the meanings of this as we go along. So, in multidimensional turbulent flows, even if inner layer universality is exploited, Rans equations must be solved in the outer layers through modeling and that means that the turbulent stresses u i prime u j prime and heat fluxes rho c p u prime i prime t prime must be modeled to recover lost information through averaging. Remember these terms arise because of Reynolds averaging. This recovery however must be carried out in a general way so that the model need not be changed from one flow situation to another. Now, it is this imparting absolute generality to a turbulence model has however been found to be quite difficult task, but nonetheless we will see as we go along what progress has been made. Thirdly, the model must be economical that is the computational expense or the ease must be must not be very much in excess of that which would be required for computation of say a laminar flow under the same situation. So, with these premises we begin to look at what are the possible turbulence models. There are two main approaches to modeling a stress. One is to say that the stress is related to the mean strain rate through a property called turbulent or ad viscosity. This is very analogous to Stokes's stress strain relationships. Alternatively, the stresses can be recovered from solution of transport equations for the quantity u i prime u j prime in which the convection and diffusion of this quantity is principally balanced by the rates of its reduction and dissipation. In this lecture, I am going to concentrate on the so called ad viscosity turbulence models the first type. The second one we will take up in the next lecture. So, the main idea is behind turbulent viscosity mu t is to model a turbulent stress minus rho u prime u i prime u j prime. Exactly analogous to the laminar stress is except that mu is replaced now by a turbulent viscosity mu t. This is the strain rate and then minus 2 3 rho kinetic energy e into delta i j. Delta i j is the chronicle of delta. Its value is 0 when i is not equal to j, but when i is equal to j it is equal to 1. Remember mu t is the property of the flow and not that of the fluid. Mu t is also isotropic that is mu t in all directions x i is identical at a point, but at different points in the flow its magnitude may vary. The term involving chronicle delta this one is necessary in turbulent modeling because the sum of the normal stresses when i is equal to j. Remember this will be d u i d x i and this will be also d u i d x i. So, this will become 2 mu d u i d x i, but that is equal to 0 in a incompressible flow. Therefore, we are left with minus rho u i prime square equal to that quantity which means it should equal to rho e. That is why that quantity is included in the equations just to make sure that the kinetic energy is balanced in this representation of a stress and strain relationship. Using this model the number of unknowns u i and p in the Rans equation equals the number of unknowns, but then the value of mu t now must be provided somehow and that is the modeling part characterizing mu t. So, from kinetic theory laminar viscosity mu is written in dimensionally correct form as mu times rho into the mean free path length l into the average molecular velocity u mole over bar. Analogously we say mu t would be modeled as rho times some length scale multiplied by a velocity scale u i u dash and these are representative length and velocity scales of turbulent fluctuations. So, then the main task now is in modeling mu t is to model l and u i prime. Let us see how we can do that in a as universal manner as possible. So, for example, the simplest model of u i prime which is called the general mixing line model is to say it is proportional to some quant length scale l m multiplied by the mean velocity gradient and in a three dimensional flow this would be given as so that mu t would become rho m times l m multiplied by l m into d u i by d x j into d u i by d x j into d u i by d x j d u j by d x i raised to point five with summation. So, that this is really phi v if you like or the average mean velocity gradient l is equal to l m is called the Prentzels mixing length and l m is given as you will recall from our lecture on near wall flows kappa times near wall distance y or n in here into 1 minus e raised to minus xi where xi is 1 by 26 n plus that is n plus by 26 and l m the is equal to point 2 into k in the kappa into r where r may be the radius of the pipe or it may be the boundary layer thickness if n plus is greater than 26. Now, the wall shear stress that is required here is calculated simply by the product of mu multiplied by the total velocity gradient at the wall that is the normal velocity gradient at the wall. So, the model is now implementable because the velocity distributions would be available from solution of momentum equations and the stresses in them would be modeled through mu t times this and therefore, tau wall can always be recovered. For two dimensional boundary layers of course, this mu t which we had represented here for the general three dimensional flows would become simply rho l m square d u by d y. Now for inner and outer layer boundary layers from lecture 24 we say l m is equal to kappa y 1 minus exponential of minus y plus by a plus and a plus was given by this quantity where a and b were sensitized to value of the suction or blowing parameter or the pressure gradient parameter p plus and of course, when y plus exceeds a plus it would simply be equal to what we have shown here point 2 times kappa times the length dimension of the flow. When we encounter free shear layers like the a turbulent jet or a turbulent wake then there is no wall present in that case how do we specify a mixing length well it is specified as l m time equal to beta times the half width of the jet. Let us say the jet emerging then this would be the velocity profile and this would be u max and we say wherever u divided by u max is equal to half that is the distance y half this is called the half jet width in a free shear layer likewise for a wake also. As I said y half is the half width of the shear layer and beta takes different values this has been observed from experiments as well as has been tuned by numerical experiments. Beta equal to 0.225 if the jet was a plane jet if it was a round jet it would be 0.1875 and if it was a plane wake it would be 0.4. These are the typical values for free shear layers this is just as an aside we are not dealing here in this course with the free shear layers but then the less it is useful to document this prescription. Remember there is no presence of wall in a free shear layer and therefore there is no notion of a y plus. So with this prescription l m is now available and therefore the equations can be solved as I said earlier. The mixing length model unfortunately predicts that rho u prime u i prime u j prime will be 0 where the strain rate S i j will be 0 as you can appreciate from this figure if the mean strain rate was 0 at any point mu t would automatically be 0. In many situations this is not found to be the case for example if I were to take the case of a annulus let us say this is the inner wall and that is the outer wall then we would have a velocity profile which will look something like and let us say the outer wall is rough and the inner wall is smooth this is the smooth wall and we would expect a velocity profile to be something of this variety something of this variety. So the plane of zero shear which is here let us say this is where d u by d r will be equal to 0 tends to be closer to and if I were to measure now in this case rho u prime v prime then I would find that the plane of zero shear tends to be closer this is where rho u prime v prime will be 0 then to the smooth wall then the plane of zero velocity gradient. So in fact the shear stress does not really go with the velocity or the velocity gradient or the strain rate and there is a separation here whereas the mixing length model would predict zero shear stress here also whereas the true zero shear turbulent stress is 0 here. Now in order to meet with such asymmetric situations that we need refinements in the models. If we for example represent velocity scale u i prime equal to square root of kinetic energy then mu t will be called rho times L square root of E where L is now has a different meaning and it would be called the integral length scale which we defined earlier and the turbulent kinetic energy equation E must be obtained from turbulent kinetic energy equation which is as you will recall given by the convection terms here that this is the turbulent diffusion due to pressure fluctuations and velocity fluctuation and then there will be the laminar diffusion. This should be the rate at which energy is extracted from the mean motion to produce turbulence and this with the rate at which the energy will be dissipated that we called rho times epsilon. The model of the one equation type goes like this it is mu t is equal to rho L square root of E where E would be obtained from a differential equation whereas L would be now specified algebraically again like in the mixing length and because we have to now solve this equation additional equation along with the Rans equation it is called the one equation model. Now the turbulent term cannot be directly measured turbulent diffusion term which contains pressure fluctuations and the triple velocity correlation these cannot be measured directly because no probe can simultaneously measure pressure fluctuations and velocity fluctuations and therefore, we need to model it. So, but noting the redistributive character of this term you will recall that when we did the energy balance we noted that this term is redistributed and therefore, we will assume a gradient diffusion for this. So, minus u j prime p prime rho u i prime u i prime by 2 would be simply mu put as mu t divided by sigma E d E by d x j where sigma E is the turbulent transfer number of turbulent kinetic energy. Now when r e t is high dissipation can be represented in terms of large scale fluctuation this is what we had observed when we looked at the formal aspects of turbulence and therefore, this mu times d u i prime by d x j prime which we had modeled as rho epsilon can be written as C d times rho times e raise to 3 by 2 divided by L. This is now the representative v dash cubed if you recall we had said that the rho times epsilon should be proportional to v dash cube divided by L and that is what I am writing here. So, v dash being square out of u we have said that that is equal to rho e raise to 3 by 2 by L. If we now replace minus rho u prime u j prime equal to mu t times s i j then the model turbulent kinetic energy equation would be simply convection terms here plus mu plus mu t divided by sigma which is into this gradient of e which would be the turbulent diffusion and then there would be the generation term energy production term by replacing this and then there would be the dissipation term where C d and sigma epsilon are expected to be universal constant when r e t is high that is away from the wall and beyond the transitional layer. How do we get C d and sigma epsilon? Recall that in the fully turbulent inner layer production is very nearly equal to dissipation as we saw in lecture 24 and the equilibrium conditions prevail. In this layer tau t is equal to mu t d v by d n and hence tau t d v by d n which would be the production term. This term would be the production term in the two dimensional boundary layer or very close to the wall minus equal to rho square root of e L d v by d n whole square equal to C d rho e raise to 3 by 2 by L which is the dissipation and where v is the velocity parallel to the wall and n is the normal distance. Now, in this layer we are also said that tau t is very much close to tau w that is the wall shear stress and hence tau w divided by rho divided by e would be u tau square by e equal to square root of C d. Now, you will recall that we had shown from the experimental data that tau w is approximately equal to 0.3 rho times e and therefore, C d takes the value about 0.09. Sigma epsilon sigma e is taken equal to 1 from numerical experiment in several flow situations like a pipe flow or a boundary layer and so on and so forth. Hence, the model turbulent kinetic energy equation can be solved. We must now specify the length scale L. How do we do that? Well, L is determined as follows. Consider equilibrium condition again production is equal to dissipation then mu t will be equal to rho e raise to and the definition of mu t is now rho e raise to half L which is the definition. If I combine these two equations and eliminate e from them, I would get mu t equal to C d raise to minus half rho L square d v by d L. Now, if I compare this equation with the equation of the mixing length model, then I would see that L would be you recall the mixing length model mu t equal to rho L m square into d u d y. So, L would now equal C d raise to 0.25 into L m which is 0.5477 L m and L m equal to kappa y times the function that we normally use. With these specifications L mu t C d n sigma epsilon d a turbulent kinetic energy solution along with the Rand's equations and in general flows however, further refinements become necessary and in the next slide we will show you why this is so. Although kinetic energy can equation can take care of the non-coincidence of 0 velocity gradient and the shear stress, it is not a very good approximation in some other situation as we shall see shortly. For the more general flows involving strong and variable pressure gradients or even recirculations or swirl or buoyancy effects, it is necessary to devise means for determining distribution of L in a multidimensional flow. So, recall the energy spectrum E k in wave number space with wave number k and if we divide this E k by k and integrate from 0 to infinity, an equation for L in the spectral space can indeed be derived since L is equal to 1 over E 0 to infinity E k by k d k. Unfortunately, this equation is not known as not tractable in physical space and therefore, we need to find some alternative means of determining the length scale. So, direct equation for a length scale is not possible which is the need in multidimensional flows with all the effects that I mentioned here, but then a length scale equation need not necessarily we have L as its itself as its dependent variable. Any combination of the form z equal to E raise to m L raise to n will suffice since E can be known from the solution of the model turbulent kinetic energy equation. There have been some proposals for equation for a z which is a composite variable E raise to m L raise to m. The model equation for z is postulated like this, there would be convection of z, there would be laminar and turbulent diffusion of z, there would be appropriately scaled production of z and there would be dissipation of z plus a additional source term as z. So, where p as you recall is the turbulent kinetic energy generation and C 1 C 2 are constants when R e t is high. Now, there have been some proposals for this for them for the first one took m equal to 3 by 2 and L equal to minus 1 and that would mean that z is equal to or proportional to E raise to 3 by 2 divided by L which is exactly would mean z would imply dissipation and therefore, it is called the dissipation equation. Equations for m equal to 1 and L equal to 1 has also been derived, but then that equation requires an additional s z term which I am not mentioned here. Likewise, you can also have a E raise to 1 divided by L square that is n equal to minus 2 and 1 which again requires another type of s z term and you can also have E raise to half multiplied by L which would be the turbulent viscosity equation itself, but that also requires an s z term. It is the dissipation rate equation is the only one which does not require s z and that is equal to 0. So, we have if we take z equal to E raise to 3 by 2 divided by L we would have a nice convection term, a diffusion term, a production term minus dissipation term with s z equal to 0 and this is the most preferred alternative because this s z equal to 0 term is an attractive feature whereas, in all other models s z has to be tuned as it were and has to be modeled separately. Therefore, the most preferred choice is the dissipation equation with E raise to 3 by 2 and L raise to minus 1. So, one solves the z equation along with the E equation and simply recovers the L that we need at each point in the flow so as to form the turbulent viscosity. How do we justify the dissipation equation about which I just mentioned? Now, at high r e t local isotropic prevails and an equation for nu d u i prime by d x j square can be derived by differentiating instantaneous terms on Navier-Stokes equation for u i dash with x j and then multiplying by 2 nu d u i dash by d x j. If you then time average this nu equation then you get an exact equation for epsilon the quantity that we have chosen for the second variable z and that equation looks like this. It is a convection term and then there is a these two terms constitute diffusion of epsilon firstly due to pressure gradient or the pressure fluctuations and secondly due to the velocity fluctuation. This is the laminar diffusion, this is the production of epsilon, this contributes to the production of kinetic energy and these two are the dissipation of epsilon itself. Now, the complex correlations can be discerned from the DNS. So, how can we model these complex terms that you see here? That has been done by looking at computations of direct numerical simulation. So, the diffusion term apparently can be modeled as d by d x k into c 3 time e square by epsilon d epsilon by d x k in which u j dash u j k dash has been taken proportional to e and this is an assumption in this modeling. The production term about which I mentioned this is the production term is modeled as simply c 1 times the u i dash u j dash multiplied by strain rate multiplied by the time scale epsilon by e and then the last two terms are simply modeled as c 2 times epsilon square by e. Remember this term would not be relevant in the highest wave numbers space, but relevant in the inertial sub range whereas, this term would be dominant in the highest wave numbers where dissipation really occurs. These two terms can be grouped together and modeled as c 2 times epsilon square by epsilon e. Remember all this modeling work has been done earlier was made possible from measured turbulence data, but now it is also possible from direct numerical simulation data. Therefore, we get the dissipation equation which reach like this. This is the convective term of the dissipation equation, this is the diffusion term, this is the turbulent diffusion, this is laminar diffusion and then we get the production term and the dissipation term itself. Mu t is rho square root e into L and therefore, this will simply become c d rho e square by epsilon because epsilon is proportional to e raise to 3 by 2 by L. Hence, the dissipation equation would be written as d by dx k of diffusion gradient of epsilon into c 1 times production minus c 2 times rho epsilon square by e and where u i dash u k the strain rate multiplied by the viscosity. This is the model that we choose to model the stress which appears here. Sigma epsilon is the turbulence frontal number for dissipation rate. It absorbs constant c d and c 3, this c 3 and having modeled u i dash u j dash by this as c d times u i dash u j dash by this c d times so the mu t is this whole term is replaced by mu t by sigma epsilon and c 3 and c d are absorbed near. So, this is then the final form of the epsilon equation which we shall choose. The so called two equation model is essentially this. There is the kinetic energy equation which you have already appreciated earlier and this would be the dissipation equation now written with u i dash u dash replaced by the strain rate and the boundary layer form of this would be this are the convection terms of e diffusion terms, generation term minus dissipation and for dissipation this would be the convection term, this will be the diffusion term and this will simply be mu t d u by d y whole square minus c 2 rho epsilon. The question now is how do we determine c 1 and c 2 and sigma epsilon? These are determined from observation of experimental data in decay of homogeneous turbulence. So, what is done is that in a wind tunnel you have you put a screen like so and the air flows through the screen and you simply because of the screen high production of energy would take place and the energy kinetic energy will go on with x it will go on declining. If you look at d e by d t in this case in this flow both production and diffusion are absent and v the normal velocity component is 0 and therefore, d e by d t would be simply you will see in this case d e by d t would be simply equal to minus rho epsilon. So, rate of change of energy is simply balanced by the dissipation and that is what I have written here and likewise d epsilon by d t would be equal to minus c 2 epsilon square by e. If we solve these two equations d e by d t with the observed fact that e is proportional to t to the power of minus n it follows from this d e by d t is equal to minus epsilon and d epsilon by d t is equal to minus c 2 times epsilon square by e and therefore, inter substitution shows that e should vary as t raise to minus n and if you look at the experimental data you find that n varies between 1 and 2 and 1 n into 1.2 at least in the early part for small times and for later times n is of the order of 1.8 or so which we shall look at a little later. So, in the initial times at least it is between n equal to 1 and 1.2 and therefore, the solution gives c 2 equal to n plus 1 divided by n simultaneous solution and therefore, taking n equal to 1.1 say we will get c 1 equal to about 1.91. To determine c 1 we consider the inner turbulent layer where convection is 0, nu t is very much greater than nu and production nu t d u by d y square is equal to dissipation and hence if epsilon equation this is 0 almost and therefore, production is equal to dissipation because that is also 0. Since convection is 0 we will get the total diffusion equal to that term and which we have modeled as e epsilon square equal to and remember production would now be nu t d u by d y whole square equal to dissipation and therefore, epsilon square by epsilon e into c 2 minus c 1. In the inner layer the logarithmic law suggests that u plus equal to l n e y plus by kappa gives d u by d y equal to u tau by kappa y and that is equal to tau all over rho nu t and that is also equal to u tau square over nu t which gives nu t equal to u tau kappa y and hence since production is equal to dissipation epsilon would be nu t d u by d y whole square equal to u tau by kappa y which is gives relationship between epsilon and y and therefore, d epsilon by d y would be minus u tau cube over k y square. So, if we now substitute d epsilon by d y here in this equation and also take mu t appropriately as shown here and d epsilon by d y as shown here. So, then you will get d by d y of u tau kappa y by sigma epsilon into this quantity equal to u tau square by sigma epsilon y square equal to all this or kappa square over sigma epsilon y square epsilon would be c 2 minus c 1 u tau square by e and therefore, c 1 would be equal to c 2 minus kappa square over sigma epsilon over c d and that would turn out to be equal to 1.44 because c d recall was equal to 0.09 kappa was 0.41 and sigma epsilon has been found from numerical experiments to be approximately 1.3. Therefore, c 1 will be equal to c 2 minus kappa square over this equal to 1.44. So, with this we conclude how c 1 and c 2 are determined. Remember we had said c 2 was equal to 1.91 and therefore, c 1 can now be determined from this and it turns out to be about 1.44. Now of course, these equations apply as I said to high turbulent Reynolds number region which means beyond the transitional layer and therefore, in a practical flow the layer up to the transitional layer must be somehow captured in a computation. That is done as I have shown here. So, in a the equations I have solved let us say by a finite difference procedure then if this is the wall then you will take the this will be the grid node at the wall, but the next grid node would be taken at a distance substantially away from the wall. So, that y p plus would be say anywhere between third and fourth t and 100 and you are sure that you are in the inner turbulent layer. So, the point p must be in the inner turbulent layer. So, that is the most important part and then you only solve for the outer layer which is beyond point p. So, all the boundary conditions must be applied at p. Now of course, we do not know the exact boundary conditions for e and epsilon at p, but they can be derived. So, for example, if we take tau wall equal to mu effective d u by d y p equal to mu effective u p by y p and hence mu effective by y p becomes simply rho u tau kappa ln e y plus p where u tau is equal to c d raised to point 25 e p raised to half, you will recall that we had shown this for the inner layer or the turbulent part of the inner layer and therefore, u tau can be estimated like this in rather than from mu d u d y equal to 0 because we are not computing up to this point. We are computing from this point onwards. u would be of course, obtained from the ranch equations at this point. So, in the e equation if we say production equal to that and equate it to an average dissipation over this length epsilon bar p would be equal to y p raised to minus 1 0 to y p epsilon d y. Substituting for epsilon equal to production we get that and that would integrate to simply this and therefore, the source term of e is simply mu effective u p by y p whole squared minus rho epsilon bar which is what I have shown here what epsilon bar is and epsilon p itself would be simply c d raised to 0.75 e p raised to half 1.5 into kappa y p. So, what is done is the source term at point p is modified in the k and epsilon equations. So, these are the source terms in the k and epsilon equation then at the point p these two source terms are simply modified to read as I have shown here. Now, with this the region where very sharp gradients of velocity kinetic energy and other quantities take place is completely circumvented and everything begins from here and we have essentially exploited the universality of the inner layer in dimensionless quantities. This gives enormous savings in computational time and therefore, it is routinely used in two equation turbulence model and this is called the wall function approach to solving the two equation turbulence model equations. In the next lecture I will consider further refinements for situations in which the eddy viscosity models do not quite apply in a precise manner and therefore, one has to turn to the modeling of the turbulent stresses themselves are obtaining turbulent stresses through transport equations for U i prime U j prime.