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

Question

Sanley Guerrier 2024-1-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2067176-solving-heat-transfer-using-crank-nicolson-method-with-time-dependent-boundary-conditions

回答： SAI SRUJAN 2024-1-22

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 个评论
显示 1更早的评论隐藏 1更早的评论

Sanley Guerrier 2024-1-7

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 2024-1-7

If it was not accurate with "pdepe", you won't succeed with a different solution method, either. Sorry.

请先登录，再进行评论。

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

Answer 1

SAI SRUJAN 2024-1-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2067176-solving-heat-transfer-using-crank-nicolson-method-with-time-dependent-boundary-conditions#answer_1394221

在 MATLAB Online 中打开

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.

https://www.mathworks.com/help/matlab/ref/interp1.html

https://www.mathworks.com/help/matlab/ref/double.unique.html

I hope this helps!