Tutorial 13 – Turbulent Flow over NACA0012

Posted on July 17, 2017 - Tutorials

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

Introduction and Problem Definition

This tutorial is a continuation of Tutorial 12 and it will be assumed that you are familiar with concepts described in the previous tutorial. Here, we are going to simulate turbulent flow around a NACA-0012 airfoil and introduce a yet another turbulence model referred to as Constant Intensity Turbulence Model (CITM), which is developed as a hybrid model which uses Van-Driest model close to the wall and in the freestream it assumes turbulence with a predefined intensity and length scale. Thus, it leaves the user space to adjust it carefully to his case.

Geometry of the computational domain is shown below along with the ids of Physical Lines. We will simulate flow at an angle of attack of 5 degrees and Reynolds number of 50 000 based on airfoil chord.

Turbulence Model

CITM (Constant Intensity Turbulence Model) is a model which combines the well known Van-Driest model (refer to Tutorial 12 for details and an example) and the idea of the VL (Chen-Xu) turbulence model used for HVAC applications (see reference [1]). Let us remind a couple of formulae from turbulence theory and modelling. The $k-\epsilon$ turbulence model determines turbulent viscosity by

    $$\nu_t = 0.09 \frac{k^2}{\epsilon}$$

Turbulent kinetic energy can be defined by means of turbulence intensity $t_i$ and flow velocity $U$ as

    $$k = \frac{3}{2}(t_i U)^2$$

Futhermore, knowing the turbulent length scale $L_t$ one can determine dissipation of turbulent kinetic energy as

    $$\epsilon = \frac{k^{3/2}}{L_t}$$

Combining the above equations yields the following formula for turbulent viscosity

    $$\nu_t = 0.1102 t_i U L_t$$

The CITM model uses distance to the closest wall as the turbulent length scale or the ‘Ltmax’ parameter specified by the user, if distance to the closest wall is greater than ‘Ltmax’. This manifests the fact that turbulent length scale should not grow unbounded in freestream. In the boundary layer the CITM model uses a blending function to switch from the assumption of constant turbulence intensity and length scale to Van-Driest model which proved to be appropriate in proximity of walls. The Van-Driest model used close to the wall is defined 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 provide the script which simulates turbulent flow around a NACA-0012 airfoil. 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('naca0012.msh');


% Add boundary layer mesh (refer to Tutorial 5 for boundary layer meshing)
layers.h1 = 0.001;
layers.q = 1.3;
layers.n = 5;
wallIds = [9];
[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 (air kinematic viscosity in m^2/s)
nu = 1.5e-5;


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


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


% Define inlet velocity profile
vinf = 0.75; % 0.75 m/s yields Re = 50 000 at airfoil chord of 1 m.
aoa = 5; % angle of attack in degrees
aoa = aoa/180*pi; % calculate angle of attack in radians
vel = vinf*[cos(aoa) sin(aoa)]; % define velocity vector


% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute turbulent viscosity (turbulent intensity of 4% and Lt = 0.5m)
    nut = turbulentViscosity(p,t,u,d,nu,indices,'CITM', 0.04, 0.5);

    % 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,7,'inlet', vel);
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,9,'wall');

    % 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
% Display x-velocity field
figure(3)
displaySolution2D(p,t,u(indices.indu),'x-velocity');


% Display pressure coefficient field
pressure = generatePressureData(u,p,t);


% Define Cp = 1 at the stagnation point (find maximum pressure value on
% airfoil surface)
airfoilNodes = extractNodeIdsOnEdges(e,9);
pmax = max(pressure(airfoilNodes));


% Compute pressure coefficient and display it
Cp = (pressure-pmax+0.5*norm(vel)^2)/(0.5*norm(vel)^2);
figure(4);
displaySolution2D(p,t,Cp,'Cp [-]');


% Plot distribution of pressure coefficient on airfoil surface
figure(5)
scatter(p(1,airfoilNodes),-Cp(airfoilNodes));
xlabel('x/c [-]');
ylabel('-Cp [-]');
title('Pressure coefficient distribution');
grid on;


% Compute force at the airfoil (Physical Line id = 9)
Force = computeForce(p,e,t,u,nu,9);


% Project the force to direction tangential and perpendicular to freestream
% velocity (drag and lift)
D = Force*[cos(aoa) sin(aoa)]';
L = Force*[-sin(aoa) cos(aoa)]';


% Compute drag and lift coeeficients by dividing the force by dynamic
% pressure (unit fluid density is assumed for incompressible flow
% formulation)
Cd = D/(0.5*norm(vel)^2)
Cl = L/(0.5*norm(vel)^2)


% Below we show the pressure coefficient and velocity field in close
% surrounding of the airfoil.

 

References

[1] Q. Chen, W. Xu, “A zero-equation turbulence model for indoor airflow simulation”, Energy and Buildings 28 (1998) 137-144.>

Tags