CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Solving k-e model

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

Reply
 
LinkBack Thread Tools Display Modes
Old   December 9, 2010, 19:54
Default Solving k-e model
  #1
New Member
 
Alexander Ruth
Join Date: Dec 2010
Posts: 1
Rep Power: 0
Picknikka is on a distinguished road
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
Picknikka is offline   Reply With Quote

Old   December 10, 2010, 03:19
Smile
  #2
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by Picknikka View Post
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.
ignat is offline   Reply With Quote

Old   December 15, 2010, 09:59
Default
  #3
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
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
jyothishkumar is offline   Reply With Quote

Old   December 15, 2010, 15:35
Default
  #4
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by jyothishkumar View Post
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.
ignat is offline   Reply With Quote

Old   December 16, 2010, 06:40
Default
  #5
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
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
jyothishkumar is offline   Reply With Quote

Old   December 18, 2010, 13:59
Default
  #6
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by jyothishkumar View Post
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 u_w=0 does not contradict with condition \tau_w =(any value). You can use either of them.

Since wall shear stress \tau_w 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 u_w=0.
But if closest node is far from the wall (y+>10) then better choice for BC is \tau_w
ignat is offline   Reply With Quote

Old   December 22, 2010, 07:30
Default
  #7
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
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
jyothishkumar is offline   Reply With Quote

Old   December 22, 2010, 10:08
Default
  #8
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by jyothishkumar View Post
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 \tau_w you should solve the equation:

\frac{u_1}{u_*}=5.75\lg{\frac{y_1 u_*}{\nu}}+5.5

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

Then \tau_w=\rho u_*^2. (1)
On the other hand

\tau_w=(\mu+\mu_{turb})\frac{u_1 - u_{-1}}{2y_1} (2)

where u_{-1} is the false velocity at the false node outside the domain which is placed at the same distance
y_1. From (1) and (2) you find u_{-1} and use it as boundary condition for momentum equation.
Look picture
Attached Images
File Type: jpg bc.JPG (7.1 KB, 15 views)

Last edited by ignat; December 22, 2010 at 15:20.
ignat is offline   Reply With Quote

Old   December 22, 2010, 15:17
Default
  #9
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by jyothishkumar View Post
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.
ignat is offline   Reply With Quote

Old   December 23, 2010, 05:45
Default
  #10
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
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
jyothishkumar is offline   Reply With Quote

Old   December 23, 2010, 10:33
Default
  #11
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Hi jyothishkumar

It was really misunderstanding, sorry. That's ok.

Best regards,
ignat is offline   Reply With Quote

Old   December 29, 2010, 02:14
Default
  #12
New Member
 
Amir
Join Date: Mar 2009
Posts: 8
Rep Power: 8
Amiragha72 is on a distinguished road
Hi
Use Teach code as a pattern.
Amir
Amiragha72 is offline   Reply With Quote

Old   December 29, 2010, 04:38
Default
  #13
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
hi amir,

what is this Teach code. Is it something to tell about turbulence model implementation ?

thanks

jyothish
jyothishkumar is offline   Reply With Quote

Old   March 31, 2011, 03:14
Default
  #14
Member
 
jk
Join Date: Jun 2009
Posts: 64
Rep Power: 8
jyothishkumar is on a distinguished road
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
jyothishkumar is offline   Reply With Quote

Old   March 31, 2011, 07:20
Default
  #15
New Member
 
Alexey
Join Date: May 2009
Location: St.Petersburg, Russia
Posts: 22
Rep Power: 8
ignat is on a distinguished road
Quote:
Originally Posted by jyothishkumar View Post
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.
ignat is offline   Reply With Quote

Old   August 10, 2012, 06:09
Default Prandtl mixing length implementation
  #16
Member
 
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 5
michaelmoor.aero is on a distinguished road
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
michaelmoor.aero is offline   Reply With Quote

Old   August 10, 2012, 13:18
Default
  #17
Member
 
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 5
michaelmoor.aero is on a distinguished road
Quote:
Originally Posted by michaelmoor.aero View Post
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?
michaelmoor.aero is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Darcy-Forchheimer law for specifying Porous Zones Ger_US OpenFOAM Running, Solving & CFD 22 November 17, 2014 04:28
Low Mach number Compressible jet flow using LES ankgupta8um OpenFOAM Running, Solving & CFD 7 January 15, 2011 14:38
reactingFoam: adiabatic flame temperature rFkhemek OpenFOAM 2 July 31, 2009 06:42
Modeling in micron scale using icoFoam m9819348 OpenFOAM Running, Solving & CFD 7 October 27, 2007 00:36
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58


All times are GMT -4. The time now is 21:33.