One dimensional k epsilon implementation from scratch

salazardetroya · July 19, 2015, 11:48

Hello,

I am trying to implement a 1D k- $\varepsilon$ turbulent model for a turbidity current, hence the conservation equation for $c$. I'm solving for the variables $u,c,k$ and $\varepsilon$. The remaining terms are parameters.
$\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + C_{\mu} \frac{k^2}{\varepsilon} \right) \frac{\partial u }{\partial z} \right] + \frac{\rho_s - \rho_w}{\rho_w} c g_x$

$\frac{\partial c}{\partial t} = v_s \frac{\partial c}{\partial z} + \frac{\partial}{\partial z} \left[ \left( \nu + C_{\mu c} \frac{k^2}{\varepsilon} \right) \frac{\partial c }{\partial z} \right]$

$\frac{\partial k}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{C_{\mu} }{\sigma_k} \frac{k^2}{\varepsilon} \right) \frac{\partial k }{\partial z} \right] + C_{\mu} \frac{k^2}{\varepsilon} \left(\frac{\partial u }{\partial z}\right)^2 - \varepsilon + \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} C_{\mu} \frac{k^2}{\varepsilon} \frac{\partial c }{\partial z}$

$\frac{\partial \varepsilon}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{C_{\mu} }{\sigma_{\varepsilon}} \frac{k^2}{\varepsilon} \right)\frac{\partial \varepsilon }{\partial z}\right] + c_{\varepsilon 1} \frac{\varepsilon}{k}$
$\left( C_{\mu} \frac{k^2}{\varepsilon} \left(\frac{\partial u }{\partial z}\right)^2 + c_{\varepsilon 3} \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} C_{\mu} \frac{k^2}{\varepsilon} \frac{\partial c }{\partial z} \right) - C_{\varepsilon 2} \frac{ \varepsilon^2}{k}$

So far, I am using finite difference. Upwind scheme for the first order derivatives. Because the paramter

is positive, concentration travels downwards at a speed

, I'm implementing the upwind scheme as:
$\frac{c^{i+1}_n - c^{i}_n}{\Delta z}$
for the node i. I also use the same upwind scheme for all the first-order derivative terms, i.e. also the terms inside the diffusion terms, which are nonlinear. For the second order partial derivatives, I apply second order central difference. For the time integration, I use a Crank-Nicholson method treating all terms implicitly.

The boundary conditions are complex. I'm using a standard wall function, using the log law at a node far from the wall. Using the log law, I also impose a gradient boundary condition because the log law adds the shear velocity as a unknown, hence the additional equation/condition. For

and $\varepsilon$ I use conditions that depend on the shear velocity. That's for one side of the domain, at the other side, sufficiently far, I use zero value conditions.

I'm not seeing much success with this implementation, should I use different discretization schemes for the diffusion and the advective terms? Should I use a different method altogether such as Finite Volume? What confuses me are the diffusion terms. Which discretization should I use there?

With regards to the time discretization. Should I treat some terms explicitely? The source terms are nonlinear as well.

Thanks in advance.

EDIT: Response to the first comment

Sorry I was not clear there. The model blows up for certain initial values of the concentration (too much acceleration), although that might have to do with the wall function implementation I think. Another issue that I had were the initial conditions for k and ε. I'll explain. The initial conditions for the velocity and concentration do not span the entire domain. They have positive values only for one tenth of the lower part of the domain, where I implement the wall function boundary condition, eventually, they will diffuse upwards. However, they will not do so if the values for k and ε are zero above where u and c are. This is because they are in the diffusion term and I think because of the upwind scheme, which assumes the wave travels downwards.

FMDenaro · July 19, 2015, 13:18

what about the problem in your results?
I don't think is a good idea to use upwind in the term inside the diffusion fluxes

