Solving k-e model

 Register Blogs Members List Search Today's Posts Mark Forums Read

 December 9, 2010, 19:54 Solving k-e model #1 New Member   Alexander Ruth Join Date: Dec 2010 Posts: 1 Rep Power: 0 Hi, I am writing a code with FORTRAN to solve the k-e model. I have a Problem with the eddy viscosity formulation. The flow I like to predict is one over a rotating disc. This flow is considered to disappear in the infinty and so at the end of the grid (far away from the disc) too. So k and e become zero far away from the disc....my problem is that in formulation of the eddy viscosity for the k-e model i have to calculate kČ/e....and this causes problems because k is in the infinity 1+10^-30 and e is 1*10^-33....so kČ/e is becoming something like 1000 but it should be zero. Has anybody an idea how I can solve this? Thanks! Alex

December 10, 2010, 03:19
#2
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by Picknikka Hi, I am writing a code with FORTRAN to solve the k-e model. I have a Problem with the eddy viscosity formulation. The flow I like to predict is one over a rotating disc. This flow is considered to disappear in the infinty and so at the end of the grid (far away from the disc) too. So k and e become zero far away from the disc....my problem is that in formulation of the eddy viscosity for the k-e model i have to calculate kČ/e....and this causes problems because k is in the infinity 1+10^-30 and e is 1*10^-33....so kČ/e is becoming something like 1000 but it should be zero. Has anybody an idea how I can solve this? Thanks! Alex
Specify eps some larger so that turbulent viscosity become much less then molecular one.

 December 15, 2010, 09:59 #3 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 Dear Mr Ignat, I have developed a 2d collocated pipe flow cfd code. I am trying to implement turbulence model into it (simple prandtl's mixing length model). But after implementation turbulence profile looks more peaky than the laminar profile at the exit. Can you suggest me some solution for this error. Thanks in advance

