%
%  Calculation of Control Valve Characteristics
%  and pressure flow curves
%
%  Bob Newell, February 1996
%
%-------------------------------------------------------
% Copy file and enter your data here
%
Psource = 450;           % source pressure or maximum
                         % pump delivery pressure (kPa)
Rpump = 0.8;             % pump internal resistance
                         % (kPa.s.s/kg.kg)
Psink = 250;             % sink pressure (kPa)
Rpipe = 2.4;             % pipe and fitting resistance
                         % (kPa.s.s/kg.kg)
Dflow = 5;               % design flowrate (kg/s)
Dmaxflow = 6.2;          % maximum design flowrate (kg/s)
Dminflow = 3.5;          % minimum design flowrate (kg/s)
inherent = 3;            % 1=linear inherent characteristic,
                         % 2=square root inherent char.,
                         % 3=equal percentage inherent char.
%
%--------------------------------------------------------
% Pressure - Flow Curve
%
mmax = sqrt( ( Psource - Psink ) / ( Rpump + Rpipe ) );
disp( [ 'Maximum flowrate = ' num2str( mmax ) ] )
if mmax < Dmaxflow,
  disp( [ 'WARNING - insufficient pressure drop ' ...
               'for maximum flow' ] );
end,
cmargin = fix( 100 * mmax / Dmaxflow );
disp( [ 'Control Margin = ' num2str(cmargin-100) ' %' ] )
if cmargin < 115,
  disp( [ 'WARNING - insufficient control margin ' ...
               'at maximum flow' ] );
end,
if cmargin > 115,
  cmargin = 115;
end;
%
dpmax = Psource - Psink - Rpump * Dmaxflow * Dmaxflow;
dp = Psource - Psink - ( Rpump + Rpipe ) ...
                               * Dmaxflow * Dmaxflow;
pc = 100 * dp / dpmax;
disp( [ 'DelP at maximum design flowrate = ' ...
      num2str( dpmax ) ' ( ' num2str(fix(pc)) '% )' ] )
if pc < 25,
  disp( [ 'WARNING - insufficient controllability ' ...
                 'at maximum design flowrate' ] );
end,
if pc > 50,
  disp( [ 'WARNING - excessive pressure drop ' ...
                 'at maximum design flowrate' ] );
end,
%
m = 0:mmax/100:mmax;
figure
plot( m, Psink * ones( size( m ) ), ...
  m, Psource * ones( size( m ) ) - Rpump * m .* m, ...
  m, Psink * ones( size( m ) ) + Rpipe * m .* m );
xlabel( 'Flow ( kg/s )' )
ylabel( 'Pressure ( kPa )' )
grid
hold on
ax = axis;
Pinc = ( ax(4) - ax(3) ) / 10;
plot( [ Dflow Dflow ], [ ax(3) ax(4) ], '--' )
plot( [ Dminflow Dminflow ], [ ax(3)  ax(4) ], '--' )
plot( [ Dmaxflow Dmaxflow ], [ ax(3)  ax(4) ], '--' )
hold off
%
%--------------------------------------------------------
% Valve Sizing Calculation
%
flow = cmargin * Dmaxflow / 100;
dp = Psource - Psink - ( Rpump + Rpipe ) ...
			                    * flow * flow;
cvmax = flow / sqrt( dp );
disp( [ 'Cvmax calculated = ' num2str(cvmax) ] )
cvsel = input( 'Selected Cvmax ? ' );
if length(cvsel) > 0,
  cvmax = cvsel;
end;
fx = 1;
dp = ( Psource - Psink ) / ( 1 + ( Rpump + Rpipe ) * ...
            cvmax * cvmax * fx .* fx );
mmax = cvmax * fx .* sqrt( dp );
disp( [ 'Maximum flowrate = ' num2str(mmax) ] )
%--------------------------------------------------------
% Valve Characteristics
%
% valve stem position
x=0:0.01:1;
if inherent == 1,
  fx = x;
elseif inherent == 2,
  fx = sqrt( x );
elseif inherent == 3,
  a = input( 'Equal Percentage design constant [30] ? ' );
  if length(a) == 0,
    a = 30;
  end;
  fx = a .^ ( x - 1 );
else
  disp( 'ERROR - unknown characteristic' )
end;
dp = ( Psource - Psink ) * ones( size( fx ) ) ...
         ./ ( ones( size( fx ) ) + ( Rpump + Rpipe ) * ...
            cvmax * cvmax * fx .* fx );
m = cvmax * fx .* sqrt( dp );
fxin = m ./ mmax;
%
figure
plot( x, fx, '--', x, fxin, '-' )
legend( 'Inherent', 'Installed' );
xlabel('Valve stem position')
ylabel('Valve characteristic')
grid
if inherent > 1,
  hold on
  line( [ 0 1 ], [ 0 1 ] )
  hold off
end;
%-----------------------------------------------------------
% the end
