Matlab scripts (scattering matrix for heat transfer)

%% This script calculates the average temperture rise of a heater on a film-substrate stack

%% using the scattering matrix formalism

%% developed by Tao Li, a PhD student in Prof. Zhen Chen's group at Southeast University of China.

%% posted Jun. 09, 2022

%% Related paper information:

%% Tao Li and Zhen Chen, Journal of Applied Physics 132,125103(2022)

%% https://doi.org/10.1063/5.0111267

clear all;

close all;

%% Thermal properties

kxy = [1,1]; %% kxy = kx/ky is the anisotropy ratio

ky = [1.5,150]; %% kx and ky are in-plane and the cross-plane thermal conductivities,[W/m-k]

C = [2,2]*1e6; %% C is the volumetric specific heat [J/m^3-K]

r_12 = 0; %% thermal contact resistance between layers 1 and 2,[m^2-K/W]

r = [0,r_12]; %% Assuming no thermal contact resistance between the heater and the first layer

%% Geometric parameters

d = [3e-6,500e-6]; %% thickness of layer, [m]

L = 2e-3; %% length of the heater, [m]

b = 50e-6; %% half width of the heater, [m]

%% the number of layers, boundary condition

n = 2; %% Number of layers

bound = 0; %% semi-infinite:0; adiabatic:1; isothermal:-1 Eq.23

%% Electrical parameters

P = 1; %% Root-mean-square value of the Joule heating power, [W]

w = logspace(0,12,50); %% Current angular frequency [rad/s]

%% Calculate the average temperature rise of the heater by different methods

%------------------- Scattering Matrix,J. Appl. Phys. 132, 125103 (2022); -----------------

for i = 1:1:length(w)

T_int_SM(i) = quadv(@(lambda)SM_integ(lambda,w(i),b,kxy,ky,C,d,n,bound,r),0,1e8,1e-10); %% lambda is integral variable

T_SM(i) = P./pi./L.*real(T_int_SM(i));

end

%------------------- Transfer matrix, J. Appl. Phys. 86, 3959-3963 (1999).---------------

for i = 1:1:length(w)

T_int_TM(i) = quadv(@(lambda)HTM_integ(lambda,w(i),b,ky,kxy,C,d,n,bound),0,1e8,1e-10);

T_TM(i) = P./pi./L.*real(T_int_TM(i));

end

%% plot

figure

semilogx(w,T_TM,'bo',w,T_SM,'r','linewidth',2.5);

h=legend('Transfer Matrix','Scattering Matrix');

set(h,'FontSize',16);

xlabel('Angular frequency of electrical current [rad/s]','FontSize',16)

ylabel('Temperature rise of heater [K]','FontSize',16)

%% SM_integ

function F = SM_integ(lambda,w,b,kxy,ky,C,d,n,bound,r)

Lambda = lambda; %% Performing the Fourier cosine transform, x-->lambda

L_Lambda = length(Lambda);

%-------------------- layers above the heater -------------------------

% Here we assume vacuum or air above the heater, as shown in the second row of Table 2

for m = 1:1:L_Lambda

M11_a = 1;

M12_a = 1;

M21_a = 0;

M22_a = 0;

M_a{m} = [M11_a,M12_a;

M21_a,M22_a]; % M-matrix of above layers, Eq.13

end

for m = 1:1:L_Lambda

S_a{m} = [1,0;0,1]; %% scattering matrix of above layers

end

%---------------------------layers below the heater ---------------------------------

for j = 1:1:n

alpha_b(j) = ky(j)/C(j); %% thermal diffusivity,

q_b(j) = sqrt(1i*2.*w./alpha_b(j)); %% Coefficient 2 comes from the conversion between current angular frequency and heat flux angular frequency

for m = 1:1:L_Lambda

u_b(j) = sqrt(kxy(j).*Lambda(m).^2+q_b(j).^2); %% Eq. 31

gamma_b(j) = ky(j).*u_b(j); %% Eq. 14

L_b(j) = u_b(j).*d(j);

f_b{j,m} = exp(-L_b(j));

M11_b(j) = 1;