December 15, 2010, 15:35
#4
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by jyothishkumar Dear Mr Ignat, I have developed a 2d collocated pipe flow cfd code. I am trying to implement turbulence model into it (simple prandtl's mixing length model). But after implementation turbulence profile looks more peaky than the laminar profile at the exit. Can you suggest me some solution for this error. Thanks in advance
hi jyothishkumar.
As I know prandtl's mixing length model is too rough.
I recommend you to use Spalart-Allmaras model without ft2 term.
Look here http://turbmodels.larc.nasa.gov/spalart.html

It is simple one-equation model for one variable with no-problem boundary conditions. For pipe flows it gives good resuts.
Be careful with boundary conditions for velocity at the wall. If closest cell size > then 10 (in wall units y+) you should use wall-function - logarithmic profile.

 December 16, 2010, 06:40 #5 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 Dear Mr ignat, Thankyou for your message. One more doubt is that you said to be cautious at the wall boundary condition for velocity, confused whether to give the wall velocity zero or not. Actually from the log law layer expression are we finding the tou_w (wall shear stress) ? Can you please clear this doubt. Thanks

December 18, 2010, 13:59
#6
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by jyothishkumar Dear Mr ignat, Thankyou for your message. One more doubt is that you said to be cautious at the wall boundary condition for velocity, confused whether to give the wall velocity zero or not. Actually from the log law layer expression are we finding the tou_w (wall shear stress) ? Can you please clear this doubt. Thanks
Hi jyothishkumar!
Sorry for delay so forum was down ! (many thanks to cfdonline team for fast restoring)

In fact, that no-slip condition on the wall does not contradict with condition =(any value). You can use either of them.

Since wall shear stress is the momentum flux it can be viewed as boundary condition (BC) for momentum equation written in conservation form.

If closest mesh node is placed inside laminar sub-layer (which very-very thin) then better choice is .
But if closest node is far from the wall (y+>10) then better choice for BC is

 December 22, 2010, 07:30 #7 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 Dear Mr Ignat Extremely thankful for your reply. Can you please check the following procedure that i am going to implement in my laminar code Procedure pertaining to turbulence -------------------------------------- 1 Assume tou_w (wall shear stress) to some initial value 2 Then find y+ according to the formula y+ = row * ut * y/meu 3 Substitute this y+ value in the tou_w expression obtained from log-law layer Then add this tou_w to the momentum equation as a source term along with the dp/dx like this (neglecting ‘a s’ term) apUp=awUw+anUn+aeUe+(dp/dx) * volume + (tou_w x area) but for v momentum equation apVp=awVw+anVn+aeVe+asVs +(dp/dy)*volume Overall procedure as follows First initialization – u, v, p and the turbulence quantities like neucap (from spallart almaras 1 equation model) and the tou_w Apply boundary condition for u,v and p (say for a 2d duct flow it will be like velocity inlet, wall boundary (no slip) and the pressure exit. Also for the 1equation from the SA we need to give values for ‘neucap’. At the inlet just extrapolate from the interior of the domain. (like neucap(x_inlet)=neucap(x_inlet+1). Similarly at the exit it will be like neucap(x_outlet)=neucap(x_outlet-1). At the wall neucap = 0. Solve the momentum equation with neglecting ‘a s’ south direction flux and adding a source term (tou_w x area) only to the u momentum eqn. as mentioned in the previous heading (point 3) Solve the pressure poison from continuity (in case of incompressible flow) do all the corrections to the u, v and the p. solve the SA 1 equation using this new value of u, v and find the value of neucap find the value of meueff from the value of meut (obtained from SA equation) and the meu. repeat the above process till convergence. thanks in advance

December 22, 2010, 10:08
#8
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by jyothishkumar Dear Mr Ignat Extremely thankful for your reply. Can you please check the following procedure that i am going to implement in my laminar code Procedure pertaining to turbulence -------------------------------------- 1 Assume tou_w (wall shear stress) to some initial value 2 Then find y+ according to the formula y+ = row * ut * y/meu 3 Substitute this y+ value in the tou_w expression obtained from log-law layer Then add this tou_w to the momentum equation as a source term along with the dp/dx like this (neglecting ‘a s’ term) apUp=awUw+anUn+aeUe+(dp/dx) * volume + (tou_w x area) but for v momentum equation apVp=awVw+anVn+aeVe+asVs +(dp/dy)*volume Overall procedure as follows First initialization – u, v, p and the turbulence quantities like neucap (from spallart almaras 1 equation model) and the tou_w Apply boundary condition for u,v and p (say for a 2d duct flow it will be like velocity inlet, wall boundary (no slip) and the pressure exit. Also for the 1equation from the SA we need to give values for ‘neucap’. At the inlet just extrapolate from the interior of the domain. (like neucap(x_inlet)=neucap(x_inlet+1). Similarly at the exit it will be like neucap(x_outlet)=neucap(x_outlet-1). At the wall neucap = 0. Solve the momentum equation with neglecting ‘a s’ south direction flux and adding a source term (tou_w x area) only to the u momentum eqn. as mentioned in the previous heading (point 3) Solve the pressure poison from continuity (in case of incompressible flow) do all the corrections to the u, v and the p. solve the SA 1 equation using this new value of u, v and find the value of neucap find the value of meueff from the value of meut (obtained from SA equation) and the meu. repeat the above process till convergence. thanks in advance
Hi!

Is this algorithm your own invention?
Sorry,but I think it does't work. Don't think up your own methods. Time is money. The better way is to study well known and robust methods for calculation of incompressible flows. I recommend you to use some staggered grid algorithm, for example this

Kim, J., Moin, P., Application of Fractional-Step Method to Incompressible Navier-Stokes Equations. // J.Comp.Phys. 1985, v.59, pp. 308-323

or some others which use staggered grid, but not collocated one because for collocated grids there is a problem of pressure-velocity link.

Some about boundary condition. I do the following way.

In order to calculate you should solve the equation:

and find where is known velocity at closest node, is the distance from the closest node to the wall and is the molecular kinematic viscosity.
This is transcendental equation and it is simply solved by Newton iterations
with initial guess for small value, for example 0.00001

Then . (1)
On the other hand

(2)

where is the false velocity at the false node outside the domain which is placed at the same distance
. From (1) and (2) you find and use it as boundary condition for momentum equation.
Look picture
Attached Images
 bc.JPG (7.1 KB, 17 views)

Last edited by ignat; December 22, 2010 at 15:20.

December 22, 2010, 15:17
#9
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by jyothishkumar At the inlet just extrapolate from the interior of the domain. (like neucap(x_inlet)=neucap(x_inlet+1).
This is wrong. At the INFLOW boundary neucap(x_inlet) must be specified by the user in accordance with upwind principle.
Your variant is correct for OUTFLOW boundary only.

 December 23, 2010, 05:45 #10 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 Hi, Thanks. I think you have mistaken me. It was just a outline only. Actually i have taken care of all those pressure-velocity coupling problems. Infact i had a systematic start in coding with 1d staggered solver then 2d staggered, 1d colocated and then 2d colocated. Now for this 2d colocated code i am trying to implement the turbulence model. (Basically mine is an incompressible flow solver). I will use your ideas to implement the boundary condition for the turbulence. thankyou jyothishkumar

 December 23, 2010, 10:33 #11 New Member   Alexey Join Date: May 2009 Location: St.Petersburg, Russia Posts: 23 Rep Power: 10 Hi jyothishkumar It was really misunderstanding, sorry. That's ok. Best regards,

 December 29, 2010, 02:14 #12 New Member   Amir Join Date: Mar 2009 Posts: 8 Rep Power: 10 Hi Use Teach code as a pattern. Amir

 December 29, 2010, 04:38 #13 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 hi amir, what is this Teach code. Is it something to tell about turbulence model implementation ? thanks jyothish

 March 31, 2011, 03:14 #14 Member   jk Join Date: Jun 2009 Posts: 64 Rep Power: 10 Dear Mr Ignat, I am trying with k epsilon turbulence implementation in my code (std k epsilon). I am just following the malalasekara book (latest edition) for this implementation. Or please tell me the ideas like what is the source term to be added in momentum eqn etc etc. how to handle the source term in k epsilon equation. Whether to solve for k equation for the first node or can i use directly like k = utou^2/sqrt(cmu). Thanks in advance jyothish

March 31, 2011, 07:20
#15
New Member

Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 23
Rep Power: 10
Quote:
 Originally Posted by jyothishkumar Dear Mr Ignat, I am trying with k epsilon turbulence implementation in my code (std k epsilon). I am just following the malalasekara book (latest edition) for this implementation. Or please tell me the ideas like what is the source term to be added in momentum eqn etc etc. how to handle the source term in k epsilon equation. Whether to solve for k equation for the first node or can i use directly like k = utou^2/sqrt(cmu). Thanks in advance jyothish
Dear jyothishkumar,

As I know standard k-e model is very stiff so both of equations have to be solved implicitly.
Both k-equation source term and e-equation one must be linearized and approximated at n+1 time level. Moreover, diffusion terms in both equations have to be solved implicitly too.
This way usually can avoid negative k & e values.

wish good luck.

 August 10, 2012, 06:09 Prandtl mixing length implementation #16 Member   Michael Moor Join Date: May 2012 Location: Ireland Posts: 30 Rep Power: 7 Hi everyone, I have succesfully created a solver in Matlab for 2d STEADY incompressible Poiseulle flow, using a backward staggered grid, FVM, SIMPLE algorithm and hybrid differencing scheme, with good results to Hagen-Poiseulle velocity profile. Boundary conditions are inlet-outlet, with the walls being brought into the domain the wall shear stress, (I have stuck rigidly to Versteeg). I am now at the stage of adding turbulence models. I understand from Patankar, that it is best to leave out the Reynolds stress terms in the momentum equations explicitly, but to the use the effective viscosity, which is the dynamic + turbulent. (mu_eff =mu + mu_t). Turbulence modelling: In the RANS and 2D flow, the only relevant reynold stress is tau_xy, which we then discretise -d(tau_xy)/dy for a cell volume, i.e. (tau_xy*A)n - (tau_xy*A)s, where n & s are the north and south FACES, and tau_xy, for the Prandtl mixing length is: rho*mixing_length^2*du/dy^2. Again, du/dy is ((u)n - (u)s)/deltayp. Similarly for the v-momentum, except now (tau_xy*A)e - (tau_xy*A)w.... du/dy is still the same. I feel ashamed having to ask about discretisation! lol! I just want to be sure and do not have time for second guessing myself!! mixing length for pipe flow is L[0.14-0.08(1-y/L)^2 - 0.06(1-y/L)^4], where y is the distance from the node to the wall , L is the radius/ half channel width. Therefore, for each u- and v- cell, i must calculate: 1. the mixing length, with y being the distance from the wall to the u-cell centre 2. the velocity gradient ((u)n - (u)s)/deltayp, where the value of (u)n or (u)s will be linearly approximated by, e.g. (u)n = (u(i,J) + u(i,J+1))/2 3. 1 & 2 allows me to calculate tau_xy 4. now, to get the value of tau_xy at cell faces, use central differencing? 5. along wall nodes, do i use u+ instead of up to calculate the shear stress and add it in the source serm b? It is SUPPOSED to be quick and easy to implement... so I appreciate any assistance. following this, i want to add k-epsilon Best regards, Michael 0565113@gmail.com

August 10, 2012, 13:18
#17
Member

Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 7
Quote:
 Originally Posted by michaelmoor.aero Hi everyone, I have succesfully created a solver in Matlab for 2d STEADY incompressible Poiseulle flow, using a backward staggered grid, FVM, SIMPLE algorithm and hybrid differencing scheme, with good results to Hagen-Poiseulle velocity profile. Boundary conditions are inlet-outlet, with the walls being brought into the domain the wall shear stress, (I have stuck rigidly to Versteeg). I am now at the stage of adding turbulence models. I understand from Patankar, that it is best to leave out the Reynolds stress terms in the momentum equations explicitly, but to the use the effective viscosity, which is the dynamic + turbulent. (mu_eff =mu + mu_t). Turbulence modelling: In the RANS and 2D flow, the only relevant reynold stress is tau_xy, which we then discretise -d(tau_xy)/dy for a cell volume, i.e. (tau_xy*A)n - (tau_xy*A)s, where n & s are the north and south FACES, and tau_xy, for the Prandtl mixing length is: rho*mixing_length^2*du/dy^2. Again, du/dy is ((u)n - (u)s)/deltayp. Similarly for the v-momentum, except now (tau_xy*A)e - (tau_xy*A)w.... du/dy is still the same. I feel ashamed having to ask about discretisation! lol! I just want to be sure and do not have time for second guessing myself!! mixing length for pipe flow is L[0.14-0.08(1-y/L)^2 - 0.06(1-y/L)^4], where y is the distance from the node to the wall , L is the radius/ half channel width. Therefore, for each u- and v- cell, i must calculate: 1. the mixing length, with y being the distance from the wall to the u-cell centre 2. the velocity gradient ((u)n - (u)s)/deltayp, where the value of (u)n or (u)s will be linearly approximated by, e.g. (u)n = (u(i,J) + u(i,J+1))/2 3. 1 & 2 allows me to calculate tau_xy 4. now, to get the value of tau_xy at cell faces, use central differencing? 5. along wall nodes, do i use u+ instead of up to calculate the shear stress and add it in the source serm b? It is SUPPOSED to be quick and easy to implement... so I appreciate any assistance. following this, i want to add k-epsilon Best regards, Michael 0565113@gmail.com
Or am i making this a lot more complicated? i.e. we just need to find mu_t = rho*mixing length^2*du/dy. du/dy is (u(norh face) - u(south face))/deltayp. mixing length is a function of height from wall. Then, if y+ <= 11.63, i make no changes. However, if y+=> 11.63, i use u+ value for wall shear stress, i.e. u+ = ln(Ey+)/kappa. and this wall shear stress only applies along the nodes adjacent to the wall?

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Ger_US OpenFOAM Running, Solving & CFD 31 August 5, 2017 03:00 ankgupta8um OpenFOAM Running, Solving & CFD 7 January 15, 2011 14:38 rFkhemek OpenFOAM 2 July 31, 2009 06:42 m9819348 OpenFOAM Running, Solving & CFD 7 October 27, 2007 00:36 msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58

All times are GMT -4. The time now is 06:39.