% % Name: Meera Garud
% % NE211: Project
% % Title Coding of effects of Thermoelastic damping
% % Date: 29th Nov '14
% % All mathematical relation are refered from : Thermoelastic damping in micro- and nanomechanical systems
% % - Ron Liftshitz and M.L. Roukes
% % - Physical review B, volume 61. number 8, 2000


clc;
clear all;
close all;

L = 100e-6;     % Length of the beam in meter
b = 20e-6;       % Width of the beam in meter
c = 2e-6;       % Thickness of the beam in meter
A = b*c;        % Cross sectional area of the beam in meter^2
I = b*(c^3)/12; % Mechanical contribution to momet of inertia in meter^4
E = 160e9;      % Young's modulus of elasticity for single crystal silicon in Pa
m = 50e-9;      % Mass of beam in Kg
rho = 2300;     % Mass density of beam in Kg/cubic-meter
alpha = 2.6e-6; % Thermal coefficient of expansion per kelvin
T = 300;        % Ambient temperature in kelvin

Elements =15;
n =  Elements + 1;
i = 1:1:n;
nodes = linspace(0,L,n);
N = [i' nodes'];

connectivity_array = zeros(Elements,3);
for j = 1:Elements
    connectivity_array(j,1) = i(j);
    connectivity_array(j,2) = i(j);
    connectivity_array(j,3) = i(j+1);
end
l = nodes(2)-nodes(1);

% Analytical solution for mode shapes

Beta_n = [4.73004074 7.85320462 10.9956078 14.1371655];
Sn = [0.982502214576238 1.000777311907269 0.999966450125409 1.000001449897657];

beta = Beta_n/L;

dx = 0.01e-6;
x = 0:dx:L;
x_normalised = x/max(x);

for j= 1:4
    MS(:,j) = cosh(beta(j).*x) + cos(beta(j).*x) - (Sn(j))*(sinh(beta(j).*x) + sin(beta(j).*x));
    
    MS_max(j) = max(max(MS(:,j),abs(min(MS(:,j)))));
    MS_Normalised(:,j) = MS(:,j)./MS_max(j);
end

% Calculating isothermal eigen values

for j = 1:4;
    w(j) = (beta(j)^2)*sqrt(E*I/(rho*A));
    f(j) = w(j)./(2*pi);
end
Natural_frequecies_for_first_four_mode_shapes = f'
 
% Plotting of mode shapes, first four mode shapes, for beam with both ends
% fixed

for j = 1:4
    figure(j)
    plot(x_normalised,MS_Normalised(:,j),'Linewidth',2);
    grid on;
    xlabel('Length (x)');
    ylabel('Mode shape function (U(x))');
    title(sprintf('Natural frequency=%0.3fHz',f(j)));
end

% Data to find out frequency dependent Young's modulus

i = sqrt(-1);
Xd = 0.86e-4;   % Thermal diffusivity at 300K for single crystal silicon

% Calculation of frequency dependent Young's modulus

for j = 1:4
    k(j) = (1+i)*sqrt(w(j)/(2*Xd));
    f_w(j) = (24/(b*c)^3)*((b*k(j)/2) - tan(b*k(j)/2));
    C = 713;    % Specific heat in J/kg-K for SCS
    delta_E = E*(alpha^2)*T/C;
    Ew(j) = E*(1+delta_E*(1+f_w(j)));
end

% Finding new eigen frequencies
 
for j = 1:4;
    w_TED(j) = ((beta(j))^2)*sqrt(Ew(j).*I/(rho*A));
    f_TED(j) = w_TED(j)./(2*pi);
end

Natural_frequecies_for_first_four_mode_shapes_with_TED = real(f_TED)'

% Calculating inverse of quality factor i.e finding thermoelastic damping

for j = 1:4
        Q_inv_withoutTED(j) = 2* abs(imag(w(j))/real(w(j)));
        Q_inv(j) = 2*abs(imag(w_TED(j))/real(w_TED(j)));
end

% Set of different widths to observe the variation in TED effects
b = [10e-6 20e-6 40e-6 80e-6];

% Calculation of inverse of quality factor using value of zeta 
for j = 1:4
    zeta(j) = b(j).*sqrt(w(1)/(2*Xd));
    Re_w(j) = w(1)*(1+(delta_E/2)*(1-(6/(zeta(1))^2)*((sinh(zeta(1))-sin(zeta(1)))/(cosh(zeta(1))+cos(zeta(1))))));
    Im_w(j) = w(1)*(delta_E/2)*((6/(zeta(1))^2)*((sinh(zeta(1))+sin(zeta(1)))/(cosh(zeta(1))+cos(zeta(1))))-(6/(zeta(1)^2)));
    Q_inv_z(j) = (E*(alpha^2)*T/C)*((6/((zeta(1))^2))-(6/(zeta(1))^3)*((sinh(zeta(1))+sin(zeta(1)))/(cosh(zeta(1))+cos(zeta(1)))));
end

% Normalising
    q_inv_max_by_d_E = max(max(real(Q_inv(j)/delta_E)),abs(min(real(Q_inv(j)/delta_E))));

    normalised_TED = real(Q_inv/delta_E)./q_inv_max_by_d_E;
    re_zeta = real(zeta);
    log_zeta = log(real(zeta));
    log_normalised_TED = log(normalised_TED);
    TED = real(Q_inv);
  
    Q_inv = 2*abs(imag(w_TED)./real(w_TED));


    figure(4)
    plot(re_zeta,TED)
    xlabel('zeta');
    ylabel('Qinv by delta E');
  
    
    figure(5)
    plot(log_zeta,log_normalised_TED)
    xlabel('zeta');
    ylabel('Qinv by delta E');

    mode = 1:1:4;
    re_w_TED = real(w_TED);
    w_TED_max = max(re_w_TED);
    w_TED_nor = re_w_TED./w_TED_max;
    w_max = max(w);
    w_nor = w./w_max;
    
% %     f = w/(2*pi);
% %     re_f_TED = re_w_TED/(2*pi);
     

% Observing effect of TED on natural frequency
figure(6)
stem(f)
hold on
stem(re_f_TED,'fill')
legend('withot TED ','With TED')
xlabel('mode');
ylabel('Natural Frequency in Hz');
title('Comparison of natural frequencies for first 4 mode shapes of euler beam\n with fixed-fixed boundary condition')

f = w/(2*pi);
f_TED = w_TED /(2*pi);
    
    
% Plot to show that TED is observed only in materials having finite thermal
% conductivity, otherwise system shows isothermal behaviour
figure(7)
plot(mode,f,mode,f_TED)
legend('withot TED ','With TED')
xlabel('mode');
ylabel('Natural Frequency Hz');
title('Effect of very high thermal conductivity on TED')