M12_b(j) = 1;

M21_b(j) = gamma_b(j);

M22_b(j) = -gamma_b(j);

M_b{j,m} = [M11_b(j),M12_b(j);

M21_b(j),M22_b(j)]; % M-matrix of below layers, Eq.13

M11_b_r(j) = 1-r(j).*gamma_b(j);

M12_b_r(j) = 1+r(j).*gamma_b(j);

M21_b_r(j) = gamma_b(j);

M22_b_r(j) = -gamma_b(j);

M_b_r{j,m} = [M11_b_r(j),M12_b_r(j);

M21_b_r(j),M22_b_r(j)]; %% Intoroduce thermal contact resistance, Eq.18

end

for j = 1:1:n-1

for m = 1:1:L_Lambda

I_b{j,m} = M_b{j,m}\M_b_r{j+1,m}; % interface matrix of above layers, Eq.17

end

for m = 1:1:L_Lambda

S_b{1,m} = [1,0;0,1]; %% S-matrix of the first layer below the heater

end

for j = 2:1:n

for m = 1:1:L_Lambda

I_0_b = I_b{j-1,m};

I11_b = I_0_b(1,1);

I12_b = I_0_b(1,2);

I21_b = I_0_b(2,1);

I22_b = I_0_b(2,2);

S_0_b = S_b{j-1,m};

S11_b = S_0_b(1,1);

S12_b = S_0_b(1,2);

S21_b = S_0_b(2,1);

S22_b = S_0_b(2,2);

f_b_j_1 = f_b{j-1,m};

f_b_j = f_b{j,m};

S11_b_j = S11_b.*f_b_j_1./(I11_b-S12_b.*I21_b.*f_b_j_1);

S12_b_j = (S12_b.*I22_b.*f_b_j_1.*f_b_j-I12_b.*f_b_j)./(I11_b-S12_b.*I21_b.*f_b_j_1);

S21_b_j = S21_b + S22_b.*I21_b.*S11_b_j;

S22_b_j = S22_b.*(S12_b_j.*I21_b+I22_b.*f_b_j); %% Table 1

S_b{j,m} = [S11_b_j,S12_b_j;

S21_b_j,S22_b_j]; %% S-matrix below the heater

end

%----------------- calculate the temperature of heater ------------------------

for m = 1:1:L_Lambda

M_a_na = M_a{m};

M11_i = M_a_na(1,1);

M12_i = M_a_na(1,2);

M21_i = M_a_na(2,1);

M22_i = M_a_na(2,2); %% M-matirx of the nearest layer above the heater

M_b_1 = M_b{1,m};

M11_i_1 = M_b_1(1,1);

M12_i_1 = M_b_1(1,2);

M21_i_1 = M_b_1(2,1);

M22_i_1 = M_b_1(2,2); %% M-matirx of the nearest layer below the heater

S_i_a = S_a{m};

S11_i_a = S_i_a(1,1);

S12_i_a = S_i_a(1,2);

S21_i_a = S_i_a(2,1);

S22_i_a = S_i_a(2,2); % Above the heater:S-matrix from the top layer to the layer closest to the heater as shown in fig.2(1-->j)

S_N_b = S_b{n,m};

S11_N_b = S_N_b(1,1);

S12_N_b = S_N_b(1,2);

S21_N_b = S_N_b(2,1);

S22_N_b = S_N_b(2,2); % below the heater:S-matrix from the layer closest to the heater to the The bottom layer as shown in fig.2(j+1-->N)

alpha = 0; %% Eq.28. here we assume a vacuum or air above the heater

f_b_n = f_b{n,m};

beta = S21_N_b+bound*(S11_N_b*S22_N_b*f_b_n)./(1-bound*S12_N_b*f_b_n); %%% Eq.29

f_b_1 = f_b{1,m};

f_a_na = 0;

H11 = M11_i_1+M12_i_1.*beta.*f_b_1;

H12 = -(M11_i.*alpha.*f_a_na+M12_i);

H21 = M21_i_1+M22_i_1.*beta.*f_b_1;

H22 = -(M21_i.*alpha.*f_a_na+M22_i);

