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

  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

  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

    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