HW 7 Problem 2(a)

Q2_NonlinearApprox.m

clear; clc; cla; clf; close all; format long;

global Ccb Chp Chl Cre zB;

% Physical constants
g   = 9.81;         % m / s^2
rho = 1000;         % kg / m^2
nu  = 0.9e-6;       % m^2 / s

% Given properties
D       = 100e-3;   % m
L       = 50;       % m
zB      = 12;       % m
epsilon	= 0.3e-3;   % m
P_pump	= 4.5e3;    % W

% Derived properties
A = pi/4 * D^2;     % m^2

% Precomputed constants
Ccb = epsilon / D / 3.72;   % Colebrook constant term
Chp = P_pump / (g*rho*A);   % Pump head constant factor
Chl = L / (D * 2*g);        % Head loss constant factor
Cre = D / nu / 2.51;        % Reynolds number as used in Colebrook

% Initial point
f0 = 0.0038;  % ASU ID: xxxxxxxx 38
U0 = 1;       % m/s

x = num2cell(fsolve(@pipeSystem, [f0, U0]));
[f,U] = x{:}


function F = pipeSystem(x)
  global Ccb Chp Chl Cre zB;
  
  f = x(1);
  U = x(2);

  Re = Cre * U;
  
  % Colebrook equation
  sqrt_f_i = 1 / sqrt(f);
  F(1) = sqrt_f_i + 2*log10(Ccb + sqrt_f_i/Re);
  
  % Pump head, head loss, and gravitational potential
  F(2) = (Chp / U) - (f * Chl * U^2) - zB;
end