%% 3.3 Flow Analysis and Graphs
% Ford Duratec 2.0 L cylinder head
% SuperFlow SF-110 test data

clear;
clc;
close all;

%% Measured valve data

inletValveDiameter_mm = 35.1;
exhaustValveDiameter_mm = 30.1;

inletValveLift_mm = ...
    [1.755 3.510 5.270 7.020 8.780 10.530 12.290 14.040];

inletFlow_Lps = ...
    [28.552 56.457 78.774 96.740 109.178 113.324 117.470 117.470];

exhaustValveLift_mm = ...
    [1.500 3.010 4.520 6.020 7.530 9.030 ...
     10.540 12.040 13.550];

exhaustFlow_Lps = ...
    [17.55975 41.625 57.000 66.000 70.5375 ...
     72.765 74.250 74.32425 75.735];

%% Quarter-diameter valve lifts

inletQuarterLift_mm = inletValveDiameter_mm / 4;
exhaustQuarterLift_mm = exhaustValveDiameter_mm / 4;

%% Valve-head areas

% Convert valve diameters from mm to cm
inletValveDiameter_cm = inletValveDiameter_mm / 10;
exhaustValveDiameter_cm = exhaustValveDiameter_mm / 10;

inletValveArea_cm2 = pi / 4 * inletValveDiameter_cm^2;
exhaustValveArea_cm2 = pi / 4 * exhaustValveDiameter_cm^2;

%% Flow rate per unit valve-head area

inletFlowPerArea = inletFlow_Lps / inletValveArea_cm2;
exhaustFlowPerArea = exhaustFlow_Lps / exhaustValveArea_cm2;

%% Figure 12 - Flow rate versus valve lift

figure;

plot(inletValveLift_mm,inletFlow_Lps, ...
    '-o','LineWidth',2,'MarkerSize',6);

hold on;

plot(exhaustValveLift_mm,exhaustFlow_Lps, ...
    '-s','LineWidth',2,'MarkerSize',6);

xline(inletQuarterLift_mm,'--', ...
    sprintf('Inlet D/4 = %.2f mm',inletQuarterLift_mm), ...
    'LabelVerticalAlignment','middle');

xline(exhaustQuarterLift_mm,'--', ...
    sprintf('Exhaust D/4 = %.2f mm',exhaustQuarterLift_mm), ...
    'LabelVerticalAlignment','bottom');

grid on;
xlabel('Valve Lift (mm)');
ylabel('Volumetric Flow Rate (L/s)');
title('Measured Flow Rate vs Valve Lift');
legend('Inlet','Exhaust','Location','southeast');

xlim([0 15]);
ylim([0 130]);

%% Figure 13 - Flow rate per unit valve area

figure;

plot(inletValveLift_mm,inletFlowPerArea, ...
    '-o','LineWidth',2,'MarkerSize',6);

hold on;

plot(exhaustValveLift_mm,exhaustFlowPerArea, ...
    '-s','LineWidth',2,'MarkerSize',6);

xline(inletQuarterLift_mm,'--', ...
    sprintf('Inlet D/4 = %.2f mm',inletQuarterLift_mm));

xline(exhaustQuarterLift_mm,'--', ...
    sprintf('Exhaust D/4 = %.2f mm',exhaustQuarterLift_mm));

grid on;
xlabel('Valve Lift (mm)');
ylabel('Flow Rate per Valve Area (L/s/cm^2)');
title('Area-Normalised Flow Rate vs Valve Lift');
legend('Inlet','Exhaust','Location','southeast');

xlim([0 15]);

%% Flow coefficient calculation

% Assumption:
% SuperFlow test pressure = 25 inches of water

testPressure_inH2O = 25;

% Convert inches of water to pascals
inH2O_to_Pa = 249.0889;
deltaP_Pa = testPressure_inH2O * inH2O_to_Pa;

% Standard air density
airDensity = 1.225;       % kg/m^3

%% Effective valve curtain area

% Curtain area = circumference x valve lift
%
% A_curtain = pi x valve diameter x valve lift

