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

clc
clear all

%% VSTUPNI HODNOTY PISTU
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 vetknuta -> Okrajove Podminky:
% r = R ... W=0 .. ve vazbe posuv 0
% r = R ... V=0 .. ve vazbe natoceni 0

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

C1 = -((p*(R^2))/(16*B));
C3 = +((p*(R^4))/(64*B));

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

%% Radialni a obvodovy moment
% Radialni M
Mr = @(r)((-p/16)*((r^2)*(3+u)-(R^2)*(1+u)));
% Obvodovy M
Mt = @(r)((-p/16)*((r^2)*(1+3*u)-(R^2)*(1+u)));


%% 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

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

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

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


% SIGMA red 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));
    
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 napětí [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;