## 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 ...

### Backward Time Centered Space Approach to Solve a Partial Differential Equation

Let us assume that there is a temperature distribution across a rectangular plate which is defined by the following unsteady one-dimensional convection-diffusion equation.

Tt + uTx = αTxx

Inside the plate, there is a flow of a fluid whose initial velocity is zero. The thickness of the plate is 1 cm, and the thermal diffusivity, α is 0.01 cm2/s. The initial temperature distribution is represented as,

T(x, 0) = 100 x / L;    0 ≤ x ≤ L

And, the boundary conditions are,

T(0, t) = 0;   T(L, t) = 100;

Assume, u = 0.1 cm/s. Develop a MATLAB program that solves this problem with the Backward Time Centered Space (BTCS) Approach. Consider, Δx = 0.1 and ΔT = 0.5. Plot the results for various times steps, such as, 0.5, 2.5, 5, 10, and 50 seconds. The Forward Time Centered Space (FTCS) approach has been shown here.

## Backward Time Centered Space (BTCS) Approach:

The following MATLAB program implements the BTCS approach to solve the unsteady one-dimensional convection-diffusion equation using the input parameters and boundary conditions.

% Backward Time Centered Space (BTCS) Approach

close all;
clc;

% Input Properties from Problem
L=1.0; % (dimension in cm)
alpha=0.01;
u=0.10;
Tt=50.0;  % Maximum simulation time
% Mesh/Grid size calculation
% Gird size definition along x and y axes
dx=0.1;
dt=0.5;
c=u*dt/dx;
d=alpha*(dt/dx^2);

% No of grid points
imax=(L/dx)+1; % Array dimension along the width
tmax=Tt/dt; % Array dimension along the height
ndim=imax-2;

% Coefficient Matrix Generation for Gauss-Seidel Method
for j=1:imax-2
for i=1:imax-2
if i==j
a(i,j)=1.0+2.0*d;end
if j==i+1
a(i,j)=c/2.0-d;end
if j==i-1
a(i,j)=-(c/2.0+d);end
end
end

for i=1:imax
x(i)=(i-1)*dx;
T(i,1)=100*x(i)/L;
end

for t=1:tmax
time(t+1)=t*dt;
for i=1:imax-2
b(i)=T(i+1,t);
end
b(1)=b(1)+(c/2.0+d)*T(1,t);
b(imax-2)=b(imax-2)-(c/2.0-d)*T(imax-2,t);

% Calling Gauss-Seidel function defined another 'm-file'
T(2:ndim+1,t+1)=GaussSeidelMethod(a,b);

T(1,t+1)=T(1,1);
T(imax,t+1)=T(imax,1);

end

hold on

plot(x,T(:,1),'LineWidth',2)
plot(x,T(:,2),'-p','LineWidth',2)
plot(x,T(:,6),'-*','LineWidth',2)
plot(x,T(:,11),'-s','LineWidth',2)
plot(x,T(:,21),'--','LineWidth',2)
plot(x,T(:,100),'^-','LineWidth',2)

xlabel('x, cm')
ylabel('Temperature, \circ C')
legend('t=0s','t=0.5s','t=2.5s','t=5s', 't=10s', 't=50s')

To solve the algebraic equations in BTCS method, the following MATLAB function develops Gauss-Seidel algorithm to solve the equations.

function x = GaussSeidelMethod(A,b,ET,maxiteration)
% x = GaussSeidelMethod(A,b): Gauss Seidel without relaxation
% Input:
% A = Coefficient Matrix
% b = Right Side Vector
% ET = Error Tolerance as Stop Criterion (default in program = 0.00001%)
% maxiteration = Maximum Iterations (default in program = 100)
% Output:
% x = Solution Vector

% The following logic determines if the given number of inputs are less
% than required
if nargin < 2
error('Minimum 2 Input Arguments Required'),
end

% The following condition checks the maximum number of iterations;
% If not given, this program considers maxiteration = 100
if nargin < 4 || isempty(maxiteration)
maxiteration = 100;
end
% The following loop is for error tolerance for solution;
% If none mentioned, this program considers it as 0.00001
if nargin<3||isempty(ET)
ET=0.00001;
end

% Gauss-Seidel Iteration Approach
[m,n] = size(A);
if m~=n
error('A Must be a Square Matrix');
end
C = A;

format short

for i = 1:n
C(i,i) = 0;
x(i) = 0;
end
x = x';

for i = 1:n
C(i,1:n) = C(i,1:n)/A(i,i);
end

for i = 1:n
d(i) = b(i)/A(i,i);
end
iteration = 0;

while (1)
xold = x;      % Value of x is being updated simultaneously
for i = 1:n
x(i) = d(i)-C(i,:)*x;
if x(i) ~= 0
r = abs(b-A*xold);   % Equation for Residual
end
end
iteration = iteration+1;     % Number of iterations updated

%End Loop Condition
if max(abs(b-A*xold)) <= ET | iteration >= maxiteration
break,
end
end

Program Output: