%% TOMAS MIVALT
%% MODEL SLOZENEHO TELESA - ANALYTICKY VYPOCET

clc
clear all

%% VSTUPNI HODNOTY

D = 81; % prumer pistu mm   
L = 61.3-9.8; % vyska pistu (delka skorepiny)
hd = 8.45; % tloustka dna pistu mm 
hs = 2.7; % tloustka steny skorepiny 

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
% Vyjadreni pouzivanych hodnot pri vypoctu
% Strednicovy polomer
Rs = D/2-hs/2;
% A maximalni polomer
R = D/2;

% Ohybova tuhost desky - dna pistu
Bd = (E*hd^3)/(12*(1-(u^2)));
% Ohybova tuhost skorepiny - plast pistu
Bs = (E*hs^3)/(12*(1-(u^2)));
% Parametr beta skorepiny - poruseni membranove napjatosti
beta = ((3*(1-(u^2)))/((Rs^2)*(hs^2)))^(1/4);

%% Vyjadreni stykove sily T1 a momentu M1 mezi skorepinou a steno-deskou

% Pozn: rovnice vygenerovany systemem MAPLE dle predchoziho vyjadreni
% deformacnich podminek Usk = Ust a Vsk = Vdesky, v pripade pocitani jineho
% telesa, nez je pist, pripadne jinak zatizeneho je potreba rovnice znovu
% vyjadrit. Viz BP.

M1 = -(1/8)*p*Rs*Bs*beta*(E*hd*hs*Rs^2+E*hd*hs*Rs^2*u-3*E*hd*hs*R^2-E*hd*hs*R^2*u+8*hd*Bd*u^2*R*beta+8*hd*Bd*u*R*beta+2*Rs^3*Bs*beta^3*hs-6*Rs*Bs*beta^3*hs*R^2+4*Rs*Bs*beta^3*hs*R^2*u-2*Rs^3*Bs*beta^3*u^2*hs+2*Rs*Bs*beta^3*u^2*hs*R^2)/(hs*(-E*hd*Bd-E*hd*Bd*u-4*Bd*Rs*Bs*beta^3+4*Bd*u^2*Rs*Bs*beta^3-2*Rs*Bs*beta*E*hd-4*Rs^2*Bs^2*beta^4+4*Rs^2*Bs^2*beta^4*u));
T1 = (1/8)*hd*p*Rs*Bs*beta^2*(E*hs*Rs^2+E*hs*Rs^2*u-3*E*hs*R^2-E*hs*R^2*u+16*Bd*u^2*R*beta+16*Bd*u*R*beta+16*Rs*Bs*beta^2*u*R)/(hs*(-E*hd*Bd-E*hd*Bd*u-4*Bd*Rs*Bs*beta^3+4*Bd*u^2*Rs*Bs*beta^3-2*Rs*Bs*beta*E*hd-4*Rs^2*Bs^2*beta^4+4*Rs^2*Bs^2*beta^4*u));

fprintf('Stykova sila T1=%f [N/mm] , stykovy moment M1=%f [N*mm/mm] \n',T1,M1);

% Funkce jednotlivych prubehu
% Prvne skorepina

% sila, pusobici na skorepinu Nz
% p*pi*R^2=nz*2*pi*R
nz = p*R/2;

% Konstanty skorepiny
skor_c1 = M1/(2*Bs*beta^2);
skor_c2 = -T1/(2*Bs*beta^3)-M1/(2*Bs*beta^2);
% Posuv skorepiny
u_sk = @(z)((exp(-beta*z))*(skor_c1*sin(beta*z)+skor_c2*cos(beta*z))-(Rs*u*p*R)/(E*hs*2));

% Natoceni skorepiny
v_sk = @(z)(-beta*(exp(-beta*z))*(skor_c1*(sin(beta*z)-cos(beta*z))+skor_c2*(sin(beta*z)+cos(beta*z))));

% Moment Mz
mz_sk = @(z)(2*Bs*(beta^2)*exp(-beta*z)*(skor_c1*cos(beta*z)-skor_c2*sin(beta*z)));
% Moment Mt
mt_sk = @(z)(u*mz_sk(z));

% Vysledna liniova sila Nt
nt_sk = @(z)(u*nz+E*hs*u_sk(z)/Rs);

% Napeti skorepiny na vnejsi plose - extrem
sigma_z_sk = @(z)(nz/hs+6*mz_sk(z)/hs^2);
sigma_t_sk = @(z)(nt_sk(z)/hs+6*mt_sk(z)/hs^2);

% DESKA
des_C1 = (-M1/(Bd*(1-u))-(p*(R^2)*(3+u))/(16*Bd*(1+u)));
v_des = @(r)(des_C1*r+(p*r^3)/(16*Bd));

