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

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',[])

[Torsten]

If you have no experience with numerical methods and just take any code from the internet to adopt it to your problem, you should use the more reliable way to use MATLAB's "pdepe". It's easy to use and trustworthy.

[Sanley Guerrier]

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]

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

Answer 1

[SAI SRUJAN]

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!