Solving heat transfer using Crank Nicolson method with time dependent boundary conditions
24 次查看(过去 30 天)
显示 更早的评论
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.
Appriciate your help.
%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
%Load the data.
load('mydata.mat', 'mydata');
%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
Hs = env(:,5);% net radiation
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 个评论
Torsten
2024-1-7
If it was not accurate with "pdepe", you won't succeed with a different solution method, either. Sorry.
回答(1 个)
SAI SRUJAN
2024-1-22
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!
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!