2 Dimensional heat equation using ADI method

Question

0 个投票

I have made the code using TDMI line by line for my problem. The code is taking too long to run and is not producing the results, also not showing any errors or warnings. I would be grateful if someone could help.

A slab with length 50 mm, height 1mm and width 1 mm. It is a stainless steel slab, having the temperature at the bottom= 90C(363.15l). The left and right sides of the salb are insulated and the top side has a flowing fluid of temperature equals to 40C(313.15k). I have used TDMI line by line.

%convective boundary layer(Tinf = 40C = 313.15K)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%i i%

%n n%

%s s%

%u u%

%1mm%l l%

%a a%

%t t%

%e e%

%d d%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Constant Temperature wall (T0 = 90C = 363.15K)

%50mm

% converting the 2d heat convection equation in the for,:

% ap*Tp = aw*Tw + ae*Te + an*Tn + as*Ts +bp , where ap, aw, as, an and ae are the coefficients of the temperatures %at point p , west, south, north and east respectively to the node selected. Point p is the point of our concern.

% By applying the energy equation for different nodes we can get the

% respective value values of aw, ae, an, as and bp

%To apply TDMA we need to have an equation in the form:

% T(i) = P(i)*T(i+1) + R(i)

%The above equation can also be written as T(i-1) = P(i-1)*T(i) + R(i-1)

% now after putting the value of T(i-1) in

% ap*Tp = aw*Tw + an*Tn + as*Ts+ ae * Te + bp

% we get an equation which gives us the value of P(i) and R(i).

% P(i,j) = ae/ (ap(i,j) - aw*P(i-1,j));

% R(i,j) = (aw*R(i-1,j) + bp**) / (ap(i,j) - aw*P(i-1,j));

% T(i,j) = P(i,j)*T(i+1,j) + R(i,j)

% for i = 1

%P(1,j) = ae/ap(i,j)

%R(i,j) = bp**/ap(i,j)

% for i = nx

% T(nx,j) = R(nx,j);

CODE

% Full Implicit Form

%defining initial condition, Y is initial temperature matrix, changes as per problem%

% nx nodes in row, ny nodes in column%

lenx = 0.001; %length of the plate

leny = 0.001; %height of the plate

dx = 0.0001; %grid length in x

dy = 0.0001; % grid length in y

z = 0.001; %width of the plate

nx = lenx/dx; %total number of grids in x

ny = leny/dy; %total number of gris in y

t = 10; %total time

dt = 0.1; %length of time step

nt = t/dt; %number of time steps

kt = 15;%conductivity

Ae = dy*z;%area for calculating heat flux in east and west directions

Aw = Ae; %Area for east direction is equal to the area for west

As = dx*z; %area for calculating heat flux in north and south directions

An = As; %area for calcuating heat flux for south is equal to area for north

rho = 11; %density of the slab

C = 0.01; %specific heat of the slab

V = lenx*leny*z; %volume of the slab (metre^3)

h = 5000; %Convection heat transfer coefficient at the upper end(Water), W/(m2*K),

Tinf = 313.15; %Temperature of the fluid flowing over the plate

f = rho*V*C/dt; % a constant used in solving TDMAae = kt*Ae/dx; %Coefficient of Temperature for East Side

ae = kt*Ae/dx; % coefficient of temperature in the east side

aw = ae; %Coefficient of Temperature for East Side is equal to Temperature Coefficient in West side

as = kt*As/dy; %Coefficient of Temperature for South Side

an = as; %Coefficient of Temperature for South Side is equal to Temperature coefficient in North side

T = zeros(nx,ny); %Temperature matrix

T(1:nx,1)=363.15; %Boundary condition

T

Y = T; %Initial temperature matrix (will change after each iteration)

D = zeros(nx,ny); % utilized to calculate error later on

e=1; % defineded as 1 to start while loop

ha=1; % to get itretion number

A = zeros(1,nt+1);% central node temperature matrix over time

a=zeros(1,nx); % a,b,c,d utilizd to use TDMA

b=zeros(1,nx);

c=zeros(1,nx);

d=zeros(1,nx);

TN = zeros(1,nx);

TS = zeros(1,nx);

%% calculation started with 1st time step

for t=1:1:nt % 1st time step:1:last time step

t

while e > 10^-4 % 4 decimal accuracy to get converged answers.

%% TDMA in column

for j = 2:ny

if j<ny

for i=1:nx

if i>1 && i<nx

TS = T(:,j-1);

TN = T(:,j+1);

a(i)= aw + an + as + an + f;

b(i)= -ae;

c(i)= -aw;

d(i)= as*TS(i)+an*TN(i)+f*Y(i); %constant changes as per dx,dy,dt

else

a(nx)=f + ae + an + as;

a(1)= f + aw + an + as;

TS = T(:,j-1);

TN = T(:,j+1);

d(1) = f*Y(1) + as*T(1) + an*TN(1);

d(nx) = f*Y(nx) + an*TN(nx) + as*TS(nx);

end

else

for i = 1: nx

if i>1 && i<nx

TS = T(:,j-1);

TN = Tinf*ones(nx,1); %Tn for the upper boundary condition

a(i) = f + aw + ae + an + As*h;

b(i) = -ae;

c(i) = -aw;

d(i)= as*TS(i)+an*TN(i)+f*Y(i);

else

TS = T(:,j-1);

a(nx) = f + aw + as + As*h;

a(1) = f + ae + as + As*h;

d(1) = as*TS(1) + As*Tinf*h + f*Y(i);

d(nx) = as*TS(nx) + As*Tinf*h + f*Y(i);

end

X = TDMA(a,b,c,d);

T(:,j)=X;

end

for i=2:nx-1

for j=2:ny-1

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

- ((T(i,j)-Y(i,j))/dt);

end

e = max(max(abs(D)));

ha=ha+1;

end

T

A(t+1)= T((nx+1)/2,(nx+1)/2);

e=1;

Y=T;

end

%%graph for final time step,

u= 0:2/(nx-1):lenx;

v= 0:2/(ny-1):leny;

meshgrid(u,v);

contour(u,v,T)

title('temperature distribution')

figure()

t=0:nt

A(1)=0;

A

plot(t,A)

xlabel('time as 0.2 time interval')

ylabel('center tamperature')

figure()

s1=surf(u,v,T);

% creating function for TDMA %

function X = TDMA(a,b,c,d)

n = length(d);

% a is diagonal,b is super diagnonal,c is sub diagonal, d is the constant

% Elimination %

for i= 2:n

a(i)= a(i)-(c(i)*b(i-1))/a(i-1);

d(i)= d(i)-(d(i-1)*c(i))/a(i-1);

end

S(n)= d(n)/a(i);

%back substitute%

for i= n-1:-1:1

S(i)= (d(i)-b(i)*S(i+1))/a(i);

end

X = S;

end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

2 Dimensional heat equation using ADI method

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

类别

标签

Community Treasure Hunt

2 Dimensional heat equation using ADI method

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

回答（0 个）

类别

标签

另请参阅

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论