Solving heat transfer using Crank Nicolson method with time dependent boundary conditions

27 views (last 30 days)
Sanley Guerrier on 7 Jan 2024
Answered: SAI SRUJAN on 22 Jan 2024
Dear experts-
I am trying to solve the heat transfer model using Crank Nicolson method, the inputs for the heat fluxes Hs, Ta, hc, and Ts are column vectors (time dependent). I don't know how to deal with that. Can someone help, please?
I updated the code hopefully it is more explicit.
I got an error with this command:
T1(1:n, 1:ny+1,1)= Tin;
Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.
I am not confident about the IC, BCs, and Thomas function.
%Diffusion equation - 1D- Crank Nicolson
% Neumann BC conditons - Time dependent heat flux at boundaries
% Unconditionally stable
% A*T_n+1 = B*T-n + 4*d*delta_x*T_BC = C
% IC : T(x, t=0) = Tin
%BCs: at x=0; fa = qin; at x=L; fb = 0
clear
clc
clf
%data = mydata; % use all of the data.
data = mydata(1:1000, :); % extract a small subset for testing.
temp=data{:, 5:8};
Ts = temp(:,1); % surface temperature
env = data{:,11:end};
w=env(:,9);% wind speed
Ta=env(:,1); % air temperature
hc=7.4+6.39*w.^0.75; %convective heat transfer coefficient (W/m2 K)
kt = 1.2; % thermal conductivity [W/(m*K)]
rho = 2200; % density [kg/m^3]
Cp = 920; % specific heat capacity [J/(kg*K)]
alpha = kt/(rho * Cp)*86400;
time = datetime(data.DAY) + (data.HOUR + data.QHR/4)/24;
nx=106;% Number of gridpoints
xlength = 8.52 * 0.0254; % convert distance from inch to meter
delta_x = xlength/nx;
ylength = 8.52 * 0.0254;
ny= 106;
delta_y = ylength/ny;
t=104.1563;%total time day
nt=1000;
delta_t = t/nt; %Timestep
t_mesh = convertTo(time, 'datenum') - convertTo(time(1), 'datenum');
%compute the heat flux fa
qconv = hc.*(Ts-Ta);
fa=qconv + Hs; % heat flux at x =0
fb=0; % heat flux at x = L
plot(time,qconv)
plot(time,Hs)
z_d = [ 0.96, 2.88, 6, 8.52] * 0.0254;
x = unique([linspace(z_d(1), z_d(end), 106), z_d]);
%initial temperature
Tin= interp1(z_sensor, temp(1, :), x, 'linear')';
%solution
n=nx+1; % total no. of nodes/ point in the domain
d = alpha*delta_t/delta_x^2; %diffusion number
x(1) = 0;
for i=1:n
x(i)=x(1) + (i-1)*delta_x;
end
x;
%creating initial and boundary conditions B vector
for k=1:1
for i=1:n
IC(i,k)=Tin(i);
end
end
IC;
for k=1:nt
for i=1:n
if i==1
BC(i,k)= -fa(k)/kt;
else
BC(i,k)=0;
end
end
end
BC
%creating B matrix
B= zeros(n,n);
for i =1:n
for j=1:n
if i==j
B(i,j)=2*(1-d);
elseif(i==j+1)&(i<n)
B(i,j)=d;
elseif(i==j+1)&(i==n)
B(i,j)=2*d;
elseif(i==j-1)&(j==2)
B(i,j)=2*d;
elseif(i==j-1)&(j>2)
B(i,j)=d;
else
B(i,j)=0;
end
end
end
B;
C=B*IC + 4*d*delta_x*BC
%creating A matrix
% A matrix diagonal, upper and lower values
for i=1:n
dA(i)=2*(1+d);
if (i==1)
uA(i)= -2*d;
elseif (i<n) && (i>1)
uA(i)= -d;
else
uA(i)=0;
end
if (i>1) && (i<n)
pA(i)= -d;
elseif (i==n)
pA(i)=-2*d;
else
pA(i)=0;
end
end
uA;
dA;
pA;
%creating T vector
% T1 = A\C
T1(1:n, 1:ny+1,1)= Tin;
for k = 2:nt+1
T = ThomasAlgorithm(dA, pA, uA, C);
T = T';
for j = 1:ny+1
T1(:,j,k)=T;
end
C = B*T + 2*d*delta_x*BC;
end
T1;
n;
Tmin = min(T1(:));
Tmax = max(T1(:));
T2 = T1
% Temperature profile
subplot(1,2,1)
%initial Temperature profile
T3 = T2(:,:,1);
subplot(2,2,1)
T3 =rot90(T3);
x=0:delta_x:xlength;
y = 0:delta_y:ylength;
imagesc(x,y,T3)
colorbar;
title(['Temperature Profile -', 'T in deg.C','@ time(t) =', num2str(0), 'day'])
xlabel(['x in m'])
set(gca, 'ytick', [])
%final Temperature profile
T4 = T2(:,:, nt+1);
subplot(2,2,3)
T4 = rot90(T4);
x = 0:delta_x:xlength;
y=0:delta_y:ylength;
imagesc(x,y,T4)
colorbar;
title(['Temperature Profile -','T in deg.C','@ time(t) =', num2str(t), 'day'])
xlabel(['x in m'])
set(gca, 'ytick',[])
3 CommentsShow 1 older commentHide 1 older comment
Sanley Guerrier on 7 Jan 2024
I already used pdepe, the solution I got is not too accurate, so I would like to compare the solutions.
If you could help it would be appriciated.
Torsten on 7 Jan 2024
If it was not accurate with "pdepe", you won't succeed with a different solution method, either. Sorry.

SAI SRUJAN on 22 Jan 2024
Hi Sanley,
I understand that you are facing an error in solving the one-dimensional heat transfer problem using the Crank-Nicolson method with time-dependent boundary conditions.The error message indicates that the sizes of the matrices on the left and right sides of the assignment do not match.
The issue stems from a mismatch in matrix dimensions in the assignment statement, The dimensions of `Tin` must match the dimensions of the submatrix `T1(1:n, 1:ny+1,1)`.
nx = 106;
n = nx + 1;
ny = 106;
z_d = [ 0.96, 2.88, 6, 8.52] * 0.0254;
x = unique([linspace(z_d(1), z_d(end), 106), z_d]);
Tin = interp1(z_sensor, temp(1, :), x, 'linear')';
T1(1:n, 1:ny+1,1) = Tin;
The varibale 'z_d' is a 1 x 4 array, while the variable 'x' is generated using 'linspace' and 'unique', resulting in a column vector whose length may vary but will typically be around 106 to 110. 'Tin' is generated using 'interp1' on 'x', therefore its length will be the same as the length of 'x'.
In the context of the given assignment, there is an expectation to populate the submatrix 'T1(1:n, 1:ny+1,1)' with a 107x107 matrix. However, 'Tin' is a column vector with a length corresponding to 'x', which presents a dimension mismatch.
To resolve this, you should ensure that `Tin` is appropriately reshaped or replicated to match the dimensions of the `T1` submatrix.
For a comprehensive understanding of the 'unique' and 'interp1' function in MATLAB, please refer to the following documentation.
I hope this helps!

Categories

Find more on Geometry and Mesh in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!