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.
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:
with on plate surface, at inlet and 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');
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
and the local Nusselt number defined by:
In the above stands for heat transfer coeeficient (defined below), for coordinate measured from the beginning of the flat plate. stands for difference of wall temperature and freestream fluid temperature. Heat transfer coefficient is defined by following relation:
Eventually, analytically derived expression for local Nusselt number in laminar heat convection from a flat plate is given by
where Prandtl number is given by
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 [-]');