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

31 views (last 30 days)

Show older comments

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 Comments

Torsten
on 7 Jan 2024

### Answers (1)

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!

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!