H_0 = H11.*H22-H12.*H21; %% Table 2

a_i_1 = (M11_i.*alpha.*f_a_na+M12_i)./H_0;

b_i_1 = (M11_i.*alpha.*f_a_na+M12_i)./H_0.*beta; %% The third and fourth rows in Table 2

T_1 = a_i_1+b_i_1.*f_b_1; %% The heater temperature is calculated by the upper interface of layer j+1

% It's equivalent to using the j-layer lower interface to calculate the heater temperature(the first and second rows in Table 2)

T_integ_SM_1(m) = T_1.*(sin(b.*Lambda(m))).^2./(b.*Lambda(m)).^2; %% 3Omega method, performing the Fourier cosine inverse transform(lambda-->x),

% Averaging the temperature across the heater width

end

F = T_integ_SM_1;

end

%% HTM_integ J. Appl. Phys. 86, 3959-3963(1999).

function F = HTM_integ(lambda,w,b,k,kxy,C,d,n,bound)

gamma_a = 0;

for j = 1:1:n

alpha(j) = k(j)/C(j);

q(j,:) = sqrt(1i*2.*w./alpha(j));

u(j,:) = sqrt(kxy(j).*lambda.^2+q(j,:).^2);

gamma_b(j,:) = k(j).*u(j,:);

L(j,:) = u(j,:).*d(j);

end

N = length(u(1,:));

U_G = [ones(1,N); zeros(1,N);zeros(1,N); ones(1,N)];

for j = n:-1:2

u11 = exp(L(j-1,:));

u12 = zeros(1,N);

u21 = zeros(1,N);

u22 = exp(-L(j-1,:));

U = [u11;u12;u21;u22];

G11 = 1/2.*(gamma_b(j-1,:)+gamma_b(j,:))./gamma_b(j-1,:);

G12 = 1/2.*(gamma_b(j-1,:)-gamma_b(j,:))./gamma_b(j-1,:);

G21 = 1/2.*(gamma_b(j-1,:)-gamma_b(j,:))./gamma_b(j-1,:);

G22 = 1/2.*(gamma_b(j-1,:)+gamma_b(j,:))./gamma_b(j-1,:);

Gamma = [G11; G12; G21;G22];

U_G0 = [u11.*G11; u11.*G12; u22.*G21; u22.*G22];

U_G = [U_G0(1,:).*U_G(1,:)+U_G0(2,:).*U_G(3,:);

U_G0(1,:).*U_G(2,:)+U_G0(2,:).*U_G(4,:);

U_G0(3,:).*U_G(1,:)+U_G0(4,:).*U_G(3,:);

U_G0(3,:).*U_G(2,:)+U_G0(4,:).*U_G(4,:)];

end

T11_j_N = U_G(1,:);

T12_j_N = U_G(2,:);

T21_j_N = U_G(3,:);

T22_j_N = U_G(4,:);

G11_a_b = 1/2.*(gamma_b(1,:)+gamma_a)./gamma_b(1,:);

G12_a_b = 1/2.*(gamma_b(1,:)-gamma_a)./gamma_b(1,:);

G21_a_b = 1/2.*(gamma_b(1,:)-gamma_a)./gamma_b(1,:);

G22_a_b = 1/2.*(gamma_b(1,:)+gamma_a)./gamma_b(1,:);

Gamma_a_b = [G11_a_b; G12_a_b; G21_a_b;G22_a_b];

alpha = 0;

beta = (T11_j_N+T12_j_N.*bound.*exp(-2.*L(n,:)))./(T21_j_N+T22_j_N.*bound.*exp(-2.*L(n,:)));

b_j_1 = 1./2./gamma_b(1,:).*(-(alpha.*G21_a_b+G22_a_b)-(alpha.*G11_a_b+G12_a_b))./((alpha.*G11_a_b+G12_a_b)-beta.*(alpha.*G21_a_b+G22_a_b));

a_j_1 = beta.*b_j_1;

F1 = a_j_1+ b_j_1;

F = F1.*(sin(b.*lambda)).^2./b.^2./lambda.^2;

end