clc
clear
close all
%COMPRESSION
% This script is for the calculation of effective compression modulus B (dimensionless)

% Step 1: defination of variables
h1_eff = 0:0.01:1;                                      % overall thickness of material 1
h1_eff = h1_eff';
[row, column] = size(h1_eff);
h2_eff = ones([row, column]) - h1_eff;                  % overall thickness of material 2

E_ratio = [0.01; 0.1; 1; 10; 100];                      % modulus ratio of materials 1/2
B_bar = zeros(row, length(E_ratio));

% Step 2: calculation of effective compression modulus B
for i = 1:5                                             % index for E_ratio

    B_bar(:,i) = 1./(h1_eff + E_ratio(i)*h2_eff);       % dimensionless effective compression modulus B

end

% Step 3: plot results
% Plot the influence of h1_eff on B_bar
figure(1)
for i = 1:5
    plot(h1_eff, B_bar(:,i), 'LineWidth',1);
    hold on
    
end
hold off

% Figure setup
xlim([0,1]);
ylim([0,10]);
xlabel('Effective h_{1} / h','FontSize',11);
ylabel('Normalized compression modulus E_{eff} / E_{1}','FontSize',11);
legend({'E_{ratio} = 0.01','E_{ratio} = 0.1','E_{ratio} = 1','E_{ratio} = 10','E_{ratio} = 100'},'FontSize',11);



%BENDING
clc
clear
close all
z_n = 1;                    %stack height
n = 4;                      %number of layers
z_vec = 0:(z_n/n):z_n;      
z_num = zeros(n);
z_den = zeros(n);
E_ratio = [0.01; 0.1; 1; 10; 100];

for i = [2:length(z_vec)]
    if mod(i,2) == 0            %even
        z_num(i) = E_ratio*((z_vec(i)^2/z_n^2)-(z_vec(i-1)^2/z_n^2));
        z_den(i) = E_ratio*((z_vec(i)/z_n)-(z_vec(i-1)/z_n));
    else                        %odd
        z_num(i) = ((z_vec(i)^2/z_n^2)-(z_vec(i-1)^2/z_n^2));
        z_den(i) = ((z_vec(i)/z_n)-(z_vec(i-1)/z_n));
    end
end
z_num
z_den
zc = sum(z_num)/(2*sum(z_den)) %location of neutral axis




E = zeros(n,1);
for i = [2:length(z_vec)]
    if mod(i,2) == 0                   % for even items
        E(i) = E_ratio*(((z_vec(i)/z_n)-(zc/z_n))^3-((z_vec(i-1)/z_n)-(zc/z_n))^3);
    else                                  % for odd items
        E(i) = (((z_vec(i)/z_n)-(zc/z_n))^3-((z_vec(i-1)/z_n)-(zc/z_n))^3);
    end
end
E

E_eff = 4*sum(E)