# 1D k-epsilon turbulence model in a turbidity current

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. $$\newcommand{\cmu}{c_{\mu}} \newcommand{\cmuc}{c_{\mu c}} \newcommand{\dervz}{\frac{\partial #1}{\partial z}} \newcommand{\eddyvis}{\cmu \frac{k^2}{\varepsilon}} \newcommand{\eddyvisC}{\cmuc \frac{\hat{k}^2}{\hat{\varepsilon}}}$$

$$\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \cmu \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 \dervz{c} + \frac{\partial}{\partial z} \left[ \left( \nu + \cmuc \frac{k^2}{\varepsilon} \right) \dervz{c} \right]$$

$$\frac{\partial k}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{\cmu}{\sigma_k} \frac{k^2}{\varepsilon} \right) \dervz{k} \right] + \eddyvis \left(\dervz{U}\right)^2 - \varepsilon + \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c}$$

$$\frac{\partial \varepsilon}{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \frac{\cmu}{\sigma_{\varepsilon}} \frac{k^2}{\varepsilon} \right) \dervz{\varepsilon} \right] + c_{\varepsilon 1} \frac{\varepsilon}{k} \left( \eddyvis \left(\dervz{u}\right)^2 + c_{\varepsilon 3} \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c} \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 $v_s$ is positive, concentration travels downwards at a speed $v_s$, 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 linear diffusion terms, 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 $k$ 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.

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 $\varepsilon$. 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 $\varepsilon$ 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.
• What are you out for? Since turbulence is a 3D phenomenon, I don't think that $\kappa-\epsilon$ models make any sense in 1D. – Jan Jul 19 '15 at 18:18
• I don't think this is enough for turbulence (models) to have a velocity gradient only in streamwise direction. E.g., the basic models for turbulence, like the $\kappa-\epsilon$ model, base on Prandtl's mixing length model, that bases on $\partial \bar u /\partial y$, i.e., the gradient of the mean flow in transversal direction. – Jan Jul 20 '15 at 6:34
• There are no velocity component in the transversal direction, but there are velocity gradients. I typed them as $\frac{\partial u}{\partial z}$. – balborian Jul 20 '15 at 13:59