%created by Khalid Kane 
%Fluid flow analytical model %Hydrostatic model
%-------------------------------------------------------------------------- 
%Fluid parameters in SI units 
g = 9.81; %g is the gravitation field in m/s2
ro = 1000; %density of water in kg/m3
mu = 0.6913*10^-3; %viscosity of water at 3°C in Pa.s %mu = 1.3*10^-3; 
%%viscosity of water at 10°C in Pa.s
%-------------------------------------------------------------------------- 
%Add more wells to the medium inlet and outlet in equal proportions, max 10 
%wells could be added on each sides 
%Change k(a positive integer from 0 to 10 wells, 0 being the current design format) 
k=0; 
%-------------------------------------------------------------------------- 
%Medium lane parameters in SI units
we1_i=0.0045*(k+1); 
%loading well width (assume it is a  square) in m
we2_i=0.0045; 
%loading well width (assume it is a  square) in m
l_i = (9.2*10^-3)+(k*2*we2_i); 
%length of  medium channel in m
%-------------------------------------------------------------------------- 
%Medium lane parameters in SI units
w_i = 200*10^-6; 
%width of the perfusion lane in m
%-------------------------------------------------------------------------- 
d1_i = 120*10^-6; 
%height of the channel in both medium and gel channels in m
%--------------------------------------------------------------------------
vl_i=d1_i*((w_i*l_i)); 
%total volume of the perfusion lane in terms of space available in m3
%-------------------------------------------------------------------------- 
%media loaded into the medium inlet and outlet, can play with these values 
%to keep the initial pressure constant while increasing the loading wells 
v1_i=80*10^-9*(1+k); 
%loaded volume in the first well in m3
v2_i=40*10^-9*(1+k); 
%loaded volume in the second well in m3
ca_i=we1_i*we2_i; 
%loading well cross-section area
z1_i=(v1_i-vl_i)/ca_i; 
%Initial height in first well
z2_i=(v2_i)/ca_i; 
%Intial height in second well
h_i=z1_i-z2_i; 
%intial driving height
r_i=(12*mu*l_i)/((w_i*(d1_i^3))*(1-((0.630*d1_i)/w_i)));
%hydrodynamic resistance rectangular cross-section  
%width w and height h, expression valid when h<w)
%Applied pressure
p_i=ro*g*h_i;
dt=1; 
%s is the time step can be adjusted if we want more accuracy, but will take longer to compute 
for n=1:86401    
t(n)=n*dt;    
q_i(n)=p_i(n)/r_i;     
dh_i=q_i(n)*dt/ca_i;     
p_i(n+1)=p_i(n)-2*ro*g*dh_i; 
%equivalent flow rate 
end

%Change k(a positive integer from 0 to 10 wells, 0 being the current design format) 
k=1; 
%-------------------------------------------------------------------------- 
%Medium lane parameters in SI units
we1=0.0045*(k+1); 
%loading well width (assume it is a  square) in m
we2=0.0045; 
%loading well width (assume it is a  square) in m
l_p=4.395*10^-3;
%length of phaseguide
l_r = (4.805*10^-3)+(k*2*we2); 
%length of  medium channel in m
l=l_p+l_r;
%-------------------------------------------------------------------------- 
%Medium lane parameters in SI units
w = 100*10^-6; 
%width of the perfusion lane (phaseguide section) in m
%-------------------------------------------------------------------------- 
d1 = 50*10^-6; 
%height of the channel in medium lan (phaseguide section) in m
%--------------------------------------------------------------------------
w_r = 100*10^-6; 
%width of the perfusion lane (phaseguide section) in m
%-------------------------------------------------------------------------- 
d1_r = 50*10^-6; 
%height of the channel in medium lan (phaseguide section) in m
%--------------------------------------------------------------------------
r_p=(12*mu*l_p)/((w*(d1^3))*(1-((0.630*d1)/w)));
%hydraulic resistance of phaseguide section of medium lane
r_r=(12*mu*l_r)/((w_r*(d1_r^3))*(1-((0.630*d1_r)/w_r)));
%hydraulic resistance for the remaining section of the medium lane
r=r_p+r_r;
%total hydraulic resistance

vl=d1*((w*l)); 
%total volume of the perfusion lane in terms of space available in m3
%-------------------------------------------------------------------------- 
%media loaded into the medium inlet and outlet, can play with these values 
%to keep the initial pressure constant while increasing the loading wells 
v1=118*10^-9*(1+k); 
%loaded volume in the first well in m3
v2=10*10^-9*(1+k); 
%loaded volume in the second well in m3
ca=we1*we2; 
%loading well cross-section area
z1=(v1-vl)/ca; 
%Initial height in first well
z2=(v2)/ca; 
%Intial height in second well
h=z1-z2; 
%intial driving height
% r=(12*mu*l)/((w*(d1^3))*(1-((0.630*d1)/w)));
%hydrodynamic resistance rectangular cross-section  
%width w and height h, expression valid when h<w)
p=ro*g*h; 
%Applied pressure
dt=1; 
%s is the time step can be adjusted if we want more accuracy, but will take longer to compute 
for n=1:86401   
t(n)=n*dt;    
q(n)=p(n)/r;     
dh=q(n)*dt/ca;     
p(n+1)=p(n)-2*ro*g*dh; 
%equivalent flow rate 
end
graph=plot(t/3600,q_i*1000*1000*1000*3600,t/3600,q*1000*1000*1000*3600);
graph(1).LineWidth = 1;
graph(2).LineWidth = 1;
xlabel('Time (hour)','FontSize',18,'FontWeight','bold','Color','black'); 
ylabel('Average flow rate (microlitre/hour)','FontSize',18,'FontWeight','bold','Color','black'); ; 
title('Average flow rate of the media in the medium lane','FontSize',20,'FontWeight','bold','Color','r');
legend('Old design','New Design')
graph=plot(t/3600,q*1000*1000*1000*3600);
graph(1).LineWidth = 1;
xlabel('Time (hour)','FontSize',18,'FontWeight','bold','Color','black'); 
ylabel('Average flow rate (microlitre/hour)','FontSize',18,'FontWeight','bold','Color','black');  
%title('Average flow rate of the media in the medium lane','FontSize',20,'FontWeight','bold','Color','r');