% 
%  Esercitazione sulla lubrifcazione del meato a facce piane
%   Meato a spessore variabile linearmente
%
close  all
clear all
%
%  Problema INVERSO
%
%  Dati del problema
%
l = 0.4  % m
V = 2.4 % m/s
h1 = 0.40 * 10^(-3)  % m spessore massimo del meato (di ingresso)
h2 = 0.16*10^(-3)  % m spessore minimo del meato (in uscita)
mu = 40* 10^(-3) % Pa s  viscosità alla temperatura di funzionamento
n = h1/h2; % coefficiente caratteristico del cuscinetto
%
%  Deduzione del diagramma dei gradienti
%
x = 0:0.01:l;
zeri = x * 0;
h= h1 - x * (h1-h2)/(l); 
pp  =  6 * mu * V *    ( h +  h*n - 2*h2*n  )   ./   (    (1+n)* h.^ 3  );
%
%
p = -6 * l * mu * V * (h2 ^ 2 * n - h2 .* h - h2 .* h * n + h .^ 2) / h2 .^ 2 ./ h .^ 2 / (n ^ 2 - 1);
%
%
%  Valore numerico di N
%
N = 0.6e1 * l ^ 2 / (h2 * n - h2) * mu * V * (-(2 * n) + 0.2e1 + log(n) + n * log(n)) / h2 / (n ^ 2 - 1)
%
%
%  Valore numerico di xN
%
%
xN = l / (h2 * n - h2) / (-(2 * n) + 0.2e1 + log(n) + n * log(n)) * h2 * (-(5 * n ^ 2) + (4 * n) + 0.1e1 + 0.2e1 * (n ^ 2) * log(n) + 0.4e1 * n * log(n)) / 0.2e1
%
%
%  Valore del coefficiente di attrito mediato
%
fm = 0.1e1 / l * (-(3 * n) + 0.2e1 * n * log(n) + 0.2e1 * log(n) + 0.3e1) * (h2 * n - h2) / (-(2 * n) + 0.2e1 + log(n) + n * log(n)) / 0.3e1
%
%  diagramma delle velocità u delle falde fluide parallele
%
xd = 0.15;
hxd = h1 - xd * (h1-h2)/(l); 
accur = 50;
deltaz = h2/accur;
z = 0:deltaz:hxd;
u = -V * z .* (2 * hxd.^ 2 + 2 * hxd ^ 2 * n - 3 * hxd * z - 3 * hxd * n * z - 6 * h2 * hxd * n + 6 * h2 * n * z) / (1 + n) / hxd^ 3;
%
%
scala_pp = 1;
scala_p = 15;
figure(1)
set(1,'Color','w');
plot(x,zeri,'-k');
hold on;
plot(x, pp *scala_pp,'-r','LineWidth',1);
hold on
plot(x,p*scala_p,'-g','LineWidth',1);
%hold on
% plot(ep,k*scala_k,'-.b','LineWidth',1);
%
title(" Dyagrams ");
ylabel("  p gradients ");
xlabel("  x axis  ");
% legend ("ascisse","gradienti","pressioni", "location","outside");
grid "minor" on;
hold off
%
%
figure(2)
set(2,'Color','w');
plot(u,z,'-r');
%
%  ---------PROBLEMA DINAMICO DIRETTO------------------------
%
%
P = 80000  % N  carico totale
omega = 800 * 2 * pi / 60  % radianti
rm = 100 * 10^(-3)  % m  raggio medio
l = 73.5 * 10^(-3) % m   lunghezza del pattino
b = 60 * 10^(-3)  % m   larghezza del pattino
x_N = 45 * 10^(-3)  % m   applicazione del pattino noto
mu = 35 * 10^(-3)  % Pa s   larghezza del pattino
%
%
%  Ecco come trovare n coefficiente caratteristico del cuscinetto
%
%  Definisco la funzione
%
xN_funz  =  @(nn)    (( l / (h2 * nn - h2) ./ (-(2 * nn) + 0.2e1 + log(nn) + nn * log(nn)) * h2 * (-(5 * nn .^ 2) + (4 * nn) + 0.1e1 + 0.2e1 * (nn .^ 2) * log(nn) + 0.4e1 * nn .* log(nn)) / 0.2e1 )   - x_N  )
%
%  Verifica
%
nab = 1:0.05:5;
for i = 1:length(nab)
	nasx = nab(i);
	xN_verif(i) = xN_funz(nasx);
end;
figure(3)
plot(nab,xN_verif)
hold on
plot(nab,nab*0,'-y')
grid "minor" on
%
%  Usa la funzione fsolve per trovare lo zero della funzione
%  ...noi ne facciamo davvero un uso banale, in realtà 
%     fsolve è in grado di fare molto di più...
%
%  Impongo la soluzione di primo tentativo
%
tentativo = 4.5
%
soluzione = fsolve(xN_funz, tentativo)
%
% In alternativa possiamo procedere con Newton Raphson
% per trovare lo zero
%
%  Ora completa tu l'esercitazione ....