%% SYDNEY UNIVERSITY ENGINEERING UNDERGRADUATE ASSOCIATION KEG CALCULATOR
% Author: Jack Naylor
% Description: The following uses some elementary iterative calculations
% derived in MECH3260 and MECH3261 to calculate the flow and heat transfer
% characteristics of SUEUA's "Magic Box" or Keg Dispensing System.

%% SETUP
clear
clc
close all

%% INPUTS
disp('ENSURE YOU SALT THE ICE IN THE MAGIC BOX!')
P_keg_gauge = input('CO2 Pressure (PSI) = ')*6.895;
T_keg = input('Approximate Keg Temperature (deg C) = ');
T_air = input('Ambient Air Temperature (deg C) = ');
percentFull = input('What percentage full is the keg? ')/100;

%% CONSTANTS
e_ss = 0.015e-3;   % Stainless steel roughness factor
e_pvc = 0.0015e-3; % PVC roughness factor

K_reentrant = 0.78; % Reentrant loss
eq_length_serve = 150; % Equivalent loss of angle valve serving
eq_length_ninety = 60;
K_coupling = 0.08; % Coupling losses



cp_w = 4.217; % kj/kgK
rho_w = 1000; % kg/m^3
nu_w = 1/rho_w;
kf_w = 569e-3;
Prf_w = 12.99; 
mu_w = 1750e-6; % Viscosity
alpha_w = kf_w/(rho_w*cp_w);
nu_w2 = mu_w/rho_w;
cp_beer = 4.05; % kj/kgK
rho_beer = 1020; % Density of beer kg/m^3
nu_beer = 1/rho_beer; % Specific volume beer m^3/kg
mu_beer = mu_w;
alpha_beer = kf_w/(rho_beer*cp_beer);


L_ss = 21; % Length of coil m
OD_ss = 9.52e-3; % Outer diameter m
t_ss = 0.6e-3;  % Tubing thickness m
ID_ss = OD_ss - 2*t_ss; % Inner diameter m

% Stainless steel
rho_ss = 8055;
cp_ss = 480;
k_ss = 15.1;
alpha_H_ss = 3.91e-6;

% PVC
k_pvc = 0.18716;

tol = 1e-6;

P_atm = 101.325;
%P_keg_gauge = 82.7371;
P_keg = P_atm+P_keg_gauge;
g = 9.81;

z1 = 0;
z2 = 1.5;

d_ex = 5e-3;
A_ex = pi*d_ex^2/4;

d_pvc = 6e-3;
A_pvc = pi*d_pvc^2/4;

d_coup = 16e-3;
A_coup = pi*d_coup^2/4;

d_coil = ID_ss;
A_coil = pi*d_coil^2/4;

d_barb = 4e-3;
A_barb = pi*d_coil^2/4;

keg_height = 0.593725;
L_spear = 4/5*keg_height;
d_spear = 18e-3;
A_spear = pi*d_spear^2/4;

L_pvc = 1.5;

% CO2 properties @ 280K
cp_co2 = 0.83;
mu_co2 = 140e-7;
nu_co2 = 7.36e-6;
k_co2 = 15.2e-3;
alpha_co2 = 9.63e-6;
Pr_co2 = 0.765;


% Air properties @ 300K
cp_air = 1.007;
mu_air = 184.6e-7;
nu_air = 15.89e-6;
k_air = 26.3e-3;
alpha_air = 22.5e-6;
Pr_air = 0.707;

K_phi = (1+1.68*(d_coil/0.205)^0.68);

R_co2 = 0.1889;
V_keg = 80e-3;
V_co2_init = 0.05*V_keg;

d_keg = 0.409575;

%% SOLVE

% Solve for flow rate
f_ss = 0.02;
f_pvc = 0.02;
f_ss_spear = 0.02;
alpha_ss = 1;
alpha_pvc = 1;
alpha_ss_spear = 1;
syms fNew

% mDot = sqrt(P_keg-P_atm+2*g*(z1-z2))/(1/rho_beer*sqrt(A_ex^2+...
%     3*A_coup^2*K_coupling+A_spear^2*K_reentrant+A_ex^2*f_ss*eq_length_serve...
%     +alpha_pvc*A_pvc^2*f_pvc*L_pvc/d_pvc+A_barb^2*f_ss*eq_length_ninety...
%     +alpha_ss*f_ss*L_ss*A_coil^2/d_coil + alpha_ss_spear*A_spear^2*f_ss_spear*L_spear/d_spear));