% pruhyb desky
% OP: r=Rs u=0
des_C3 = (-((p*Rs^4)/(64*Bd))-(1/2)*des_C1*(Rs^2));
u_des = @(r)((1/2)*des_C1*(r^2)+((p*Rs^4)/(64*Bd))+des_C3);

% derivace desky
dv_des = @(r)(-M1/(Bd*(1-u))-(p*(R^2)*(3+u))/(16*Bd*(1+u))+(3*p*r^2)/(16*Bd));

% Radialni M
des_mr = @(r)(-Bd*(dv_des(r)+(u*v_des(r)/r)));
% Obvodovy M
des_mt = @(r)(-Bd*((v_des(r)/r)+u*dv_des(r)));

des_sr = @(r)((12*des_mr(r))/(hd^3))*(-hd/2);
des_st = @(r)((12*des_mt(r))/(hd^3))*(-hd/2);
% osove napeti je zanedbatelne oproti sr a st a je mozno psat priblizne -
% horni strana desky
des_sz = -p;

%% Vypis dat
% Rozliseni
n = 3000;

% Skorepina

z_misto = linspace(0,L,n);
nula = 0;

for i=1:n
    
    data_sk_u(i) = u_sk(z_misto(i));
    data_sk_v(i) = v_sk(z_misto(i));
    
    data_sk_mz(i) = mz_sk(z_misto(i));
    data_sk_mt(i) = mt_sk(z_misto(i));
    
    data_sk_sz(i) = sigma_z_sk(z_misto(i));
    data_sk_st(i) = sigma_t_sk(z_misto(i));
    
    
    % HMH
    A = [data_sk_sz(i) data_sk_st(i) nula];
    temp1 = sort(A);
    
    sigma1 = temp1(3);
    sigma2 = temp1(2);
    sigma3 = temp1(1);
    
    data_sk_hmh(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));

end;

% Deska
r_misto = linspace(-R,R,2*n);

for i=1:(2*n)

    data_des_u(i) = u_des(r_misto(i));
    
    data_des_v(i) = v_des(r_misto(i));
     
    data_des_sr(i) = des_sr(r_misto(i));
    data_des_st(i) = des_st(r_misto(i));
    data_des_sz(i) = des_sz;
    
    A = [data_des_sr(i) data_des_st(i) data_des_sz(i)];
    temp1 = sort(A);
    
    sigma1 = temp1(3);
    sigma2 = temp1(2);
    sigma3 = temp1(1);
    
    data_des_hmh(i) = ((2^0.5)/(2))*((((sigma1-sigma2)^2)+((sigma2-sigma3)^2)+((sigma3-sigma1)^2))^(1/2));

end;


%% Kresleni grafu

figure(1);
plot(z_misto, data_sk_u);
title('[SKOŘEPINA] Radiální u [mm]');
xlabel('z [mm]');
ylabel('u(z)');
grid on;
legend('Analyticky', 'Numericky');

figure(2);
plot(z_misto, data_sk_v);
title('[SKOŘEPINA] Natočení skořepiny [rad]');
xlabel('z [mm]');
ylabel('v(z)');
grid on;

figure(3)
plot(z_misto, data_sk_mz, z_misto, data_sk_mt);
title('[SKOŘEPINA] Průběh momentů [N*mm/mm]');
xlabel('z [mm]');
ylabel('M(z)');
legend('Osovy moment Mz','Tečný moment Mt');
grid on;

figure(4)
plot(z_misto, data_sk_sz, z_misto, data_sk_st, z_misto, data_sk_hmh);
title('[SKOŘEPINA] Průběh napětí [MPa]');
xlabel('z [mm]');
ylabel('\sigma(z)');
legend('\sigma_z','\sigma_t','\sigma_{RED}');
grid on;

% Deska
figure(5)
plot(r_misto, data_des_v);
title('[DESKA] Průběh natočení [rad]');
xlabel('r [mm]');
ylabel('v(r)');
grid on;

figure(6)
plot(r_misto, data_des_sr, r_misto, data_des_st,r_misto, data_des_sz, r_misto, data_des_hmh);
title('[DESKA] Průběh napětí [MPa]');
xlabel('r [mm]');
ylabel('\sigma(r)');
legend('Radiální napětí \sigma_z','Tečné napětí \sigma_t','Osové napětí \sigma_z','Redukované napětí \sigma_{RED}');
grid on;

figure(7)
plot(r_misto, data_des_u);
title('[DESKA] Průběh průhybu [mm]');
xlabel('r [mm]');
ylabel('W(r)');
grid on;

