% Input data
x = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160]; % D
H_elevation = [50, 120, 200, 250, 320, 400, 250, 200, 180, 210, 270, 325, 375, 425,
% Pump characteristics
Q = [0, 100, 200, 300, 400, 500, 600, 700]; % Flow rate (m^3/hr)
H_pump = [525, 520, 510, 495, 475, 450, 420, 385]; % Head (m)
% Convert flow rate to m^3/s
Q_m3s = Q / 3600; % Convert m^3/hr to m^3/s
% Two pumps in series: combined head
H_combined = 2 * H_pump;
% Pipeline and crude oil properties
D_o = 0.324; % Outside diameter (m)
epsilon = 0.06e-3; % Roughness (m)
L1 = 70e3; % Length from 0 to 70 km (m)
L2 = 90e3; % Length from 70 to 160 km (m)
nu = 10e-6; % Kinematic viscosity (m^2/s)
SG = 0.86; % Specific gravity
gamma = SG * 9.81; % Specific weight (N/m^3)
g = 9.81; % Acceleration due to gravity (m/s^2)
% Calculate static head
H_static1 = H_elevation(8) - H_elevation(1); % Static head from 0 to 70 km
H_static2 = H_elevation(end) - H_elevation(8); % Static head from 70 to 160 km
% Friction factor and k values (use Darcy-Weisbach)
f = 0.02; % Assumed friction factor (can be refined with Moody chart later)
D_i = D_o - 2 * 0.006; % Internal diameter (subtract wall thickness)
k1 = 0.8 * f * L1 / (g * D_i^5);
k2 = 0.8 * f * L2 / (g * D_i^5);
% System curves
H_system1 = @(Q) H_static1 + k1 * Q.^2; % System curve for 0-70 km
H_system2 = @(Q) H_static2 + k2 * Q.^2; % System curve for 70-160 km
% Generate flow range for intersection search
Q_range = linspace(0, max(Q_m3s), 1000); % Fine flow range for interpolation
H_pump_interpolated = interp1(Q_m3s, H_combined, Q_range); % Pump curve
H_sys1_values = H_system1(Q_range); % System curve 1 values
H_sys2_values = H_system2(Q_range); % System curve 2 values
% Find intersection points
[~, idx1] = min(abs(H_pump_interpolated - H_sys1_values)); % Closest match for 0-70
[~, idx2] = min(abs(H_pump_interpolated - H_sys2_values)); % Closest match for 70-1
Q_intersection1 = Q_range(idx1); % Flow rate at intersection for 0-70 km
Q_intersection2 = Q_range(idx2); % Flow rate at intersection for 70-160 km
H_intersection1 = H_pump_interpolated(idx1); % Head at intersection for 0-70 km
H_intersection2 = H_pump_interpolated(idx2); % Head at intersection for 70-160 km
1
% Convert intersection points from m^3/s to m^3/hour
Q_intersection1_m3h = Q_intersection1 * 3600; % Flow rate in m^3/hour for 0-70 km
Q_intersection2_m3h = Q_intersection2 * 3600; % Flow rate in m^3/hour for 70-160 km
% Display intersection results in m^3/hour
disp('Intersection Points:');
disp(['0-70 km: Flow Rate (Q) = ', num2str(Q_intersection1_m3h), ' m^3/hour, Head (
disp(['70-160 km: Flow Rate (Q) = ', num2str(Q_intersection2_m3h), ' m^3/hour, Head
% Plot pump and system curves
figure;
plot(Q_range * 3600, H_pump_interpolated, 'r-', 'LineWidth', 1.5); hold on; % Pump
plot(Q_range * 3600, H_sys1_values, 'b--', 'LineWidth', 1.5); % System curve 0-70 k
plot(Q_range * 3600, H_sys2_values, 'g--', 'LineWidth', 1.5); % System curve 70-160
xlabel('Flow Rate (m^3/hour)');
ylabel('Head (m)');
legend('Pump Curve (2 Pumps in Series)', 'System Curve (0-70 km)', 'System Curve (7
title('Pump and System Curves');
grid on;
% Mark intersection points on the plot
plot(Q_intersection1_m3h, H_intersection1, 'ko', 'MarkerSize', 10, 'MarkerFaceColor
plot(Q_intersection2_m3h, H_intersection2, 'mo', 'MarkerSize', 10, 'MarkerFaceColor
text(Q_intersection1_m3h, H_intersection1, ...
sprintf(' Q = %.1f m^3/hour, H = %.1f m', Q_intersection1_m3h, H_intersection1
'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
text(Q_intersection2_m3h, H_intersection2, ...
sprintf(' Q = %.1f m^3/hour, H = %.1f m', Q_intersection2_m3h, H_intersection2
'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
% Hydraulic gradient points calculation
% Internal diameter
D_i = D_o - 2 * 0.006; % Internal diameter (subtract wall thickness)
A = pi * (D_i / 2)^2; % Cross-sectional area (m^2)
% Velocity and losses for 0-70 km
v1 = Q_intersection1 / A; % Velocity for 0-70 km
losses1 = f * (L1 / D_i) * (v1^2 / (2 * g)); % Head loss for 0-70 km
% Velocity and losses for 70-160 km
v2 = Q_intersection2 / A; % Velocity for 70-160 km
losses2 = f * (L2 / D_i) * (v2^2 / (2 * g)); % Head loss for 70-160 km
% Hydraulic gradient points
H1 = 50 + H_intersection1; % Point at X(0)
H2 = H1 - losses1 + 200; % Point at X(70)
H3 = H2 + H_intersection2; % Point at X(70) after second pump
H4 = 100 + H3 - losses2; % Point at X(160)
% Elevation profile
x = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160]; % D
H_elevation = [50, 120, 200, 250, 320, 400, 250, 200, 180, 210, 270, 325, 375, 425,
% Hydraulic gradient line
x_hg = [0, 70, 70, 160]; % X-coordinates for hydraulic gradient
2
H_hg = [H1, H2, H3, H4]; % Y-coordinates for hydraulic gradient
% Plot elevation and hydraulic gradient
figure;
plot(x, H_elevation, 'b-', 'LineWidth', 1.5); hold on;
plot(x_hg, H_hg, 'r--', 'LineWidth', 1.5);
scatter(x_hg, H_hg, 50, 'ro', 'filled'); % Mark hydraulic gradient points
xlabel('Distance (km)');
ylabel('Height/Elevation (m)');
legend('Elevation Profile', 'Hydraulic Gradient Line', 'Location', 'Best');
title('Hydraulic Gradient and Elevation Profile');
grid on;
% Annotate points
text(x_hg, H_hg, ...
arrayfun(@(h) sprintf(' %.1f m', h), H_hg, 'UniformOutput', false), ...
'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');
% Given data
utilization_factor = 0.913; % Utilization factor
rho = 860; % Density of crude oil (kg/m^3)
g = 9.81; % Gravitational acceleration (m/s^2)
% Pump data
Q_pump = [0, 100, 200, 300, 400, 500, 600, 700]; % Flow rates in m^3/hour
Efficiency = [0, 36, 60, 74, 80, 78, 70, 58]; % Efficiency in %
% Convert Q_intersection2 to m^3/hour for throughput calculation
Q_intersection2_m3h = Q_intersection2 * 3600; % Flow rate in m^3/hour
% Annual throughput
annual_throughput = Q_intersection2_m3h * 24 * 365 * utilization_factor; % m^3/year
% Interpolated pump efficiency at Q_intersection2
efficiency2 = interp1(Q_pump, Efficiency / 100, Q_intersection2 * 3600, 'linear', '
% Absorbed horsepower at each station
H_absorbed_station1 = H_intersection1; % Head at intersection
H_absorbed_station2 = H_intersection2; % Head at intersection
HP_station1 = (rho * g * H_absorbed_station1 * Q_intersection1) / (efficiency2 * 10
HP_station2 = (rho * g * H_absorbed_station2 * Q_intersection2) / (efficiency2 * 10
% Display results
disp('Results:');
disp(['Annual Throughput: ', num2str(annual_throughput, '%.2f'), ' m^3/year']);
disp(['Absorbed Horsepower at Station 1: ', num2str(HP_station1, '%.2f'), ' HP']);
disp(['Absorbed Horsepower at Station 2: ', num2str(HP_station2, '%.2f'), ' HP']);
Intersection Points:
0-70 km: Flow Rate (Q) = 501.7017 m^3/hour, Head (H) = 898.979 m
70-160 km: Flow Rate (Q) = 509.4094 m^3/hour, Head (H) = 894.3544 m
Results:
Annual Throughput: 4074195.33 m^3/year
Absorbed Horsepower at Station 1: 1368.29 HP
Absorbed Horsepower at Station 2: 1382.16 HP
3
4
Published with MATLAB® R2014a