mDotPrev = 1;
difference = 5;


% While the difference > O(10^(-5))
while difference>=tol
    
    % Solve energy equation using assumed values
    mDot = sqrt(P_keg-P_atm+2*g*(z1-z2))/(sqrt((1/A_ex^2+3*1/A_coup^2*K_coupling...
            +1/A_spear^2*K_reentrant+1/A_ex^2*f_ss*eq_length_serve+alpha_pvc*1/A_pvc^2*...
            f_pvc*L_pvc/d_pvc+1/A_barb^2*f_ss*eq_length_ninety+alpha_ss*f_ss*K_phi*L_ss*1/A_coil^2/d_coil...
            + alpha_ss_spear*1/A_spear^2*f_ss_spear*L_spear/d_spear+1/A_coil^2*(f_ss*eq_length_ninety*(20*4)*2))/rho_beer^2));
    
    
    % Limit solution to be positive and real
    mDot = mDot(mDot>0&imag(mDot)==0);
   
    % Calculate Reynold's number of flow
    Re_ss = 4*mDot/(pi*mu_beer*d_coil);
    Re_pvc = 4*mDot/(pi*mu_beer*d_pvc);
    Re_ss_spear = 4*mDot/(pi*mu_beer*d_spear);
    
    % Solve for friction factor using solved diameter and Re number
    if alpha_ss == 1
        fSolved1 = vpa(solve(1/sqrt(fNew)==-2*log10((e_ss/d_coil)/3.7+2.51/(Re_ss*sqrt(fNew))),fNew));
    else
        fSolved1 = 64/Re_ss;
    end
    
    if alpha_pvc == 1
        fSolved2 = vpa(solve(1/sqrt(fNew)==-2*log10((e_pvc/d_pvc)/3.7+2.51/(Re_pvc*sqrt(fNew))),fNew));
    else
        fSolved2 = 64/Re_pvc;
    end
    
    if alpha_ss_spear == 1
        fSolved3 = vpa(solve(1/sqrt(fNew)==-2*log10((e_ss/d_spear)/3.7+2.51/(Re_ss_spear*sqrt(fNew))),fNew));
    else
        fSolved3 = 64/Re_ss_spear;
    end
    % Set this as new f value
    f_ss = fSolved1;
    f_ss_spear = fSolved3;
    f_pvc = fSolved2;
    
    % Check difference as condition of convergence
    difference = abs(mDot - mDotPrev);
    
    % Set solution for next loop and print to console
    mDotPrev = mDot;
    
    if Re_ss < 2300 && alpha_ss == 1
        alpha_ss = 2;
        difference = 1;
        d = 0;
        
    end
    
    if Re_pvc < 2300 && alpha_pvc == 1
        alpha_pvc = 2;
        difference = 1;
        d = 0;
        
    end
    
    if Re_ss_spear < 2300 && alpha_ss_spear == 1
        alpha_ss_spear = 2;
        difference = 1;
        d = 0;
        
    end
end

disp(['Beer flows at ' num2str(double(mDot)) ' kg/s'])
disp(['It takes ' num2str(double(532e-6/(mDot*nu_beer))) ' seconds to fill a red cup.'])
%% THERMALS
T_ice = -0.5;
T_beer = T_keg;
%T_air = 30;

T_co2 = (T_keg+T_air)/2;
% Heat transfer to fluid from re-entrant, heat transfer to CO2, convection
% in spear, convection in PVC, conduction in coil, natural covection to ice


% Assume 50/50 air fluid in Keg - No self transfer!!!!! - this can be
% adjusted but as trials have shown, there is little heat transfer to the
% CO2.
%Nu = 3.66; % Constant surface temperature, laminar flow

if Re_ss_spear > 2300
    Nu_spear = 0.023*double(Re_ss_spear)^(4/5)*Prf_w^0.3;
else
    Nu_spear = 3.66;
end
h_spear_i = Nu_spear*kf_w/(L_spear*percentFull);
if Re_pvc > 2300
    Nu_pvc = 0.023*double(Re_pvc)^(4/5)*Prf_w^0.3;
