This tutorial is intended for the full version of the toolbox.
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.
Mathematical formulation of the Prandtl k-l turbulence model is given below. The RANS equations read:
In the above stands for turbulent kinetic energy, which is defined in the present model by
where a model constant , stands for turbulence kinetic energy and for distance to the closest wall. The equation in turn reads
where and , where
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 ( 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:
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.
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 = ; 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,,'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 [%]');