July 19, 2015, 11:48	One dimensional k epsilon implementation from scratch	#1
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	Hello, I am trying to implement a 1D k- $\varepsilon$ turbulent model for a turbidity current, hence the conservation equation for $c$. I'm solving for the variables $u,c,k$ and $\varepsilon$. The remaining terms are parameters. $\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + C_{\mu} \frac{k^2}{\varepsilon} \right) \frac{\partial u }{\partial z} \right] + \frac{\rho_s - \rho_w}{\rho_w} c g_x$ $\frac{\partial c}{\partial t} = v_s \frac{\partial c}{\partial z} + \frac{\partial}{\partial z} \left[ \left( \nu + C_{\mu c} \frac{k^2}{\varepsilon} \right) \frac{\partial c }{\partial z} \right]$ $\frac{\partial k}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{C_{\mu} }{\sigma_k} \frac{k^2}{\varepsilon} \right) \frac{\partial k }{\partial z} \right] + C_{\mu} \frac{k^2}{\varepsilon} \left(\frac{\partial u }{\partial z}\right)^2 - \varepsilon + \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} C_{\mu} \frac{k^2}{\varepsilon} \frac{\partial c }{\partial z}$ $\frac{\partial \varepsilon}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{C_{\mu} }{\sigma_{\varepsilon}} \frac{k^2}{\varepsilon} \right)\frac{\partial \varepsilon }{\partial z}\right] + c_{\varepsilon 1} \frac{\varepsilon}{k}$ $\left( C_{\mu} \frac{k^2}{\varepsilon} \left(\frac{\partial u }{\partial z}\right)^2 + c_{\varepsilon 3} \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} C_{\mu} \frac{k^2}{\varepsilon} \frac{\partial c }{\partial z} \right) - C_{\varepsilon 2} \frac{ \varepsilon^2}{k}$ So far, I am using finite difference. Upwind scheme for the first order derivatives. Because the paramter is positive, concentration travels downwards at a speed , I'm implementing the upwind scheme as: $\frac{c^{i+1}_n - c^{i}_n}{\Delta z}$ for the node i. I also use the same upwind scheme for all the first-order derivative terms, i.e. also the terms inside the diffusion terms, which are nonlinear. For the second order partial derivatives, I apply second order central difference. For the time integration, I use a Crank-Nicholson method treating all terms implicitly. The boundary conditions are complex. I'm using a standard wall function, using the log law at a node far from the wall. Using the log law, I also impose a gradient boundary condition because the log law adds the shear velocity as a unknown, hence the additional equation/condition. For and $\varepsilon$ I use conditions that depend on the shear velocity. That's for one side of the domain, at the other side, sufficiently far, I use zero value conditions. I'm not seeing much success with this implementation, should I use different discretization schemes for the diffusion and the advective terms? Should I use a different method altogether such as Finite Volume? What confuses me are the diffusion terms. Which discretization should I use there? With regards to the time discretization. Should I treat some terms explicitely? The source terms are nonlinear as well. Thanks in advance. EDIT: Response to the first comment Sorry I was not clear there. The model blows up for certain initial values of the concentration (too much acceleration), although that might have to do with the wall function implementation I think. Another issue that I had were the initial conditions for k and ε. I'll explain. The initial conditions for the velocity and concentration do not span the entire domain. They have positive values only for one tenth of the lower part of the domain, where I implement the wall function boundary condition, eventually, they will diffuse upwards. However, they will not do so if the values for k and ε are zero above where u and c are. This is because they are in the diffusion term and I think because of the upwind scheme, which assumes the wave travels downwards. Last edited by salazardetroya; July 19, 2015 at 13:22. Reason: add details

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	16	March 4, 2017 08:30
calculation stops after few time steps	sivakumar	OpenFOAM Running, Solving & CFD	7	March 17, 2013 06:37
Calculation of k and epsilon freezes	Nigirim	OpenFOAM Running, Solving & CFD	1	November 14, 2012 07:52
epsilon and K blowing up.	sivakumar	OpenFOAM Running, Solving & CFD	1	October 25, 2012 04:50
SimpleFoam k and epsilon bounded	nedved	OpenFOAM Running, Solving & CFD	1	November 25, 2008 20:21

July 19, 2015, 13:18		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 Rep Power: 71	what about the problem in your results? I don't think is a good idea to use upwind in the term inside the diffusion fluxes