inletDiameter_m = inletValveDiameter_mm / 1000;
exhaustDiameter_m = exhaustValveDiameter_mm / 1000;

inletLift_m = inletValveLift_mm / 1000;
exhaustLift_m = exhaustValveLift_mm / 1000;

inletCurtainArea_m2 = ...
    pi .* inletDiameter_m .* inletLift_m;

exhaustCurtainArea_m2 = ...
    pi .* exhaustDiameter_m .* exhaustLift_m;

%% Ideal volumetric flow

% Q_ideal = A x sqrt(2*deltaP/rho)

idealVelocity = sqrt(2 * deltaP_Pa / airDensity);

inletIdealFlow_m3s = inletCurtainArea_m2 .* idealVelocity;
exhaustIdealFlow_m3s = exhaustCurtainArea_m2 .* idealVelocity;

% Convert measured flow from L/s to m^3/s
inletMeasuredFlow_m3s = inletFlow_Lps / 1000;
exhaustMeasuredFlow_m3s = exhaustFlow_Lps / 1000;

%% Flow coefficients

inletCF = inletMeasuredFlow_m3s ./ inletIdealFlow_m3s;
exhaustCF = exhaustMeasuredFlow_m3s ./ exhaustIdealFlow_m3s;

%% Figure 14 - Flow coefficient versus valve lift

figure;

plot(inletValveLift_mm,inletCF, ...
    '-o','LineWidth',2,'MarkerSize',6);

hold on;

plot(exhaustValveLift_mm,exhaustCF, ...
    '-s','LineWidth',2,'MarkerSize',6);

xline(inletQuarterLift_mm,'--', ...
    sprintf('Inlet D/4 = %.2f mm',inletQuarterLift_mm));

xline(exhaustQuarterLift_mm,'--', ...
    sprintf('Exhaust D/4 = %.2f mm',exhaustQuarterLift_mm));

grid on;
xlabel('Valve Lift (mm)');
ylabel('Flow Coefficient, C_F');
title('Flow Coefficient vs Valve Lift');
legend('Inlet','Exhaust','Location','best');

xlim([0 15]);

%% Summary results

fprintf('FLOW ANALYSIS RESULTS\n');
fprintf('---------------------\n');

fprintf('Inlet valve area: %.3f cm^2\n',inletValveArea_cm2);
fprintf('Exhaust valve area: %.3f cm^2\n',exhaustValveArea_cm2);

fprintf('Inlet quarter-diameter lift: %.3f mm\n', ...
    inletQuarterLift_mm);

fprintf('Exhaust quarter-diameter lift: %.3f mm\n', ...
    exhaustQuarterLift_mm);

fprintf('Maximum inlet flow: %.3f L/s\n', ...
    max(inletFlow_Lps));

fprintf('Maximum exhaust flow: %.3f L/s\n', ...
    max(exhaustFlow_Lps));

fprintf('Maximum inlet flow coefficient: %.3f\n', ...
    max(inletCF));

fprintf('Maximum exhaust flow coefficient: %.3f\n', ...
    max(exhaustCF));

fprintf('Mean inlet flow coefficient: %.3f\n', ...
    mean(inletCF));

fprintf('Mean exhaust flow coefficient: %.3f\n', ...
    mean(exhaustCF));

fprintf('Maximum-flow exhaust/intake ratio: %.1f %%\n', ...
    max(exhaustFlow_Lps)/max(inletFlow_Lps)*100);

%% Results table

inletResults = table( ...
    inletValveLift_mm', ...
    inletFlow_Lps', ...
    inletFlowPerArea', ...
    inletCF', ...
    'VariableNames', ...
    {'Lift_mm','Flow_Lps','FlowPerArea_Lps_cm2','FlowCoefficient'});

exhaustResults = table( ...
    exhaustValveLift_mm', ...
    exhaustFlow_Lps', ...
    exhaustFlowPerArea', ...
    exhaustCF', ...
    'VariableNames', ...
    {'Lift_mm','Flow_Lps','FlowPerArea_Lps_cm2','FlowCoefficient'});

disp('Inlet results');
disp(inletResults);

disp('Exhaust results');
disp(exhaustResults);