Tutorial 24 – RANS Prandtl k-l Turbulence Model

Posted on July 17, 2017 - CFD Numerics, Turbulent Flow, Tutorials

This tutorial is intended for the full version of the toolbox.

Introduction

In this tutorial we are going to introduce a one-equation turbulence model, which is more general and accurate than simple zero-equation algebraic turbulence models. The model falls within the class of RANS turbulence models and models transport of turbulence kinetic energy. The second turbulent quantity which needs to be specified to compute turbulent viscosity is after Prandtl assumed to be length scale of turbulence which in this model is equal to the distance to the closest wall.

Turbulence Model

Mathematical formulation of the Prandtl k-l turbulence model is given below. The RANS equations read:

    $$\vec{u}\cdot\nabla\vec{u}=-\nabla p+\nabla\cdot((\nu+\nu_{t})(\nabla\vec{u}+(\nabla\vec{u})^{T}))$$

In the above $nu_{t}$ stands for turbulent kinetic energy, which is defined in the present model by

    $$\nu_{t}=C_{\mu}k^{\frac{1}{2}}l$$

where a model constant $C_{\mu} = 0:09$ , $k$ stands for turbulence kinetic energy and $l$ for distance to the closest wall. The $k$ equation in turn reads

    $$\vec{u}\nabla k-\nabla\cdot((\nu+\frac{\nu}{\sigma_{k}})\nabla k)=\nu_{t}S^{2}-\frac{k^{\frac{3}{2}}}{l_{mix}}$$

where $l_{mix}=l$ and $S=\sqrt{2S_{ij}S_{ij}}$ , where $S_{ij}=0.5(\nabla\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})$

Geometry, Mesh and Boundary Conditions

In the present case we are going to solve turbulent backward facing step flow with Re = 50000 based on inlet velocity and inlet height. Furthermore, we need to specify turbulence intensity level as fraction of unity ( $t_{i} = 0:04$ denotes turbulence intensity level of 4). All other settings of the solver like convergence settings or boundary layer mesh generation settings should already be known to the user. Otherwise, please refer to introductory tutorials on QuickerSim CFD Toolbox and tutorials on boundary layer mesh generation. The respective geometry with ids of physical lines is depicted in figure below.

The mesh is displayed by the solver during simulation run and is shown farther in the tutorial. One more comment should be made on the boundary conditions for the turbulence kinetic energy transport equation. The model requires that each inlet has specified values of velocity and turbulence intensity. The solver calculates Dirichlet boundary condition for k at given inlet through the well known formula in turbulence modeling:

    $$k=\frac{3}{2}(U\cdot t_{i})^{2}$$

Outlet assume homogenuous Neumann condition and kinetic turbulence energy at all material walls is identically zero due to the fact that turbulent fluctuations vanish on a physical wall.

Code

Below we show complete code that solves the above specified problem. Clear MATLAB Command Window and Workspace

clc;
clear;

% Read and display mesh
[p,e,t] = importMeshGmsh('bfs.msh');

% Add boundary layer mesh
layers.h1 = 0.015;
layers.q = 1.5;
layers.n = 3;
wallIds = [11 12];
[p,e,t] = extrudeLayers2D(p,e,t,wallIds,layers);


% Convert mesh to second order and display mesh
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);
displayMesh2D(p,t);


% Define fluid
nu = 1.5e-5;


% Init solution and define convergence and under-relaxation parameters
[u, convergence] = initSolution(p,t,[0 0],0);
maxres = 1e-4;
maxiter = 20;
alpha = 0.7;


% Init turbulence model
turbulenceModel = 'k-l';
ti = 0.04;
[d,us,k,M] = initTurbulenceModel(p,e,t,u,wallIds,layers,turbulenceModel,0,ti);


% Define inlet velocity profiles
vel = [0.75 0]; % (Re = 50 000)
inletIds = [9];
turbData = [norm(vel) ti];


% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute turbulent viscosity (take care of additional arguments
    % depending on turbulence model)
    [nut,yplus] = turbulentViscosity(p,t,u,d,nu,indices,turbulenceModel,k);

    % Under-relax velocity and assemble matrix and right hand side
    u = alpha*u+(1-alpha)*us;
    % [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu+nut,u(indices.indu),u(indices.indv),'nosupg');
    [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu+nut,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic');

    % Apply boundary conditions
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,9,'inlet', vel);
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[11 12],'wall');
    % [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[13],'slipAlongX');

    % Compute and plot residuals
    [stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres);
    plotResiduals(convergence,2);

    % Break if solution converged
    if(stop)
        break;
    end

    % Solve equations
    us = u;
    u = NS\F;

    k = solveTurbulence(p,e,t,u,nut,k,d,nu,M,indices,wallIds,inletIds,turbulenceModel,turbData);
end


% Display Solution

figure(5)
displaySolution2D(p,t,u(indices.indu),'x-velocity')


figure(6)
displaySolution2D(p,t,k,'Turbulence kinetic energy [m2/s2]');


figure(7)
displaySolution2D(p,t,nut/nu,'Turbulent viscosity ratio');


figure(8)
ti = 100*sqrt(2/3*k)/norm(vel);
displaySolution2D(p,t,ti,'Turbulence intensity [%]');

Tags