Tutorial 19 – Turbulent HVAC Simulation

Introduction and Problem Definition

This tutorial introduces into simulation of a benchmark HVAC (heating, ventilation, air-conditioning problem). The test case will be based on reference [1], where a 2-D geometry has been proposed along with detailed experimental results that allow us to compare simulation results with the experiment. Quasi 2-D geometry has been reproduced from [1] in picture below.

We will be dealing in the benchmark with ventilation case of a model room with one inlet, through each air enters with velocity of 0.4554 m/s. One outflow with zero average pressure is located in the bottom right corner of the room.

Actual 2-D geometry, which will be taken for simulation is shown below with all Physical Line ids that will be important for specification of boundary conditions.

Turbulence Modelling

The flow under investigation is a turbulent flow. Thus, we will generate mesh for boundary layer in this tutorial (please refer to Tutorial 5 on boundary layer meshing for details of parameters) and we will exploit one of the turbulence models implemented in QuickerSim CFD Toolbox. A fast turbulence model developed and calibrated for HVAC applications which is implemented in QuickerSim CFD Toolbox is that proposed by Chen and Xu and named as ‘ChenXu’ model in this toolbox. Please refer to [2] – original paper, which introduces the model.

RANS turbulence model require determination of turbulent viscosity, which in case of ‘ChenXu’ model is given by explicit formula:

where stands for local velocity magnitude and for distance to the closest wall.

Code

In this section we provide the complete code which simulates the turbulent flow for the above mentioned HVAC benchmark.

% Clear MATLAB Command Window and Workspace
clc;
clear;
% Read mesh
[p,e,t] = importMeshGmsh('hvac.msh');


% Add boundary layer mesh (refer to Tutorial 5 for details)
layers.h1 = 0.005;
layers.q = 1.5;
layers.n = 3;
wallIds = [16];
[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; % air kinematic viscosity
% Init solution and define convergence and under-relaxation parameters

[u, convergence] = initSolution(p,t,[0 0],0);
maxres = 1e-4;
maxiter = 20;
alpha = 0.6;

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


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


% Iterate nonlinear terms
for iter = 1:maxiter
    % Compute turbulent viscosity
    nut = turbulentViscosity(p,t,u,d,nu,indices,'ChenXu');

    % 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,14,'inlet', vel);
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,16,'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 x-velocity field
figure(3)
displaySolution2D(p,t,u(indices.indu),'x-velocity');

Display solution

Now, we compare computed velocity profiles to the experimental ones. We will carry out out validation based on the x-velocity profiles which have been measured in experiment along two vertical lines at x = 3 m and x = 6 m.

First, we extract these data from the solution with with ‘extractDataAlongLine’ function. Please refer to command line help to this function for detailed description of its input and output arguments.

% Extract streamwise velocity profiles at x = 3 and x = 6
[vx3, p3] = extractDataAlongLine(p,t,u,[3 3],[-3 0],50,20);
[vx6, p6] = extractDataAlongLine(p,t,u,[6 6],[-3 0],50,21);


% Load experimental data
exp_vx3 = load('exp_vx3.txt');
exp_vx6 = load('exp_vx6.txt');


% Plot x-velocity profiles at x = 3 m
figure(4)
scatter(vx3,p3(2,:),'+')
hold on
scatter(exp_vx3(:,2),exp_vx3(:,1))
hold off
legend('Simulation','Experiment','Location','Best');
title('Validation of horizontal velocity profile at x = 3 [m]')
xlabel('x-velocity [m/s]')
ylabel('y [m]')
grid on;

% Plot x-velocity profiles at x = 6 m
figure(5)
scatter(vx6,p6(2,:),'+')
hold on
scatter(exp_vx6(:,2),exp_vx6(:,1))
hold off
legend('Simulation','Experiment','Location','Best');
title('Validation of horizontal velocity profile at x = 6 [m]')
xlabel('x-velocity [m/s]')
ylabel('y [m]')
grid on;

References

[1] P. V. Nielses, “Specification of a two-dimensional test case”, Instituttet for Bygningsteknik, Aalborg Universitetscenter, November 1990.

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