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.

**T**_{t}+ uT_{x}= αT_{xx}*α*is 0.01

*cm*

^{2}

*/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:__

__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:__

## No comments:

## Post a Comment