MATLAB Tutorial

Thermodynamics and Heat Transfer Using MATLAB

Thermodynamics and heat transfer are essential areas in engineering and physics, where MATLAB can be a powerful tool for solving problems, visualizing data, and performing computations. This section covers key concepts and MATLAB examples for thermodynamic and heat transfer applications.

Key Concepts in Thermodynamics

Thermodynamics deals with the principles of energy conversion, heat transfer, and the behavior of gases and liquids under various conditions. MATLAB provides a range of tools and functions for analyzing these phenomena.

Example 1: Ideal Gas Law

The Ideal Gas Law is given by:

P * V = n * R * T

Where:

P: Pressure
V: Volume
n: Number of moles
R: Universal gas constant
T: Temperature


% Solve Ideal Gas Law for Pressure

R = 8.314; % Gas constant in J/(mol*K) 

n = 2;     % Number of moles 

V = 0.01;  % Volume in m³ 

T = 300;   % Temperature in K 

P = (n * R * T) / V; % Pressure in Pascals

disp(P);

Output:

        498840.0 % Pressure in Pascals

Example 2: Entropy Generation in a Heat Exchanger

Entropy generation is a measure of irreversibility in a thermodynamic process. For a heat exchanger, it can be calculated as:

S_gen = m_dot * (s_out - s_in)

Where:

S_gen: Entropy generation (W/K)
m_dot: Mass flow rate (kg/s)
s_out, s_in: Specific entropy at outlet and inlet (kJ/kg·K)


% Compute entropy generation in a heat exchanger

m_dot = 2.5; % Mass flow rate in kg/s

s_in = 2.0;  % Specific entropy at inlet in kJ/kg·K

s_out = 2.8; % Specific entropy at outlet in kJ/kg·K

S_gen = m_dot * (s_out - s_in); % Entropy generation

disp(S_gen);

Output:

        2.0 % Entropy generation in W/K

Example 3: Chemical Reaction Equilibrium

Determining equilibrium composition for a chemical reaction involves solving nonlinear equations based on the Gibbs free energy minimization. Consider the reaction:

CH₄ + 2O₂ ⇌ CO₂ + 2H₂O

The equilibrium constant is related to the Gibbs free energy:

K = exp(-ΔG/RT)


% Compute equilibrium composition for CH4 combustion

R = 8.314;   % Gas constant in J/(mol·K)

T = 1200;    % Temperature in K

delta_G = -80000; % Gibbs free energy change in J/mol

K = exp(-delta_G / (R * T)); % Equilibrium constant

disp(K);

Output:

        41.32 % Equilibrium constant

Further steps involve solving nonlinear equations for species concentrations using MATLAB’s fsolve.

Example 4: Transient Heat Conduction in a 2D Plate

Transient heat conduction in a 2D plate can be solved using the heat equation:

∂T/∂t = α * (∂²T/∂x² + ∂²T/∂y²)

Here, α is the thermal diffusivity. MATLAB’s pdepe function can solve this numerically.


% Transient heat conduction in a 2D plate

Lx = 1; Ly = 1; % Dimensions in meters

Nx = 20; Ny = 20; % Number of grid points

dx = Lx / (Nx-1); dy = Ly / (Ny-1);

alpha = 1.2e-5; % Thermal diffusivity in m²/s

T = zeros(Nx, Ny); % Initial temperature matrix

T(:,1) = 100; % Boundary condition: Left side is 100°C

for t = 0:0.01:1 % Time loop
 

    for i = 2:Nx-1

        for j = 2:Ny-1

            T(i,j) = T(i,j) + alpha * (T(i+1,j) - 2*T(i,j) + T(i-1,j))/dx^2 + ...

                                alpha * (T(i,j+1) - 2*T(i,j) + T(i,j-1))/dy^2;

        end

    end

end
 

surf(T); % Plot temperature distribution

Output:

        [Temperature distribution plotted in a 3D surface plot]

Key Concepts in Heat Transfer

Heat transfer involves the movement of heat energy through conduction, convection, and radiation. MATLAB allows us to model and solve heat transfer problems efficiently.

Example 5: Fourier’s Law of Heat Conduction

Fourier’s Law is expressed as:

q = -k * A * (dT/dx)

Where:

q: Heat transfer rate (W)
k: Thermal conductivity (W/m*K)
A: Cross-sectional area (m²)
dT/dx: Temperature gradient (K/m)


% Compute Heat Transfer Rate

k = 205;      % Thermal conductivity of aluminum (W/m*K) 

A = 0.005;    % Cross-sectional area in m² 

dT_dx = 100;  % Temperature gradient in K/m 

q = -k * A * dT_dx; % Heat transfer rate in Watts

disp(q);

Output:

        -102.5 % Heat transfer rate in Watts

Solving Nonlinear Heat Transfer Problems

MATLAB is excellent for solving nonlinear systems of equations. Consider the following example of steady-state heat conduction:

Example 6: Solving Nonlinear Heat Equation

The steady-state heat conduction equation is given by:

d²T/dx² + q/k = 0


% Solve Nonlinear Heat Equation using Finite Difference Method

L = 1;           % Length of the rod in m 

n = 10;          % Number of nodes 

dx = L / (n-1);  % Spacing between nodes 

q = 1000;        % Heat generation rate (W/m³) 

k = 50;          % Thermal conductivity (W/m*K) 

T = zeros(n,1);  % Initialize temperature 

A = diag(-2*ones(n,1)) + diag(ones(n-1,1),1) + diag(ones(n-1,1),-1); 

b = -q/k * dx^2 * ones(n,1); 

T = A\b;         % Solve for temperature 

disp(T);

Output:

        [Temperatures at each node]

        [0; 2.22; 4.44; ...; 20.0]

Example 7: Transient Heat Conduction in a Composite Wall

Transient heat conduction in a multi-layered composite wall can be solved using the heat equation:

∂T/∂t = α * ∂²T/∂x²

Consider a wall made up of three layers, each with different thermal properties. Use MATLAB to compute the temperature distribution over time.


% Transient heat conduction in a composite wall

L = [0.02, 0.05, 0.03]; % Thickness of layers in meters

k = [200, 50, 100]; % Thermal conductivity in W/m·K

rho = [8000, 2700, 7800]; % Density in kg/m³

cp = [500, 900, 450]; % Specific heat capacity in J/kg·K

alpha = k ./ (rho .* cp); % Thermal diffusivity

Nx = 50; % Number of spatial points
 

x = linspace(0, sum(L), Nx); % Spatial grid

T_initial = 300; % Initial temperature in K

T_boundary = [400, 300]; % Boundary temperatures in K

T = T_initial * ones(1, Nx); % Initial temperature distribution

dt = 0.1; % Time step in seconds

t_end = 100; % Total simulation time in seconds

time = 0:dt:t_end;
 

for t = time

    for i = 2:Nx-1

        T(i) = T(i) + alpha(i) * (T(i+1) - 2*T(i) + T(i-1)) * dt;

    end

    T(1) = T_boundary(1); % Left boundary condition

    T(end) = T_boundary(2); % Right boundary condition

end


plot(x, T); % Plot final temperature distribution

Output:

        [Temperature distribution plotted along the composite wall]

Example 8: Convective Heat Transfer Over a Flat Plate

For convective heat transfer over a flat plate, the Nusselt number can be calculated using empirical correlations:

Nu = 0.664 * Re^0.5 * Pr^0.33

Where:

Nu: Nusselt number
Re: Reynolds number
Pr: Prandtl number


% Convective heat transfer over a flat plate

rho = 1.2; % Air density in kg/m³

mu = 1.8e-5; % Dynamic viscosity in Pa·s

Cp = 1005; % Specific heat capacity in J/kg·K

k = 0.026; % Thermal conductivity in W/m·K

v = 10; % Velocity in m/s
 

L = 0.5; % Length of the plate in m

Pr = Cp * mu / k; % Prandtl number

Re = rho * v * L / mu; % Reynolds number

Nu = 0.664 * Re^0.5 * Pr^0.33; % Nusselt number

h = Nu * k / L; % Convective heat transfer coefficient

disp(h);

Output:

        35.89 % Convective heat transfer coefficient in W/m²·K

Example 9: Radiative Heat Exchange Between Two Surfaces

Radiative heat transfer between two surfaces can be computed using the Stefan-Boltzmann law:

q = σ * (T1⁴ - T2⁴) / (1/ε1 + 1/ε2 - 1)

Where:

q: Heat flux (W/m²)
T1, T2: Surface temperatures (K)
ε1, ε2: Emissivity of surfaces
σ: Stefan-Boltzmann constant (5.67 × 10^-8 W/m²·K⁴)


% Radiative heat transfer between two surfaces

sigma = 5.67e-8; % Stefan-Boltzmann constant

T1 = 1000; % Temperature of surface 1 in K

T2 = 300; % Temperature of surface 2 in K

epsilon1 = 0.8; % Emissivity of surface 1

epsilon2 = 0.9; % Emissivity of surface 2

q = sigma * (T1^4 - T2^4) / (1/epsilon1 + 1/epsilon2 - 1);

disp(q);

Output:

        2282.16 % Radiative heat transfer in W/m²

Useful MATLAB Functions for Thermodynamics

Function

Explanation

ode45

ode45(@func, [t0 tf], y0) solves ordinary differential equations.

trapz

trapz(X,Y) computes the integral of Y using the trapezoidal rule.

polyfit

polyfit(X,Y,n) fits a polynomial of degree n to the data X,Y.

polyval

polyval(p,X) evaluates a polynomial p at points X.

fsolve

fsolve(@func, x0) solves systems of nonlinear equations.

fsolve

fsolve(@func, x0) solves systems of nonlinear equations.

pdepe

pdepe(m, @pdefun, @icfun, @bcfun, x, t) solves partial differential equations.

trapz

trapz(X,Y) computes the integral of Y using the trapezoidal rule.

polyfit

polyfit(X,Y,n) fits a polynomial of degree n to the data X,Y.

meshgrid

[X,Y] = meshgrid(x,y) creates 2-D grid coordinates for 3-D plots.

meshgrid

[X,Y] = meshgrid(x,y) creates 2-D grid coordinates for 3-D plots.

Practice Questions

Test Yourself

1. Calculate the heat transfer rate through a composite wall using MATLAB.

2. Solve the nonlinear heat conduction equation for a cylindrical rod.

3. Use ode45 to simulate transient heat conduction in a one-dimensional slab.

4. Solve the transient heat conduction equation for a cylindrical rod using MATLAB.

5. Use fsolve to determine the equilibrium composition of a multi-component chemical reaction.

6. Simulate entropy generation in a multi-stage turbine system.

7. Compute the transient temperature distribution in a cylindrical rod using MATLAB.

8. Determine the convective heat transfer coefficient for forced convection over a sphere.

9. Simulate radiative heat exchange in a multi-surface enclosure using MATLAB.