Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Hyperbolic Equation Solution by a First Order Upwind Approximation by MATLAB

Let's assume that the temperature distribution inside a pipe is governed by the following one dimensional unsteady hyperbolic partial-differential-equation:

Tt + u Tx = 0

Consider negligible heat diffusion and the initial velocity u is, 0.1 cm/s. The boundary conditions are described below:

T(x, 0) = 200 x;       0 ≤ x ≤ 0.5

T(x, 0) = 200 (1 – x);       0.5 ≤ x ≤ 1

Develop an algorithm in MATLAB to solve this problem using the first-order approximation in time and the first-order upwind approximation in space. Consider, Δx = 0.05 cm; Δt = 0.05. 

The following MATLAB program develops the first-order upwind method to solve the given one-dimensional unsteady hyperbolic convection equation.


% Hyperbolic equation - Convection equation: Upwind Method
close all;
clc;

% Input Properties from Problem
L=1.0; % (dimension in cm)
Lmax=5.0;
alpha=0.01;
u=0.1; % (dimension in cm/s)
Tt=10; % Simulation Time Period

% Mesh/Grid size
dx=0.05;
dt=0.05;

c=u*dt/dx; % Depends on dt values
d=alpha*((dt/dx)^2.0);

% Matrix Parameters

imax=(Lmax/dx)+1; % Array dimension along the width
tmax=(Tt/dt); % Array dimension along the height
ndim=imax-2;
tdim=tmax-2;

% Total no of iterations
iteration=10000;

% Convergence Criterion
tolerance=0.000001;

% Array and Variables
% a = Matrix of Coefficients A(i,j)
% b = Right Side Vector, b(i)
% x = Solution Vector, x(i)

% Initial Condition
f(:,1)=0.0;
for i=1:imax-2
    x(i)=(i-1)*dx;
   
    if x(i)<=1.0
        f(i,1)=0.0;
    end
    if x(i)>1.0 && x(i)<=1.5
        f(i,1)=200*(x(i)-1.0);
    end
    if x(i)>1.5 && x(i)<=2.0
        f(i,1)=200*(L-x(i)+1.0);
    end
    if x(i)>2.0
        f(i,1)=0;
    end
end

fprintf('imax = %f\n',imax);
fprintf('c=%f\n',c);

% Upwind Method

for t=1:tmax
    time(t+1)=(t)*dt;
    for i=2:imax-2
        % First Order Approximation
        f(i,t+1)=f(i,t)-c*(f(i,t)-f(i-1,t));
    end
end

hold on
plot(x,f(:,1))
plot(x,f(:,find(time==1)),'o-')
plot(x,f(:,find(time==5)),'o-')
plot(x,f(:,find(time==10)),'o-')

xlabel('x, cm')
ylabel('Temperature, C')
legend('t=0s','t=1s','t=5s','t=10s')
title ('c=0.1')


Program Output:


Temperature Distribution inside Pipe for Various c and t

Temperature Distribution inside Pipe for Various c and t

Temperature Distribution inside Pipe for Various c and t

Temperature Distribution inside Pipe for Various c and t

Temperature Distribution inside Pipe for Various c and t

Temperature Distribution inside Pipe for Various c and t

No comments:

Post a Comment