Tutorial 14 – Laminar Heat Convection over a Plate

Posted on July 17, 2017 - Tutorials

Introduction and Problem Definition

In this tutorial we are going to simulate laminar heat convection from a flat plate. 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, so we will specify inlet velocity of 0.0075 m/s to obtain Reynolds number of 2000 based on plate length. For heat transfer, we will assume zero temperature at inlet, equal to one on the plate surface and no conductive flux elsewhere.

Code

Below we present the whole computational script, which performs convective heat transfer simulation in QuickerSim CFD Toolbox. In the last section, we compare numerical results to analytical solutions.

% Clear MATLAB Command Window and Workspace
clc;
clear;


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


% Add boundary layer mesh (refer to Tutorial 5 for boundary layer meshing)
layers.h1 = 0.01;
layers.q = 1.3;
layers.n = 5;
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; % air kinematic viscosity (m^2/s)


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


% Define inlet velocity profile
vel = [0.0075 0]; % Re = 2000 based on plate length

% Iterate nonlinear terms
for iter = 1:maxiter
    % Assemble Navier-Stokes matrix (without or with stabilization)
    [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,u(indices.indu),u(indices.indv),'nosupg');
    % [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,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
    u = NS\F;
end

Solving Heat Transfer Problem

Solution of the heat transfer problem does not influence any properties of the flow field. That is why, we can decouple these two problems and solve them one after another. This would not be the case, if we wanted to incorporate calculation of natural convection or if temperature changes were big enough to modify fluid viscosity. Below, we define thermal properties of the fluid, which result in Prandtl number of 0.7 (ratio of kinematic viscosity to thermal diffusivity). The heat transfer problem is defined by the following equation:

    $$\vec{u} \cdot \nabla T = \nabla \cdot ( \frac{\lambda}{\rho c_p} \nabla T)$$

with $T = 1$ on plate surface, $T = 0$ at inlet and $\vec n \cdot \nabla T = 0$ elsewhere.

% Define thermal fluid properties
lambda = 0.0257;
rho = 1.205;
cp = 1005;
k = lambda/(rho*cp);
Pr = nu/k;


% Assemble problem matrix (convection diffusion equation)
[D,F] = assembleDiffusionMatrix2D(p,t,k);
C = assembleScalarConvectionMatrix2D(p,t,k,u(indices.indu),u(indices.indv),'nosupg');
D = D + C;

% Apply boundary conditions
Tin = 0;
Tw = 1;
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,8,'value',Tin);
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,12,'value',Tw);


% Solve equations
T = D\F;


% Display solution
figure(3);
displaySolution2D(p,t,T,'Temperature field');

Display solution

We will now verify simulation results against analytical solution. To do that, we need to get to know, which nodes (in the sense of their ids) are lying on plate surface and we will compute distribution of heat flux density from the plate

    $$\vec q = -\lambda \nabla T$$

and the local Nusselt number defined by:

    $$Nu_x = \frac{\alpha x}{\Delta T}$$

In the above $\alpha$ stands for heat transfer coeeficient (defined below), $x$ for coordinate measured from the beginning of the flat plate. $\Delta T = T_{wall} - T_{\infty}$ stands for difference of wall temperature and freestream fluid temperature. Heat transfer coefficient is defined by following relation:

    $$|\vec q| = \alpha (T_{wall} - T_{\infty})$$

Eventually, analytically derived expression for local Nusselt number in laminar heat convection from a flat plate is given by

    $$Nu_{ref} = 0.332 Pr^{1/3} \sqrt{\frac{U_\infty x}{\nu}}$$

where Prandtl number is given by

    $$Pr = \frac{\nu}{k}$$

with

    $$k = \frac{\lambda}{\rho c_p}$$

The code below computes all quantities described above and plots comparison of the simulated and reference Nusselt number distribution.

% Extract ids of nodes on the flat plate
wallNodes = extractNodeIdsOnEdges(e,12);


% Compute temperature gradient and heat flux denstiy
gradT = solutionGradient2D(p,t,T);
htc = -lambda*gradT(:,2)/(Tw-Tin);


% Compute numerical Nusselt number
Nu = htc.*(p(1,:)'-1)/lambda;


% Compute reference Nusselt number
Nuref = 0.332*Pr^(1/3)*sqrt(norm(vel)*(p(1,wallNodes)-1)/nu);


% Plot figure
figure(4)
hold on
scatter(p(1,wallNodes),Nuref,'bo');
scatter(p(1,wallNodes),Nu(wallNodes),'r+');
legend('Reference solution','Numerical solution');
hold off;
grid on
title('Comparison of computed and reference Nusselt number');
xlabel('x [m]');
ylabel('Nu [-]');

Tags