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:

Temperature Distribution across the Plate

No comments:

Post a Comment