%% TOMAS MIVALT
%% VYPOCET DELENE VETKNUTE DESKY - ANALYTICKY VYPOCET

clc
clear all;


%% Vstupni udaje
D = 81; % prumer vnejsi desky [mm]   
Dr = 56.4; % prumer vnitrni desky [mm]
h2 = 9.8; % tloustka dna pistu - vnejsi deska
h1 = 7.1; % tloustka dna pistu - vnitrni deska


p = 6; % MPa - zatizeni

% Materialove charakteristiky
% Slitina hliniku - AlSi12CuNiMg   
E = 0.73*(10^5); % [MPa] Modul pruznosti    
u = 0.33; % Poissonova konstanta 

%% Vypocet

% Pouzivany max polomer 
Rmax = D/2;
% Pouzivany delici polomer 
Rdel = Dr/2;

% Ohybova tuhost desky 2 - vetsi
B2 = (E*h2^3)/(12*(1-(u^2)));
% Ohybova tuhost desky 1 - mensi
B1 = (E*h1^3)/(12*(1-(u^2)));


% !!! Moment, ktery vznika na delicim polomeru vlivem pusobeni tlaku na
% svislou plochu osazeni
M1 = -p*(h2-h1)*(h1/2);

% konstany - pocitano v maple
% DESKA 1
C1 = -(1/16)*(-B1*Rmax^2*u*p*Rdel^2+16*B1*Rdel^2*M1+p*Rdel^4*B1-u*p*Rdel^4*B2-p*Rdel^2*B2*Rmax^2-p*Rdel^4*B2+u*p*Rdel^4*B1-2*B1*p*Rmax^4-16*B1*Rmax^2*M1+B1*Rmax^2*p*Rdel^2+p*Rdel^2*B2*u*Rmax^2)/((-Rdel^2*B2-u*Rdel^2*B2+u*Rmax^2*B2-Rmax^2*B2+Rdel^2*B1-Rmax^2*B1+Rdel^2*u*B1-u*B1*Rmax^2)*B1)
C2 = 0

% DESKA 2
C01 = -(1/16)*(B2*u*p*Rmax^4-u*p*Rdel^4*B2+16*B2*M1*Rdel^2-u*B1*p*Rmax^4-p*Rdel^4*B2-B2*p*Rmax^4-B1*p*Rmax^4+p*Rdel^4*B1+u*p*Rdel^4*B1)/(B2*(-Rdel^2*B2-u*Rdel^2*B2+u*Rmax^2*B2-Rmax^2*B2+Rdel^2*B1-Rmax^2*B1+Rdel^2*u*B1-u*B1*Rmax^2))
C02 = (1/16)*Rmax^2*Rdel^2*(-u*p*Rdel^2*B2+16*B2*M1-p*Rdel^2*B2+p*Rdel^2*B1+u*p*Rdel^2*B1+p*B2*Rmax^2+p*B2*u*Rmax^2-B1*p*Rmax^2-B1*Rmax^2*u*p)/(B2*(-Rdel^2*B2-u*Rdel^2*B2+u*Rmax^2*B2-Rmax^2*B2+Rdel^2*B1-Rmax^2*B1+Rdel^2*u*B1-u*B1*Rmax^2)) 


% Natoceni 
v1 = @(r)(C1*r+(p*(r^3))/(16*B1));
v2 = @(r)(C01*r+C02/r+(p*(r^3))/(16*B2));

% Der. natoceni
dv1 = @(r)(C1+(3*p*(r^2))/(16*B1));
dv2 = @(r)(C01-C02/(r^2)+(3*p*(r^2))/(16*B2));

% Momenty 
Mr1 = @(r)(-B1*(dv1(r)+u*(v1(r)/r)));
Mt1 = @(r)(-B1*((v1(r)/r)+u*dv1(r)));

Mr2 = @(r)(-B2*(dv2(r)+u*(v2(r)/r)));
Mt2 = @(r)(-B2*((v2(r)/r)+u*dv2(r)));

