%% TOMAS MIVALT
%% MODEL PODEPRENA DESKA - ANALYTICKY VYPOCET

clc
clear all

%% VSTUPNI HODNOTY
D = 81; % prumer pistu mm   
h = 8.45; % tloustka dna pistu
p = 6; % [MPa] pracovni tlak pri expanzi = 5-8 MPa % = 6

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

%% Vypocet

% Pouzivany max polomer 
R = D/2;

% Ohybova tuhost skorepiny
B = (E*h^3)/(12*(1-(u^2)));

% Posouvajici sila T(r) v obecnem rezu
% ze silove podminky Sum(Fz)=0 .... 2*pi*r*T(r) + pi*p*r^2 = 0
T = @(r)(-p*r/2);

%% Partikularni reseni pruhybove cary je Vpart = -1/r*B * ( int( r*int(T(r))dr )dr )
% Po dvou integracich Vpart = -1/rB * (-p*r^4/16)
% Po uprave Vpart = p*r^3/16B
Vpart = @(r)((p*(r^3))/(16*B));

% Vzhledem k symetrickemmu zatizeni na r=0 je v=0 => C2 = 0

% Posuv je pouze integral z natoceni 
% W = int V dr + C3

% Zbyva zjistit konstanty C1 a C3
% Deska je podeprena -> Okrajove podminky:
% r = R ... W=0 .. na podpore je posunuti 0
% r = R ... Mr=0 .. na podpore musi byt radialni moment rovny 0

% W = (C1*r^2/2 + p*r^4/64B) + C3
% V = C1*r+p*r^3/16B

% Z vztahu pro radialni moment a okrajove podminky lze vyjadrit tento vztah
% a tim urcit rovnici natoceni
C1 = -((p*(R^2)*(3+u))/(16*B*(1+u)));
% Natoceni
V = @(r)(C1*r+Vpart(r));

% Druhou konstantu C3 urcime z druhe okrajove podminky - r = R ... W=0 .. na kraji je posunuti 0
C3 = (p*(R^4)*(u+5))/(64*B*(u+1));

% Posuv
W = @(r)((C1)*((r^2)/2)+(p*(r^4))/(64*B)+C3);

%% Radialni a obvodovy moment

% Nachystani derivace r
der_V = @(r)(C1+(3*p*(r^2))/(16*B));

% Radialni M
Mr = @(r)(-B*(der_V(r)+(u*V(r)/r)));
% Obvodovy M
Mt = @(r)(-B*((V(r)/r)+u*der_V(r)));


%% Napeti - maximalni v danem miste .. z = -h/2 (horni strana pistu)

sr = @(r)((12*Mr(r))/(h^3))*(-h/2);
st = @(r)((12*Mt(r))/(h^3))*(-h/2);
% horni strana desky - osove napeti rovno tlaku p
sz = -p;

%% Vypis hodnot a grafy

% rozliseni
n = 1000;
misto_rezu = linspace(-R,R,n);


for i=1:n

% Deformace
data_T(i) = T(misto_rezu(i));
data_W(i) = W(misto_rezu(i));
data_V(i) = V(misto_rezu(i));

% Momenty
data_mr(i) = Mr(misto_rezu(i));
data_mt(i) = Mt(misto_rezu(i));

% napeti
data_sr(i) = sr(misto_rezu(i));
data_st(i) = st(misto_rezu(i));
data_sz(i) = sz;

end;


% SIGMA redukovane dle energeticke hypotezy HMH
t1 = 0;
t2 = 0;
t3 = 0;

for i=1:n
 
    t1 = sr(misto_rezu(i));
    t2 = st(misto_rezu(i));
    t3 = sz;
    
    A = [t1 t2 t3];
    B = sort(A);
    
    sigma1 = B(3);
    sigma2 = B(2);
    sigma3 = B(1);
    
    sigma_red(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));
    sigma_red_maxt(i) = sigma1-sigma3;
    
end;

%% Kresleni grafu 

figure(1);
plot(misto_rezu, data_W);
title('Posuv [mm]');
xlabel('r');
ylabel('W(r)');
grid on;

figure(2);
plot(misto_rezu, data_V);
title('Natoceni [rad]');
xlabel('r');
ylabel('V(r)');
grid on;


figure(3);
plot(misto_rezu, data_sr, misto_rezu, data_st, misto_rezu, data_sz,misto_rezu, sigma_red);
title('Prubeh napeti [MPa]');
xlabel('poloměr r');
ylabel('Napětí \sigma(r)');
legend('\sigma_r(r) radiální','\sigma_t(r) obvodové','\sigma_z(r) osové','\sigma_{RED}(r) redukované');
grid on;

figure(4);
plot(misto_rezu, data_mr, misto_rezu, data_mt);
title('Prubeh momentu [N*mm/mm]');
xlabel('r');
ylabel('M(r)');
legend('Mr(r) radialni','Mt(r) obvodovy');
grid on;