else
    Nu_pvc = 3.66;
end
h_pvc_i = Nu_pvc*kf_w/L_pvc;
if Re_ss > 2300
    Nu_ss = 0.023*double(Re_ss)^(4/5)*Prf_w^0.3;
else
    Nu_ss = 3.66;
end
h_ss_i = Nu_ss*kf_w/L_ss;

T_so = T_co2;
t_diff = 5;
while t_diff >tol
    Tf = 273+(T_so+T_co2)/2;
    RaD = abs(g*1/Tf*(T_so-T_co2)*(d_spear/4)^3/(nu_co2*alpha_co2));
    Nu_spear_o = (0.825+0.387*RaD^(1/6)/(1+(0.492/Pr_co2)^(9/16))^(8/27))^2;
    h_spear_o = k_co2*Nu_spear_o/(L_spear*percentFull);
    
    R_tot = 1/(h_spear_i*pi*d_spear*L_spear/2)+1/(h_spear_o*pi*d_spear*L_spear/2)+log((d_spear+0.5e-3)/d_spear)/(2*pi*k_ss*L_spear/2);
    T_out_spear = double(T_co2-(T_co2-T_beer)*exp(-(1/R_tot)/(mDot*cp_beer)));
    
    qPrime_tot = double((T_beer-T_out_spear)/(R_tot*L_spear*percentFull));
    
    T_soN = double(T_co2+qPrime_tot*1/(h_spear_o*pi*d_spear));
    t_diff = double(abs(T_soN-T_so));
    T_so = T_soN;
end 
t_diff = 5;
T_so = T_air;
while t_diff >tol
    Tf = 273+(T_so+T_air)/2;
    RaD = abs(g*1/Tf*(T_so-T_air)*(d_pvc)^3/(nu_air*alpha_air));
    Nu_pvc_o = (0.6+0.387*RaD^(1/6)/(1+(0.559/Pr_air)^(9/16))^(8/27))^2;
    h_pvc_o = k_air*Nu_pvc_o/(L_pvc);
    
    R_tot = 1/(h_pvc_i*pi*d_pvc*L_pvc)+1/(h_pvc_o*pi*d_pvc*L_pvc)+log((d_pvc+2e-3)/d_pvc)/(2*pi*k_pvc*L_pvc);
    T_out_pvc = double(T_air-(T_air-T_out_spear)*exp(-(1/R_tot)/(mDot*cp_beer)));
    
    qPrime_tot = double((T_out_spear-T_out_pvc)/(R_tot*L_pvc));
    
    T_soN = double(T_air+qPrime_tot*1/(h_pvc_o*pi*d_pvc));
    t_diff = double(abs(T_soN-T_so));
    T_so = T_soN;
end 
t_diff = 5;
T_so = T_ice;
while t_diff >tol
    Tf = 273+(T_so+T_ice)/2;
    RaD = abs(g*1/Tf*(T_so-T_ice)*(d_coil)^3/(nu_w2*alpha_w));
    Nu_ss_o = (0.6+0.387*RaD^(1/6)/(1+(0.559/Prf_w)^(9/16))^(8/27))^2;
    h_ss_o = kf_w*Nu_ss_o/(L_ss);
    
    R_tot = 1/(h_ss_i*pi*d_coil*L_ss)+1/(h_ss_o*pi*d_coil*L_ss)+log(OD_ss/d_coil)/(2*pi*k_ss*L_ss);
    T_out_coil = double(T_ice-(T_ice-T_out_pvc)*exp(-(1/R_tot)/(mDot*cp_beer)));
    
    qPrime_tot = double((T_out_pvc-T_out_coil)/(R_tot*L_ss));
    
    T_soN = double(T_ice+qPrime_tot*1/(h_ss_o*pi*d_coil));
    t_diff = double(abs(T_soN-T_so));
    T_so = T_soN;
end 

disp(['Pour temperature = ' num2str(T_out_coil) ' deg C'])
if T_out_coil < 9 && T_out_coil > 6
    disp('This is optimum temperature for enjoyment!')
elseif T_out_coil < 6
    disp('The temperature is too low, you risk freezing in the pipes!')
else
    disp('The temperature is too high, nobody likes warm beer!')
end