% Napeti - maximalni v danem miste .. z = -h/2  - horni strana
sr1 = @(r)((12*Mr1(r))/(h1^3))*(-h1/2);
st1 = @(r)((12*Mt1(r))/(h1^3))*(-h1/2);
sr2 = @(r)((12*Mr2(r))/(h2^3))*(-h2/2);
st2 = @(r)((12*Mt2(r))/(h2^3))*(-h2/2);
sz = -p;

%% Vypis dat
% rozliseni
n = 1000;

int2 = linspace(-Rmax,-Rdel,n);
int1 = linspace(-Rdel,0,n);
int01 = linspace(0,Rdel,n);
int02 = linspace(Rdel,Rmax,n);

for i=1:n
  dat_v1(i) = v1(int1(i));
  dat_mr1(i) = Mr1(int1(i));
  dat_mt1(i) = Mt1(int1(i));
  
  dat_sr1(i) = sr1(int1(i));
  dat_st1(i) = st1(int1(i));
  dat_sz1(i) = sz;
  
  A = [dat_sr1(i) dat_st1(i) dat_sz1(i)];
  B = sort(A);
  sigma1 = B(3);
  sigma2 = B(2);
  sigma3 = B(1);
    
  dat_sred1(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));

end;

for i=1:n
  dat_v2(i) = v2(int2(i));
  dat_mr2(i) = Mr2(int2(i));
  dat_mt2(i) = Mt2(int2(i));
  
  dat_sr2(i) = sr2(int2(i));
  dat_st2(i) = st2(int2(i));
  dat_sz2(i) = sz;
  
  A = [dat_sr2(i) dat_st2(i) dat_sz2(i)];
  B = sort(A);
  sigma1 = B(3);
  sigma2 = B(2);
  sigma3 = B(1);
    
  dat_sred2(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));
end;

for i=1:n
  dat_v01(i) = v1(int01(i));
  dat_mr01(i) = Mr1(int01(i));
  dat_mt01(i) = Mt1(int01(i));
  
  dat_sr01(i) = sr1(int01(i));
  dat_st01(i) = st1(int01(i));
  dat_sz01(i) = sz;
  
  A = [dat_sr01(i) dat_st01(i) dat_sz01(i)];
  B = sort(A);
  sigma1 = B(3);
  sigma2 = B(2);
  sigma3 = B(1);
    
  dat_sred01(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));

end;

for i=1:n
  dat_v02(i) = v2(int02(i));
  dat_mr02(i) = Mr2(int02(i));
  dat_mt02(i) = Mt2(int02(i));
  
  dat_sr02(i) = sr2(int02(i));
  dat_st02(i) = st2(int02(i));
  dat_sz02(i) = sz;
  
  A = [dat_sr02(i) dat_st02(i) dat_sz02(i)];
  B = sort(A);
  sigma1 = B(3);
  sigma2 = B(2);
  sigma3 = B(1);
    
  dat_sred02(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));
end;


%% Slozeni finalnich vektoru pro vypis grafu

  dat_v = [dat_v2 dat_v1 dat_v01 dat_v02];
  dat_mr = [dat_mr2 dat_mr1 dat_mr01 dat_mr02];
  dat_mt = [dat_mt2 dat_mt1 dat_mt01 dat_mt02];
  
  dat_sr = [dat_sr2 dat_sr1 dat_sr01 dat_sr02];
  dat_st = [dat_st2 dat_st1 dat_st01 dat_st02];
  dat_sz = [dat_sz2 dat_sz1 dat_sz01 dat_sz02];
  dat_sred = [dat_sred2 dat_sred1 dat_sred01 dat_sred02];
  
  int = [int2 int1 int01 int02];


%% Kresleni grafu

figure(1);
plot(int,dat_v);
xlabel('r [mm]');
ylabel('Natoceni [Rad]');
grid on;

figure(2);
plot(int,dat_mr,int,dat_mt);
xlabel('r [mm]');
ylabel('Liniovy moment [N*mm/mm]');
legend('Radialni', 'Tecny');
grid on;

figure(3);
plot(int,dat_sr,'--',int,dat_st,'--', int,dat_sz,'--',int,dat_sred);
xlabel('r [mm]');
ylabel('Napětí \sigma [MPa]');
legend('\sigma_r', '\sigma_t', '\sigma_z', '\sigma_{red}');
grid on;


