%%-------------------------------------------------------------------------
% Lubrificazione
%%-------------------------------------------------------------------------
% soluzione analitica delle eq. NS
%
% p = (6*mu*W)/(alpha*(Y0+Yz)*Y^2)*(-Y^2+(Y0+Yz)*Y-Y0*Yz);
% dpdz = (6*mu*W)/(Y.^2)*(1-(2*Y0*Yz)/((Y0+Yz)*Y));
% w=W*(1-y/Y)-1/(2*mu)*dpdz*y*(Y-y);
% sigmayz=-mu*W/Y+dpdz*(y-Y/2);
%%-------------------------------------------------------------------------
% Data: 2018-10-16
% Autore: Lucia Parussini
%%-------------------------------------------------------------------------

clear all;
close all;
clc;

% parametri del flusso
rho=1; %kg/m^3 
nu=1e-3; %m^2/s 
mu=rho*nu; 
W=0.1; %m/s
Z=5e-1; %m 
Y0=1e-2; %m
Yz=1.5e-3; %m


alpha=(Y0-Yz)/Z; % inclinazione della lastra superiore
Re= W*Z/nu; % Reynolds number

display(strcat(['inclinazione ',num2str(atan(alpha)*180/pi),'°']));
display(strcat(['Re ',num2str(Re)]));

%% plottaggio geometria
figure
fill([0 Z Z 0],[0 0 Yz Y0],'b');
axis equal
title('Geometry');
xlabel('z [m]');
ylabel('y [m]');

%% plottaggio andamento pressione e gradiente di pressione
% punti in cui calcolare pressione e gradiente di pressione
z=linspace(0,Z,100);

% altezza del meato nei punti z
Y=Y0-alpha*z;

% pressioe e gradiente di pressione nei punti z
p = (6*mu*W)./(alpha*(Y0+Yz)*Y.^2).*(-Y.^2+(Y0+Yz)*Y-Y0*Yz);
dpdz = (6*mu*W)./(Y.^2).*(1-(2*Y0*Yz)./((Y0+Yz)*Y));

% grafico
figure;
subplot(2,1,1)
plot(z,p);
grid on;
title('pressure')
xlabel('z [m]');
ylabel('p [Pa]');
subplot(2,1,2)
plot(z,dpdz);
grid on;
title('pressure gradient')
xlabel('z [m]');
ylabel('dp [Pa/m]');

zmax=Y0*Z/(Y0+Yz);
display(strcat(['z di massima pressione ',num2str(zmax),' m']));

%% profili di velocità e tensioni tangenziali su 11 sezioni + quella di massima pressione
zz=linspace(0,Z,11); % posizione delle sezioni
zz=sort([zz,zmax]);
YY=Y0-alpha*zz; % altezza del meato nelle sezioni
dpdz = (6*mu*W)./(YY.^2).*(1-(2*Y0*Yz)./((Y0+Yz)*YY)); % gradiente di pressione nelle sezioni

for i=1:length(zz) % calcolo di velocità e tensione in ciascuna sezione
    yy(i,:)=linspace(0,YY(i),20);
    w(i,:)=W*(1-yy(i,:)/YY(i))-1/(2*mu)*dpdz(i)*yy(i,:).*(YY(i)-yy(i,:)); 
    sigmayz(i,:)=-mu*W/YY(i)+dpdz(i)*(yy(i,:)-YY(i)/2);
end

ZZ=mat2cell(zz,ones(size(zz,1),1),ones(size(zz,2),1));
labels=cellfun(@(s)sprintf('%.2f',s),ZZ,'UniformOutput',false);

% grafico della velocità 
% ho usato un fattore di scala Lw per w per motivi di visualizzazione
Lw=1;
figure;
hold on;
grid on;
line([0 Z Z 0 0], [0 0 Yz Y0 0],'Color',[0 0 0],'LineWidth',2); %geometria
for i=1:length(zz)
    c=rand(1,3);
    plot(w(i,:)/Lw+zz(i),yy(i,:),'Color',c);
    plot([zz(i),zz(i)],[0,YY(i)],'Color',c)
    for j=1:length(yy(i,:)) % traccio le linee orrizontali del profilo
        plot([zz(i),w(i,j)/Lw+zz(i)], [yy(i,j),yy(i,j)],'Color',c)
    end
end
xticks(zz);
xticklabels(labels);
% axis equal
    
%grafico delle tensioni
% fattore di scala Ls per sigmayz per motivi di visualizzazione
Ls=1e0;
figure;
hold on;
grid on;
line([0 Z Z 0 0], [0 0 Yz Y0 0],'Color',[0 0 0],'LineWidth',2); %geometria
for i=1:length(zz)
    c=rand(1,3);
    plot(sigmayz(i,:)/Ls+zz(i),yy(i,:),'Color',c);
    plot([zz(i),zz(i)],[0,YY(i)],'Color',c)
    for j=1:length(yy(i,:))% traccio le linee orrizontali del profilo
        plot([zz(i),sigmayz(i,j)/Ls+zz(i)], [yy(i,j),yy(i,j)],'Color',c)
    end
end
xticks(zz);
xticklabels(labels);
% axis equal