Tutorial 15 – Laminar Heat Exchanger

Posted on July 17, 2017 - First steps, Tutorials

Introduction and Problem Definition

This tutorial is intended to be a yet another example of simulation of laminar heat transfer. Please refer also to Tutorial 14 on laminar heat convection from a flat plate. These two tutorials are closely linked with each other and it might be beneficial for you to study both of them to get more insights.

Here, we will simulate a heat exchanger with the geometry and settings depicted in figure below.

We will cut out representative geometry to make the simulation less computationally intensive. By doing this, we end up with the following geometry.

The fluid flow problem is defined by incompressible Navier Stokes equations and the heat convection problem is given by following convection diffusion equation

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

with $T = 20$ on inlet (id = 22), $T = 50$ at pipes (id = 26, 27, 28, 29 and 30) and $\vec n \cdot \nabla T = 0$ elsewhere. In the above equation the following denotations hold:

  • $\lambda$ – fluid thermal conductivity [W/(mK)],
  • $\rho$ – fluid density [kg/m3],
  • $c_p$ – fluid heat capacity [J/(kgK)].

Futhermore, it should be noted that we neither account for any viscosity changes due to temperature change nor for buoyancy effects (natural convection). This makes the flow problem decoupled flow heat transfer problem and they can be solved one after another. The coupling is only one-directional. Fluid flow solution influences the heat transfer problem, but heat transfer does not affect fluid flow. We will exploit these observation in the following section.

Code

The code provided in this section solves the above described heat exchanger problem. In the next section we do simple postproccessing of results to get global characteristics of the designed heat exchanger.

% Clear MATLAB Command Window and Workspace
clc;
clear;


% Read and display mesh
[p,e,t] = importMeshGmsh('heatExchanger.msh');
displayMesh2D(p,t);


% Convert mesh to second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);


% Define fluid properties
nu = 1e-6; % kinematic viscosity
lambda = 0.6; % thermal conductivity
rho = 1000; % fluid density
cp = 4190; % heat capacity
k = lambda/(rho*cp); % resulting thermal diffusivity


% Init solution and convergence criteria
[u, convergence] = initSolution(p,t,[0.05 0],0);
maxres = 1e-9;
maxiter = 25;


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


% Iterate nonlinear terms
for iter = 1:maxiter
    % Assemble matrix and right hand side
    [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,22,'inlet', vel);
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[24 25],'slipAlongX');
    [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,[26 27 28 29 30],'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
    u = NS\F;
end


% Display solution
figure(3);
displaySolution2D(p,t,u(indices.indu),'x-velocity');


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


% Apply boundary conditions
Tinlet = 20;
Tpipe = 50;
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,22,'value',Tinlet);
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,[26 27 28 29 30],'value',Tpipe);


% Solve equations
T = D\F;


% Display temperature field
figure(4);
displaySolution2D(p,t,T,'Temperature field');

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

% Display solution
% First, we compute the heat flux density field with formula given above 
% and calculate the total heat flux from pipe with id 26 and pipe 27.


% Total heat flux on pipe 26
gradT = solutionGradient2D(p,t,T);
qx = -lambda*gradT(:,1);
qy = -lambda*gradT(:,2);
Qpipe26 = boundaryFlux2D(p,e,[qx; qy],26)


% Total heat flux on pipe 27
Qpipe27 = boundaryFlux2D(p,e,[qx; qy],27)

Now, let us plot the temperature profile at heat exchanger outflow. To do that we will use MATLAB’s ‘scatter’ function. Hence, we need to know positions of nodes lying on the outlet and their temperatures. We first figure out, which nodes (with what numbers) are lying on the outlet (Physical Line 23) by using $\tt{extractNodeIdsOnEdges}$ function. Then, we extract temperatures of these nodes with $\tt{T(outletNodes)}$ and their y-coordinates with $\tt{p(2,outletNodes)}$ (array p stores x-coordinates of all nodes in its first row and y-coordinates in its second row).

% Outlet temperature profile
outletNodes = extractNodeIdsOnEdges(e,23);
figure(5);
scatter(T(outletNodes),p(2,outletNodes));
grid on;
xlabel('Temperature [deg]');
ylabel('y coordinate [m]');
title('Outlet temperature profile');

Next, we compute the total force exerted by the fluid on each pipe and define drag coefficients by taking the first component of that force. Drag coefficient is given by

    $$C_D = \frac{D}{\frac{1}{2} \rho v_\infty^2 A}$$

where $A$ stands for pipe projected area (in our 2-D case it is its diameter).

% Drag coefficient
% In our flow simulation we assume rho = 1, so Cd = Drag/(0.5*A*v_in^2);
v_in = 0.05;
A = 0.002;
Drag26 = computeForce(p,e,t,u,nu,26);
Drag27 = computeForce(p,e,t,u,nu,27);
Cd26 = Drag26(1)/(0.5*A*v_in^2)
Cd27 = Drag27(1)/(0.5*A*v_in^2)

We can also determine total pressure drop using $\tt{boundaryIntegral2D}$ function, which can both compute integral of a quantitiy over a boundary and average value of that quantity. We are interested in the latter case. Please remember that pressure that is result of computation is always physical pressure divided by density (common definition in incompressible flow simulations), so we have to multiply pressure field by fluid density to get real, physical pressure defined in Pascals.

% Total pressure drop
pressure = generatePressureData(u,p,t);
[~,p_inlet] = boundaryIntegral2D(p,e,pressure,22);
[~,p_outlet] = boundaryIntegral2D(p,e,pressure,23);
deltap = rho*(p_inlet-p_outlet);
disp(['deltap = ' num2str(deltap) ' Pa'])

Eventually, we are interested in temperature distribution on a vertical line in a middle of the domain. Obviously, no nodes are located there (or at least it is highly unusual for them to be aligned with any particular line). So, we first need to extract these data from the temperature field. Refer to $\tt{help}$ $\tt{extractDataAlongLine}$ for details of the subsequent arguments of the function, which does this job for us. After extraction of these data, the function puts temperature values into ‘Tl’ argument and coordinates of subsequent sampling point int ‘lp’ output arguments. Again, we use ‘scatter’ to plot the profile.

% Temperature distribution at x = 4e-3
[Tl,lp,le] = extractDataAlongLine(p,t,T,1e-3*[4 4],1e-3*[1. 3],100,50);
figure(6)
scatter(Tl,lp(2,:));
grid on;
title('Temperature distribution at x = 0.004 m');
xlabel('Temperature [deg]');
ylabel('y [m]');

Tags