Tutorial 12 – Turbulent Flow over a Plate

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

Introduction and Problem Definition

In this tutorial we are going to cover simulation of a turbulent flow over a flat plate. We will apply one of the algebraic, zero-equation turbulence models available in QuickerSim CFD Toolbox.

Geometry of the computational domain is shown below along with the ids of Physical Lines, where we will be specifying boundary conditions.

We will assume no shear at Physical Line 11, inlet at Line 8, outflow at Line 9, no shear slip at Line 10 and wall (flat plate) at Line 12. The length of the plate is 4 meters, which for inlet velocity of 10 m/s and kinematic viscosity of 1.5e-5 m^2/s yields Reynolds number based on plate length equal to 2.67e6.

We will study this flow with Van-Driest model (which is improvement of the Prandtl mixing length model). Both models are implemented in QuickerSim CFD Toolbox and the Van-Driest model is given by:

    $$\nu_t = l_{mix}^2 | \frac{\partial u}{\partial y}|$$

with

    $$l_{mix} = \kappa y (1-e^{\frac{-y^+}{A_0^+}})$$

where $\kappa = 0.41$, $A_0^+ = 26$, $y$ stands for wall distance and

    $$y^+ = \frac{u_\tau y}{\nu}$$

where $u_\tau$ is the friction velocity defined as

    $$u_\tau = \sqrt{\nu \left|\frac{\partial u}{\partial y}\right| |_{wall}}$$

Code

In this section we present the whole script to simulate a turbulent flow in QuickerSim CFD Toolbox. In contrary to laminar flow simulation we will have a couple of additional elements involved:

  • computation of wall distance,
  • definition of under-relaxation parameter, which enhances stability of convergence process.

The majority of description is provided in comments to the source code – details of subsequent functions can be found in documentation.

% Clear MATLAB Command Window and Workspace

clc;
clear;

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

% Add boundary layer mesh (for details on boundary layer meshing refer to
% Tutorial 5)
layers.h1 = 0.0001;
layers.q = 1.7;
layers.n = 12;
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-3;
maxiter = 40;
alpha = 0.8;


% choosing a lower value of under-relaxation parameter should
% prevent big oscilations in convergence process


% Compute wall distance and duplicate velocity vector
d = computeWallDistance(p,e,t,wallIds);
us = u;


% Define inlet velocity profile
vel = [10 0];


% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute turbulent viscosity (you can change 'VanDriest' to 'Prandtl')
    [nut, yplus] = turbulentViscosity(p,t,u,d,nu,indices,'VanDriest');

    % 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),'supgDoublyAsymptotic');

    % Apply boundary conditions
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,8,'inlet', vel);
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,12,'wall');
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[10 11],'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;
end

Display solution

In this section we are going to display the results and compare computed outlet velocity profile with the reference solution given by:

    $$w = w_\infty(\frac{y}{\delta})^{1/7}$$

where $\delta$ – boundary layer thickness can be calculated from

    $$\delta(x) = 0.376 (\frac{\nu}{w_\infty})^{0.2} x^{0.8}$$

% Display field of streamwise velocity component
figure(4)
displaySolution2D(p,t,u(indices.indu),'x-velocity');

% Dispay y plus field
figure(6)
displaySolution2D(p,t,yplus,'y+');

% Compare computed velocity profile with analytical one
% Extract indices of outlet nodes
outletNodes = extractNodeIdsOnEdges(e,9);


% Compute boundary layer thickness and reference velocity profile
deltaout = 0.376*(nu/norm(vel))^0.2*(4-1).^(0.8);
wout = norm(vel)*(p(2,outletNodes)/deltaout).^(1/7);


% Plot reference and computed profiles to the figure
figure(5)
scatter(wout,p(2,outletNodes),'+');
hold on;
scatter(u(indices.indu(outletNodes)),p(2,outletNodes))
hold off;
grid on;
ylim([0 deltaout])
xlabel('w [m/s]');
ylabel('y [m]');
title('Turbulent boundary layer profile')
legend('Reference solution', 'Numerical solution');