AN INVESTIGATION INTO THE EFFICIENT
NUMERICAL INTEGRATION OF STIFF
HYDRAULIC SYSTEMS CONTAINING
DISCONTINUITIES
By
CRAlG EDWARD WENZEL
Bachelor ofScience
Milwaukee School of Engineering
Milwaukee, Wisconsin
1990
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
December, 1999
AN INVESTIGATIONlNTO THE EFFICIENT
NUMERICAL INTEGRATION OF STIFF
HYDRAULIC SYSTEMS CONTAINING
DISCONTINUITIES
Thesis Approved:
Dea of the Graduate College
11
ACKNOWLEDGMENT
I wish to express my sincere gratitude to my principle advisor Dr. 1. T. Hong, for
his patience, expert guidance, and support during my thesis effort and throughout the
course ofmy graduate study. My sincere appreciation also extends to Dr. R. K.
Tessmann whose insightful suggestions and encouragement were invaluable. In addition
I wish to thank my committee members, Dr. Delahoussaye and Dr. Pagilla, for their time
and recommendations. I also wish to thank FES\Bardyne Inc., Stillwater, Oklahoma, for
providing excellent laboratory facilities and computing equipment. Particular
appreciation is given to the faculty and staff of the Department of Mechanical and
Aerospace Engineering for providing a quality educational experience. Finally, I wish to
express my gratitude to my family and friends for their support throughout all ofmy
academic pursuits.
III
TABLE OF CONTENTS
Chapter Pag
1. WTRODUCTION 1
Problem Statement. 1
Objective of Study 3
2. LITERATURE REVIEW 5
Introduction 5
Mathematical Stiffness 5
Numerical Integration of Stiff ODEs 12
Discontinuities and StiffODE Solvers 14
Summary , 17
3. MATHEMATICAL MODELS AND CUBIC SPLINES t 9
Introduction t9
Background 19
TwoPoint Cubic Splines 20
FourPoint Cubic SpLines , 24
Comparison of TwoPoint and Four Point Splines 29
Application Examples 30
4. THE EVENT SWITCHING ALGORITHM 45
Background 45
Hydraulic System State Events 46
Development of the Event Switching Algorithm 51
Application Example 56
5. CASE STUDY: PrLOTOPERATED RELIEF VALVE 67
Introduction 67
Overview ofStitfODE Solvers 67
The PilotOperated Relief Valve 69
Numerical Integration Evaluation 75
Discussion of Computation Time Results 83
Error Analysis 85
Discussion of Error Analysis Results 87
Summary 88
IV
Chapter Page
6. EXPERIMENTAL VERIFICATION 89
Introduction 89
Test Component. 89
Test System 91
Mathematical Models and Differential Equations 92
Results of Computer Analysis 92
Test Results 93
Discussion ofResults 95
7. CONCLUSION AND RECOMMENDATrONS 96
Conclusion 96
Suggestions for Further Study 98
Works Cited ]00
Appendix A: Procedure for Generating Cubic Splines to Interpolate
Between Four Points 102
Appendix B: Experimental Determination of the Pressure and
Flow Relationship ofan Orifice 105
Appendix C: MATLAB Script File for Example 3. ] 109
Appendix D: MATLAB Script File for Example 3.2 I 13
Appendix E: MATLAB Script File for Example 4. 1.. 120
Appendix F: MATLAB Script Files for PilotOperated Rehef Valve
Case Study from Chapter 5 127
LIST OF TABLES
Table Page
5.1 Simulation Results tor Example 4.1 Using Different
Variations of Gear's Method , 64
5.2 Compiled Test Results for Variations of the PilotOperated
Relief Valve Problem , , 83
5.3 Error Analysis Results for Inlet Pressure Response ofTest #4 86
B1 Test Results for 0.012 Diameter Orifice 108
vi
Figure
2.1
2.2
2.3
2.4
2.5
2.6
3.1
3.2
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
LIST OF FIGURES
Page
Orifice PressurelFlow Relationship Using Turbulent Flow Equation 7
Fluid Volume Discharge to the Atmosphere through an Orifice 7
Orifice Flow Model with Linear Laminar Flow Region 9
Orifice Flow Model with Cubic Spline Laminar Flow Region 11
Volume Changes for an Extending Hydraulic Cylinder.. 12
Cubic Spline Used to Smooth a Discontinuity 15
Cubic Spline Interpolation of Data Points 20
Two Functions Requiring a Smooth Connecting Cubic Spline 20
Cubic Spline Providing Slope Continuity Between Two Functions 21
Effect of Moving a Single Data Point (A to A') 25
Cubic Spline End Points and their Control Points 26
Cubic Spline Generated with fourPoint Form 27
Flow Chart for Generating Control Points 28
Comparison ofTwoPoint and FourPoint Cubic Splines 30
Experimental Results of Orifice Test. 3 I
Comparison ofTest Results and Turbulent Flow Equation 32
Cubic Spline Based Orifice ModeL 34
Refined Orifice Model. 35
3.13 Typical Pressure ReliefYalve 39
3.14 General Static Relief Valve Model. 40
3.15 Smoothed Static Relief Valve ModeL ..fJ
4.1 SpringMassDamper System with Mechanical Stops 46
4.2 Three States for SpringMassDamper System with Mechanical Stops .47
4.3 Flow Chart for Event Switching Algorithm as Applied to
General SpringMassDamper System 55
vii
Figure
4.4
4.5
4.6
4.7
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.1 I
6.1
6.2
6.3
6.4
6.5
B1
Page
Hydraulic Circuit Schematic for Example 4.1 56
Cylinder Rod Displacement Response for Example 4.1 65
Cylinder Rod Velocity Response for Example 4.1.. 65
Pressure Responses for Example 4.1 66
Pilot0perated Relief Valve 69
PilotOperated Relief Valve Test Circuit.. 75
Pressure Response for Test #] 79
Displacement Response for Test #1 79
Pressure Response for Test #2 80
Displacement Response for Test #2 80
Pressure Response for Test #3 81
Displacement Response for Test #3 8]
Pressure Response for Test #4 82
Displacement Response for Test #4 82
Comparison of Error Analysis Results 87
Typical DirectActing Pressure Relief Valve 90
DirectActing Relief Valve Test Circuit. 91
Simulation Results for DirectActing Relief Valve 93
Test Results for DirectActing ReliefValve 94
Simulation and Test Results Superimposed 94
Hydraulic System for Orifice Testing 106
VlIJ
A
b
b
BDF
Cd
D
ESA
f
F
Fpreload
h
k
m
m
n
NDF
NOMENCLATURE
Area
Independent Axis Intercept of General Straight Line
Damping Ratio
Backward Differentiation Formula
Discharge Coefficient
Diameter
Event Switching Algorithm
Generic Function Notation
Force
Spring Preload Force
Interval Size
Spring Constant
Flow Pressure Coefficient
Lumped Orifice Coefficient
Slope of General Straight Line
Mass
Integer Value Notation
Numerical Differentiation Formula
IX
NRr Transition Reynolds Number
ODE Ordinary Differential Equation
P Pressure
Per Relief Valve Cracking Pressure
Pmax Relief Valve Maximum Pressure
Q Volumetric Flow Rate
Q,?Per Relief Valve Flow Rate at Per
Qmax Relief Valve Flow Rate at Pmar
RK RungeKutta Method
S Second Derivative (Curvature)
SffiK SemiImplicit RungeKutta Method
t Time
V Fluid Volume
w Weight
a Poppet HalfAngle
f3 Bulk Modulus ofElasticity
8 Laminar Flow Coefficient
t1P Differential Pressure
r1Ptf Differential Pressure for Fully Turbulent Flow
E Computational Precision
A. Increment Factor
v Kinematic Fluid Viscosity
x
p Fluid Density
Time Constant
Discontinuity Function
Xl
CHAPTER 1
INTRODUCTION
Problem Statement
Computer modeling and simulation techniques are important to designers in every
engineering discipline. These techniques are particularly valuable to engineers
concerned with the dynamic response of physical systems. A number of numerical
integration methods are available to solve the ordinary differential equations (ODEs)
which describe a dynamical system. In the computer environment the effect of design
changes on dynamic response can be determined without the financial consequences of
complicated "breadboard" testing. However, the time required to simulate the dynamic
response of a system does have some inherent costs. In addition, the results ofthese
simulations are only as accurate as the mathematical models used to represent the
physical system.
As with engineers in many other fields, fluid power system desif,rners have turned
to computer modeling and simulation to predict the dynamic response characteristics of
hydraulic circuits. A hydraulic circuit uses pressurized fluid to provide a desired
mechanical output. In addition, the various pressure and flow control components of a
hydraulic system are mechanical in nature. A typical hydraulic circuit, therefore, is
comprised ofcomponents with widely varying time constants. As a result, the
differential equations which describe this combination of hydraulics and mechanics are
often numerically stiff. Small time constants force the use of relatively small time steps
during numerical integration. Because the differential equations are linked, the presence
2
of large time constants requires a great number ofthese small steps to capture the overall
response. As a result, computation time becomes a concern. In the past complex
hydraulic circuits often required hours or days to complete a single simulati.on.
Computation time must be reduced to enhance the value of computer analysi .
A number of numerical integration methods have been developed to decrease the
computation time required to solve stiff differential equations. Gear's method for stiff
ODEs and Rosenbrock's semiimplicit RungeKutta (SIRK) method are the most
prominent among the stiff solvers[l]. The current stateoftheart for both solvers
features an adaptive step size algorithm to further improve computation time. Both of
these methods require the calculation of Jacobian matrix for each time step. This
Jacobian matrix is not required for standard explicit solvers like Euler's method and the
RungeKutta (RK) method. The introduction ofthe Jacobian matrix creates a
compatibility problem with hydraulic systems.
Physical constraints in a hydraulic system often cause stiff ODE solvers to fail.
These constraints may be fluidrelated or mechanical. For example, the turbulent flow
equation has an infinite slope at the origin. This infinite slope will produce an illconditioned
Jacobian matrix when a stiffODE solver is used. Similarly, a limitation on
the displacement of a mechanical component will produce a mathematical discontinuity
in the Jacobian matrix during the integration process. These discontinuities often cause a
stiff ODE solver to fail or cause the adaptive step size routine to slow the solver
dramatically.
A few attempts have been made to apply stiffODE solv rs to th umqu
boundary constraints of hydraulic systems. Two different approaches have b n
proposed to eliminate the infinite stope at the origin ofthe turbulent flow equation.
Ellman considers both laminar and turbulent orifice discharge coefficients with a cubic
spline used to create a mathematicaJJy smooth transition between the two types of
flow[2]. Piche', Ellman, and Vilenius also used a cubic spline approach in dealing
directly with the laminar and turbulent flow equations[l]. Both of these methods are
limited to orifice flow and cannot be adapted to model other flow related discontinuities.
Bowns, Tomlinson and Dugdale have suggested a numerical integration technique to deal
with mechanical discontinuities[3]. This technique involves step size halving and
restarting procedures for the ODE solver when a discontinuity is encountered. However,
the proposed method is limited to fixed step ODE solvers and requires a modification to
the core solver algorithm.
Objective of Study
The purpose of this study is to develop a unique procedure that eliminate the
numerical integration problems caused by the fluidrelated and mechanical boundary
conditions found in hydraulic systems. A unique curve fitting technique which uses
control points is used to create mathematically smooth fluid flow models. The concept
of control points involves positioning a point in such close proximity to the end point of a
curve that the difference is physicallly negligible. However, the relative position of the
end point and the control point determines the shape of the resulting curve. As such, the
control point is physically insignificant but mathematically important. The versatility of
4
this technique allows hydraulic engineers to create a greater variety of fluid flow model.
More importantly, the continuity and smoothness of the resulting mathematical models
provide compatibility with stiffODE solvers and significantly reduce computation time.
Numerical integration problems caused b mechanical boundary condition will
be eliminated by an Event SWitch;n~ AI~orithm (ESA). The ESA tracks the displacement
of any given mass. If a physical limit is encountered by the mass, the initial conditions
are reset and a revised set ofODEs are solved numerically. The sum ofthe forces acting
on the mass is monitored while the mass is, at the physical limit. If this sum of forces
becomes unbalanced, the algorithm reverts back to the original set of ODEs and the
solver is automatically restarted. By continuing in this manner numerical integration
proceeds unimpeded by the mathematical discontinuities associated with mechanical
boundaries. In addition, the ESA requires no modification to the core ODE solver code.
As a result, it is compatible with commercially available numerical integration software.
CHAPTER 2
LITERATURE REVIEW
Introduction
As with any dynamic system, hydraulic circuits are mathematically d fined by a
series ofordinary differential equations (ODEs). Unlike most systems, however,
hydraulic systems are generally dominated by nonlinearities [3]. These highly nonlinear
ODEs are, for all practical purposes, impossible to solve with conventional
analytical techniques r4]. Computeri.zed numerical integration is the only practical
approach to solving a set ofODEs accurately describing a hydraulic system.
Unfortunately, the numerical integration of hydraulic systems is a notoriously slow
process. The nature of hydraulic systems requires relatively small time steps and a great
deal of computational effort.
Mathematical Stiffness
Hydraulic systems have been identified as mathematically stiflby a number of
researchers [1,3,4,5,6,7,8]. Bowns et. al. [3] and Krus [6] define stiffness as a set of
ODEs containing widely varying time constants. In other words, the system contains
both rapidly and slowly varying transient solutions [8]. The effect of mathematical
stiffness on numerical integration is well documented [1,3,4,5]. Piche' and Ellman [5]
summarized this effect for hydraulic systems by stating, "Conventional explicit numerical
integration methods such as classical RungeKutta schemes become numerically unstable
unless a very small time increment is used, which leads to excessively long computation
times".
Ellman [2] identified two sources of stiffness in hydraulic systems. The first
source is stiffuess inherent 1n turbulent flow through an orifice. An orifice 1S the most
basic element in hydraulic control. Pressure and flow control in a hydraulic system are
performed by components which use fixed and variable orifices. Orifice flow is
dominated by turbulence and, therefore, laminar flow is typically ignored. Turbulent
flow through an orifice is governed by the following relationship:
6
(2.1)
Q= Flow Through the Orifice
Cd = Discharge Coefficient ofthe Orifice
A = Flow Area of the Orifice
p = Fluid Density
t1P = Pressure Drop Across the Orifice
Discharge coefficient, flow area, and fluid density are normally consi.dered to be
constant. As such, the flow through an orifice is proportional to the square root of the
pressure drop across it. A graph ofthis relationship for an arbitrary orifice is shown in
Figure 2.1. The slope of the curve approaches infinity as the pressure drop approaches
zero. This characteristic is responsible for mathematical stiffness. An example
presented by Krus and Palmberg [7] effectively illustrates the problem. A volume of
o
Turbulent Flow
Equation
Pressure
7
Figure 2. 1: Orifice PressurelFlow Relationship Using Turbulent Flow Equation
fluid is connected to the environment through an orifice as depicted in Figure 2.2. The
0 .( J ~
/'..... ~ Pressure Fluid Orifice
Source Volume Reservoir
Figure 2.2: Fluid Volume Discharge to the Atmosphere through an Orifice
time constant for emptying the volume is described by the following equation:
8
(2.2)
r= Time Constant
v = Fluid Volume
p= Bulk Modulus of the Fluid
Kc = qj = Flow Pressure Coefficient ofthe Orifice
cP
If the turbulent orifice flow equation (2.1) is used, the tenn Kc becomes large when the
pressure drop is small. The time constant, in turn, becomes small. In a numerical
integration situation, relatively small time constants will lead to mathematical stiffness
and long computation times. If the pressure drop is zero, Kc becomes infinite and the
time constant goes to zero. A zero time constant creates infinite stiffness and causes all
classical numerical integration methods to fail [4].
A second source of stiffness identified by Ellman [2] involves widely varying oil
volumes within a single system. The elasticity of pressurized hydraulic fluid is
dependent upon the volume ofthe trapped fluid. Flow passages within hydraulic control
valves are orders of magnitude smaller than the hoses, tubes, housings, actuators, etc.
used to contain and utilize the maj'ority of the working fluid. The low elasticity of the
fluid contained in these small passages is directly related to small time constants. Small
time constants among relatively large time constants within a single system lead to
mathematical stiffness.
Various researchers have developed methods to eliminate the stiffuess problem
caused by the turbulent orifice equation (2.1). Technically the turbulent orifice equation
is not appropriate at low pressure differences because the flow through an orifice is
predominantly laminar under these conditions. This fact is usually ignored by designers
because very low pressure drops seldom occur in steady state analysis. Dynamically,
however, it is advantageous to introduce a laminar flow mechanism at low differential
pressures to eliminate the infinite slope problem. Bowns, Tomlinson, and Dorey [8]
proposed a small linear region about the origin as depicted in Figure 2.3. This region
eliminates the infinite slope at the origin and accurately models the linearity of laminar
flow. However, the authors admitted an inherent difficulty in detennining the width of
the laminar range.
4 Turbulent Flow
I Region
I
(1)
;a~
~o LL
oI I
~Laminar
I IFlow Region
Turbulent Flow
Equation
Pressure
Figure 2.3: Orifice Flow Model with Linear Laminar Flow Re!:,rlon
10
If the region is not wide enough, mathematical stiffitess results. An excessively wide
linear region will lead to gross inaccuracy.
Based on studies conducted at the University of Bath [8], an improved laminar
flow mechanism was developed at Tempere University of Technology in Finland.
Ellman and Vilenius [2,9] proposed the use of a cubic spline to model laminar flow and
the region of transition between laminar and turbulent flow. This model was later
refined by Piche' and EHman [5]. The resulting polynomial is presented below:
(2.3)
where,
IF Kinematic Viscosity ofthe Fluid
D = Diameter ofthe Orifice
NR1 = Transition Reynolds Number
&>/f= Pressure Difference for Fully Turbulent Flow
The transition to the turbulent flow equation (2. I) occurs at U /fas defined below:
(2.4 )
The above laminar/transition flow model for an arbitrary orifice is depicted graphically in
Figure 2.4. Piche' suggests a transition Reynolds number (NR1) of 1000. This model
eliminates the infinite slope at the origin and provides slope continuity at the transition
11
Turbulent Flow
Equation
':. Turbulent Flow
I Region
I
Cubic
Spline
0,1: Pressure
~Laminar Flow
I IRegion
(1)
+o~
~o LL
Figure 2.4: Orifice Flow Model with Cubic Spline Laminar Flow Region
point (L1P1/)'
Bowns and Wang [4] proposed an incompressible flow model to eliminate the
stiffness problem caused by small oil volumes. Pressure changes occur very rapidly
under small volume conditions. The incompressible flow model considers these pressure
changes to be instantaneous. In the absence of compressibility, the flow into a
pressurized volume is equal to the flow exiting. Bowns and Wang implemented the
incompressible flow model with an iterative approach based on this flow continuity
property.
Several researchers pointed out an inherent problem with the incompressible flow
model [1,5,6]. The dynamic nature of hydraulic systems requires oil volumes to change
during the course of simulation. As a result the incompressible flo mod I is not valid
at all times. For example, the volume of trapped oil within a hydraulic c Iinder will
increase as the actuator extends (see Figure 2.5). The complexity ofthe model must be
L
VI)} :(V2
V1t [J g=;V2
V1f 'ITc
V2
Figure 2.5: Volume Changes for an Extending Hydraulic Cylinder
increased to accommodate varying volumes by switching between compressible and
incompressible flow models. As stated by Piche' and Ellman [5]," The problem remain
as to how and when to make these transitions smoothly and automatically."
Numerical Integration of StiffODEs
In the absence of a sound modeling technique to eliminate stiffness created by
widely varying oil volumes, researchers have investigated various numerical integration
methods to improve computation time. Simulation efficiency can be improved by using
a numerical ODE solver specifically designed for stiff systems of equations. Further
efficiency gains may be realized by employing an adaptive step size algorithm.
Gear's method for solving stiff ODEs has been generally accepted as the most
efficient algorithm for simulating hydraulic systems [3,8]. The stability features ofthis
13
ODE solver are superior to those of classic nonstiff solvers like RungeKutta.
MathematicaUy, Gear's method for stiff ODEs is A table. Astabili means the
numerical integration method is stable for any step size as long as the set ofODEs is
stable [I]. This added stability allows the solver to use a much larger step size and
therefore, reduces processing time. Recently, Piche' and Ellman [1,5] proposed an Lstable
RungeKutta method for use with fluid power systems. This method provides an
even higher degree of stability for use with extremely stiff systems. Although Gear's
method is Astable, the degree of stability becomes small when the ODE system is
extremely stiff and problems arise in the form ofnumerical oscillations. An Lstable
method drives modal amplification to zero as the time constant approaches zero. As a
result, the Lstable RungeKutta method eliminates numerical oscillations associated
with certain class of hydraulic systems.
The efficiency gains realized by Gear's method for stiffODEs and the Lstable
RungeKutta method do not come without a price. Both ofthe methods require a
significant amount of computational effort. A general system of ODEs may be written as
follows.
(2.5)
In order to achieve improved stability, the numerical integration method must have some
knowledge offat each step [6]. This information is stored in a Jacobian matrix which
contains partial derivatives of/with respect to the state variables. The presence ofthis
Jacobian reintroduces the term CQ. If the turbulent flow equation (2.1) is used as an
iP
orifice model, the trend toward an infinite slope at the origin will produce an il1
14
conditioned Jacobian matrix at low pressure differences. This problem may be solved b
using a laminar flow mechanism hke the one developed by Piche' and Ellman (Eq .2.3
and 2.4).
To further improve processing time, an adaptive step size algorithm may be used
in conjunction with a stiffODE solver. The adaptive step size algorithm automatically
adjusts the step length as integration proceeds. The size of the step is controlled by a
preset error limit. Piche' and Ellman [5] states that adaptive step size control essentially
eliminates numerical stability problems because the algorithm automatically selects a
step length small enough to give an accurate solution. Improved efficiency is achieved
because a small step size is used only when necessary. Relatively large steps may used
after the dynamics associated with small time constants are damped out.
Discontinuities and StiffODE Solvers
The physical nature of hydraulic systems, unfortunately, is not directly
compatible with stiffODE solvers or adaptive step size algorithms. Numerical
integration problems arise in the form ofdiscontinuities. Discontinuities affect both the
core ODE solver and the adaptive step size algorithm. Abrupt changes to elements in the
Jacobian matrix may produce a convergence problem within the integration routine and
cause the method to fail [1). An adaptive step size algorithm will "hunt" around the
discontinuity until the step size is reduced sufficiently to cross the discontinuity within
the preset error limit. This hunting involves a large number of unsuccessful function
evaluations and results in considerable processing time [10].
Discontinuities may be present in the classical mathematical mod I ofa ph ical
entity. A simple solution to this problem is to change the mathematical mod I. Bowns,
Tomlinson, and Dorey [8] have made extensive use of cubic splines to create smooth,
continuous models. A generic example is depicted in Figure 2.6. The onginaI model
f(X)
B C
x
Figure 2.6 Cubic Spline Used to Smooth a Discontinuity
contains a hard nonlinearity at point B. A cubic spline connected between points A and
C is used to provide a smooth transition across the discontinuity. Some amount of
accuracy is sacrificed to improve computational efficiency by allowing a more gradual
change to the Jacobian matrix.
Discontinuities may also be encountered during the simulation process. Bown
Tomlinson, and Dugdale [8) identified two types of discontinuities likely to cause
problems for an ODE solver:
1) Discontinuities which occur at a known time.
2) Discontinuities which occur when a variable reaches a critical value.
Ofthese, the second discontinuity type is more difficult to handle because the exact time
of threshold crossing is unknown. Ellison [11] has labeled this type of discontinuity as a
slale event. A state event involves a switching function which changes terms in the
original set ofODEs and defines a new integration problem starting exactly at the
switching point [10]. State events in hydraulic systems occur when actuators, loads, or
internal valve parts encounter mechanical travel limits. For example, the velocity and
acceleration of a hydraulic cylinder rod almost Instantaneously drop to zero when the
stroke limit is reached. Chaney [12) has recommended restarting the ODE solver when a
state event is encountered. This procedure divides the original problem into continuous
sections and solves them in a piecewise manner. Preston and Berzins [13] later endorsed
this restarting procedure as absolutely necessary "whenever parts of a network uddenly
become (in)active."
A difficulty arises in locating the exact time at which the critical value is reached.
The critical value invariably falls between two successive time steps. Special measures
must be taken to locate the time ofthe discontinuity, within acceptable error limits.
before integration proceeds with a new set of ODEs. It is first necessary to identify a
discontinuity by defining a discontinuity function. In practice, discontinuities are located
17
by monitoring sign changes [14]. For a set ofODEs ofthe fonni i =!(X"X2tX3....X,,,f). a
discontinuity function has the form:
¢II = !(X"X2,X3, ...X",I) (1.6)
where the discontinuity occurs at,
(1.7)
In order to locate the precise time of a zero crossing, a discontinuity handling algorithm
is necessary. The goal of this algorithm is to control the step size such that the
discontinuity occurs at the end of the step [14]. Chaney [12] proposed an iterative
interval halving procedure (bisection) to locate the discontinuity within allowable error
limits. More advanced techniques involve interpolation or regula falsi (false position)
[11,13,14,15,16,17,18].
Summary
A review of the available literature has shown the value of using cubic splines.
Potentially, cubic splines may be used to eliminate infinite slope problems and to smooth
hydraulic model discontinuities. As a result, numerical integration of hydraulic ystem
can be made compatible with ODE solvers designed specifically for mathematically stitT
ODEs. However, none of the available literature provides a simple, physically
significant, method to generate a variety of cubic splines. Similarly, a physically
significant algorithm for handling state events was not contained in the available
literature. The complexity of stiffODE solvers containing adaptive step size and event
locating algorithms has forced designers to use commercially available numerical
integration packages. These packages are generic and are not specifically designed for
18
hydraulic systems. An algorithm for handling typical hydraulic system state e nts
would be extremely valuable to hydraulic system designers using commercial integrators.
In order to increase the value of computer analysis, steps must be taken to
facilitate compatibility between hydraulic systems and stiffODE solvers. Aver atile
cubic spline modeling tool is required to eliminate discontinuities in a variety of
functionspecific component models. In addition, a procedure for handling state events
must be developed to realize the potential gains in computational efficiency offered by
stiffODE solvers.
CHAPTER 3
MATHEMATICAL MODELS AND CUBIC SPLINES
Introduction
Infinite slopes and discontinuities in mathematical models are a source of
problems for stateoftheart numerical integration packages. Numerical ODE solvers
specifically designed for mathematically stiff systems require a Jacobian matrix.
Unfortunately, an infinite slope leads to an illconditioned Jacobian matrix and causes
the sol ver to fail. Abrupt changes to the Jacobian matrix, in the fonn of model
discontinuities, also cause the ODE solver to fail. In addition, adaptive step size
algorithms become extremely inefficient when a discontinuity is encountered and
considerable computation time is wasted.
Cubic splines have been used to model problem causing hydraulic system
components with some success [5,8]. However, existing cubic spline models are
function specific and cannot be used as a general modeling tool. It is necessary to
develop a versatile method that can be used to eliminate a variety of infinite slope
problems and to create a smooth bridge between any two discontinuous functions. In this
way, the dynamic analysis of stiff hydraulic systems can be advanced.
Background
Typically, cubic splines are used as interpolating polynomials for a set of data
points [19]. For n+ f data points, n third order polynomials are generated to interpolate
between the data points as shown in Figure 3.1. At each intenor point, the polynomials
are continuous in position, slope (I st derivative), and curvature (2nd derivative). At the
20
end points ofthe data set, no joining polynomial exists. As a result, the slope and
curvature are not constrained.
y
Point 2
Point 1
x
Figure 3. 1: Cubic Spline Interpolation of Data Points
TwoPoint Cobie Splines
For modeling purposes, it is often necessary to "bridge" two discontinuous
functions. A typical example is shown in Figure 3.2. It is desirable to generate a single
y
X
Figure 3.2: Two Functions Requiring a Smooth Connecting Cubic Spline
21
cubic polynomial with predetennined slope properties at the end points to provide slope
continuity as depicted in Figure 3.3. It is possible to force the end point slopes to assume
y
Cubic Spline
x
Figure 3.3: Cubic Spline Providing Slope Continuity Between Two Functions
any desired value. The twopoint fonn of a cubic spline with specified end point slopes
can be derived from the general case as follows:
A cubic spline connecting two points (XI, Yl) and (X2, Y2) takes the fonn:
Where,
(3.1)
The coefficients G, h, c, and d are expressed in terms of the second derivatives
(curvatures) at the end points [19]:
h=~
2
(3.2)
(3.3)
Where,
c _ Y2 YI  2hSi +h 2 '~
h 6
SI = Second Derivative at (Xl, YI)
S2 = Second Derivative at (Xl, Y2)
h = Interval Size: (XlX,)
(3.4
(3.5)
For cubic spline interpolation on n+ 1 points, the following formulas for forcing end point
slopes apply [19]:
Where,
I S+ 2h ('  1,,1 n 1"f'J,,+1  6(YI'I+I  YII+Ih" Y" )
h}/ = nlSJ Interval Size: (XnI XII)
S\ = Second Derivative at (Xl, yd
SnI = Second Derivative at (XnI,Yn/)
Y; = First Derivative at (x J, YI)
Y;,+J = First Derivative at (Xtl ·1, YII' J)
(3.6)
(3.7)
These equations may be simplified fOT two points as follows:
25 +5 = 6 (Y2YI _Y') 1 2 h h I (3.8)
(3.9)
Let, (3.10)
Substituting gives:
Solving these equations simultaneously produces:
(3.\1 )
(3.12)
(3.13)
(3.14)
Using this result, the cubic spline coefficients from Equations 3.2, 3.3, and 3.4 may be
computed to force the slopes at the end points to the specified values of y; and y; .
In order to link two discontinuous functions with a cubic spline and maintain
slope continuity, the first derivative of each function at the juncture points must be
detennined. Unfortunately, hydraulic systems oft~n contain highly nonlinear
relationships between operating parameters. As a result, obtaining the first derivative
analytically is time consuming and impractical. Numerical calculation of the first
derivative is the most practical alternative. This task can be effectively accomplished
using the standard forward difference approximation. This technique uses a finite
approximation of the infinitesimally small change in the independent variable that
24
defines a derivative. Details regarding the forward difference method are available in
standard Numerical Analysis texts [19].
FourPoint Cubic Splines
In order to eliminate the first derivative requirement ofthe twopoint method it is
necessary to use a natural spline. A natural spline imposes no slope or curvature
requirements at the end points. As a result, the end cubics approach linearity at their
extremes. Unfortunately, a natural spline connecting only two points is a straight line.
Obviously, a straight line is not a useful modeling tool when dealing with the inherent
nonIinearities of hydraulic systems.
To correct this problem while still using natural splines a modeling tool based on
the fourpoint fonn of cubic spline interpolation can be created. The basic concept of
this modeling tool is best illustrated with an example. Consider tne four points shown in
Figure 3.4A. Assume points Band C are fixed and points A and D may be moved.
Moving points A and 0 effects the slopes at points Band C and changes the shape of the
curve between points Band C as shown in Figure 3.4B. Therefore, points A and D may
be used as "handles" to control the slopes at points Band C.
Returning to the original goal, a cubic spline for modeling purposes connects two
points with predetennined slope values at each point. Using the fourpoint fonn, the
modeler can fix the position of the endpoints with two of the four points. The other two
points may be used to control the slope val ues at these endpoints. This technique is
conveniently implemented by locating control points in such close proximity to the true
end points that the difference is physically insignificant. Mathematically, however, the
2
......co. B
(A)
",orI x
y
IB)
Figure 3.4: Effect ofMoving a End Points (A to A' and D to D')
di fTerence is instrumental in detennining the slope at the endpoints. This di fference
needs only to be larger than the computational precision of the computer or program used
to generate the spline.
The fourpoint technique is best illustrated by returning to the example depicted
in Figure 3.2. Again, suppose it is necessary to provide a smooth connection between the
two functions. The first endpoint of the connecting spline is identified as (x /.y /). This
point lies onfi(x). In order to generate a control point for this endpoint, the original
function equation (fi(x» is used. An initial estimate on the order of computational
precision is added to x I as follows:
where, X Ie = Xcoordinate of the control point
(3. 15)
26
XI = Xcoordinate ofthe original end point
L1x = Increase in x on the order of computational pr cision
The functionIi is then evaluated at X Ie and compared to fi(x I) as shown below:
where, Yldif/ = Difference inYI as calculated by computer program
(3.16)
[fthe computer cannot discernYldif/fi"om zero, Ax is gradually increased until a
calculable difference exists. This iterative process is easily accomplished in a computing
environment. After this process is complete, the original end point (Xf,YI) and the control
point (XlcoYIJ wi]] be discernibly different in x and y but the difference will be physicaJly
negligible. Th.is procedure must be repeated using/ix) to obtain a control point for the
other end point. A graphical depiction ofthe endpoints and their control points is shown
in Figure 3.5. The relative distances have been exaggerated for visibility. (NOTE: If
jj(x) or/ix) is a horizontal line, the corresponding iteration procedure ony is not
necessary because Y does not vary with x.)
y
• Control Point
• Original End Point
x
Figure 3.5: Cubic Spline End Points and their Control Points
27
The original end points, along with their control points. may now be used to
generate three cubic spline polynomials using the procedure outlined in Appendix A.
The X and Y matrices take the form:
XI
X=
x lc
X
X 2c
(3.17)
Because the true end points are so close to their corresponding control points, the cubics
connecting them may be ignored. As a result, !he cubic bridging the middle interval is
assumed to connect the original endpoints. This assumption is true for all practical
purposes. Using the functionsjj(x) andh(x) to create the control points forces the slopes
at the end points to match the original functions. The resulting cubic spline model is
depicted in Figure 3.6
y
Cubic Spline
x
Figure 3.6: Cubic Spline Generated with the FourPoint Fonn
A complete algorithm for generating control points is displayed as a flow cha,rt in
Figure 3.7.
Identify the FL¥1C1ion
!(Xl on which the
EM Foi1t lies
2
DeltaX =E.
where E is
Computational
Precision
Xc = X .... DeltoX
Control
Point is:
(Xc, y)
Xc = X + DeltoX
Compute fIX cl
Cannol
Point is:
(Xc,Yc)
, Delto)(;
'_DelloX }., x DeltoX
v.nereklsa1
wement foetor
Figure 3.7: Flow Chart for Generating Control Points
Comparison of TwoPoint and, ,FourPoint Splines
The twopoint and fourpoint methods for generating cubic splines rna both be
used as modellng tools. Each ofthese methods has advantages over the other in tenns of
computation, The twopoint method requires knowledge of end point first derivatives
while the fourpoint method does not require derivatives. Unlike the twopoint method,
however, the fourpoint method requires the solution of a set of equations.
A preferred method may be detennined by evaluating the resulting cubic spline .
The twopoint and fourpoint methods may not produce the same cubic function. The
difference between the methods is best illustrated with an example. The following
straight line equations are given:
j(x,) = 0
j(X2) = 0.9x 720
(3.18)
(3.19)
It is desired to generate a cubic spllne between the points (720,0) onj(xf) and (825, 22.5)
onj(x}). Applying both methods to this problem produces the results shown in Figure
3.8. Slope continuity at the juncture points can be maintained adequately by using either
method. However, the cubic spline created using the twopoint method "dips" below
zero before heading in a positive direction. The potential for this behavior exists when
"tight turns" are involved. This characteristic is undesirable when dealing with flow rate
on the Yaxis. The negative sign reverses the flow direction. Flow reversal does not
exist in practical hydraulic systems. The fourpoint spline creates a more direct bridge
between the two discontinuous functions. As a result, the fourpoint method is preferred
o
as a modeling tool when the relatively sharp comers of a typical hydraulic model are
involved.
/
/
/
/
FourPoint
A
Spline "... //
'&..,V/
.,/ V
~ ~ ~TwoPoint 
Spline 
50
45
40
35
30
25
~O
15
10
5a
5
700 725 750 775
X
800 825 850
Figure 3.8: Comparison of TwoPoint and FourPoint Cubic Splines
Application Examples
The usefulness of cubic splines as a modeling tool will be illustrated by two
examples from the world of hydraulics. In the first example, an empirically based orifice
model is developed. The second example involves smoothing discontinuities in a relief
valve model.
Example 3.1  An Empirically Based Orifice Flow Model
PART A
An orifice with a diameter (D) ofa.0] 2 inch is known to have a discharge
coefficient (Cd )ofO.63. Laboratory testing was perfonned to determine the
pressureltlow relationship of this orifice using SAE lOW oil at 73°F (See Appendix B).
The measured data have been plotted in Figure 3.9. It is desired to develop a
mathematical model based on this experimental data.
"'J
v ~
.'
/
,...
V V
V V
V
V
0.12
0.14
0.02
o
o 10 20 30 40 50 60 70 80 90 100 110 120
Pressure (paid)
u
..; 0.1
c
; D.DB
~
.; 0.06
a::
~ 0.04 u:
Figure 3.9: Experimental Results of Orifice Test
In order to develop a representative model, it is first necessary to determine where
the turbulent flow equation is valid. This task may be accomplished by plotting the
experimental data along with the turbul.ent flow equation. The turbulent flow equation is
[20]:
(3.20)
Q= Flow Through the Orifice (Dependent Variable) (in3/sec)
Cd = Discharge Coefficient = 0.63
A = Flow Area of the Orifice = D
2
Jr (in2
)
4
p = Fluid Density = 8.171x 105 (lbfsec2Iin4
) for lOW oil at 73°F
t1P = Pressure Drop Across the Orifice (Independent Variable) (psid)
The results are depicted in Figure 3.10. From Figure 3.10, the turbulent flow equation
I I ~
,....
Twbudent Flow ........
Equall00 ~
~~V ~
~V
./
:,/"~ ....
~
,., Test Results
/.
~V
V o
o 10 20 30 40 50 60 70 80 90 100 110 120
Pressure (paid)
0.14
0.02
u:
0.1
~
"; 0.08
~
~ 0.06
u:
~ 0.04
iL
0.12
Figure 3.10: Comparison of Test Results and Turbulent Flow Equation
and the test data match almost exactly at 85 psid and above. Using the turbulent flow
equation at low pressure differences would obviously be inaccurate. More importantly,
the infinite slope at the origin would cause a stiffODE solver to fail. Therefore, the
turbulent flow equation will be used to model the orifice flow at pressures above 85psid.
The turbulent flow model is not accurate at pressures lower than 85 psid because laminar
and transitional flow are dominant. In order to model the portion of the flow/pressure
relationship, a cubic spline must be used. It is first necessary to detennine the end points
ofthe cubic spline. From Figure 3.10, the first end point is obviously at the origin. The
second end point can be detennined by calculating the flow rate at 85 psid with the
turbulent flow equation (3.20).
(3.21 )
Q= Flow Through the Orifice (in)/sec)
,1P = 85 psid
With the two end points in place, the control points may be generated. To obtain
a smooth transition to turbulent flow at 85 psid, the control point is detennined using the
turbulent flow equation (3.20) and the algorithm outlined in Figure 3.7. There are no
predetennined slope requirements at the origin because no joining function exists.
Physically, laminar flow is dominant near the origin. The equation for laminar flow is
linear [20]. Therefore, it is reasonable to choose a linear relationship to create the
control point at the origin. The slope of this line is detennined using the original test
data. From Appendix B, the data point (,1P, Q) nearest the origin is:
(2.5 psid, 0.00693 in. 3/sec)
The equation of a line passing through this point and the origin is simply:
Q= 0.00693 .1P
2.5
Simplifying gives:
Q = O.00277(LtP)
(3.22)
(3.23)
The algorithm described in Figure 3.7 along with equation 3.23 may now be used to
generate a control point at the origin.
Using a computer to calculate the end point and control point values gives:
End Point 1: (0.0)
Control Point 1: (2.220446049250313 x 10'13.6.150635556423367 lO'16)
End Point 2: (8.500000000000000 x 101, 1.027731761455279 x 10")
Control Point 2: (8.500000000000222 x IO J
, 1.027731761455292 x 10,1)
The number of significant figures has been carried to an extreme to show the difference
between the original end points and their control points. Obviously, the differences are
physically insignificant. However, these differences are critical when generating the
cubic splin~ because they force the slopes at the end points to the desired values.
Using these fOUT points and the procedure outlined in Appendix A, the following
cubic spline was created:
Q = 1.327 X 10,7 9 3 2.965 X 10'5 ~ + 2.770 X 10'3 t1P (3.24)
The procedure actually produces three third order polynomials. However, the cubics for
the intervals between the end points and their control points may be ignored because the
points are so close together. The resulting mathematical model is displayed in Figure
3.11.
~
.............
~
~
~
Cu ic Sp me....... L,.....:::;
~~~T, Q~o ,It.
~,,;'
V
/ o
o 10 20 30 40 50 60 70 80 90 100 110 120
Pressure (psid)
Figure 3.11: Cubic Spline Based Orifice Model
0.14
0.02
u..=..., 0.1
.5
;; 0.08
~
!J 0.06
a::
~ 0.04
Li:
0.12
3
The model may be refined by experimenting with the departing slope at the
origin. If the slope in equation 3.23 is changed from 0.00277 to 0.00220, the model is
improved as shown in Figure 3.12.
v V .......
V ~
~
~~ ~
Test esult:/V
i'oo."t,.
~ . ',".
/ o
a 10 20 30 40 50 60 70 80 90 100 110 120
Pressure (paid)
0.14
0.02
u
~ 0.1
.5
:i 0.08
~
! 0.06
a::
~ 0.04
u:
0.12
Figure 3. 12: Refined Orifice Model
PARTS
Using the results from Part A, it is desired to create a generalized mathematical
model for flow through any orifice. This task may be accomplished by relating the
results of Part A to Reynolds Number (Nr ). Flow characteristics at a particular Reynolds
Number are consistent for the flow of any fluid through any size orifice. Therefore,
Reynolds Number provides a tool to generalize the orifice model developed in Part A.
The function used to force the departing slope at the origin was based on the
linearity oflaminar flow. The laminar flow equation is [20]:
(3.25)
3
Q = Flow Through the Orifice
8 = Laminar Flow Coefficient
D = Hydraulic Orifice Diameter
D2
J[
A = Flow Area ofthe Orifice = 4
p= Fluid Density
v= Kinematic Fluid Viscosity
AP = Pressure Drop Across the Orifice
For Reynolds Numbers less than 100, the following relationship is generally accepted
[20]:
Rearranging gives:
Substituting into Equation 3.25 produces:
Q= 2Cd
2
DA AP
pvNR
The slope ofthis line is simply:
2C/DA
slope =
pvNR
(3.26)
(3.27)
(3.28)
(3.29)
From Part A, a slope of 0.00220 was used to generate the final mathematical model. [f
this slope and the rest of the constants from Part A are substituted into Equation 3.29, it
is possible to solve for Reynolds Number. The only new parameter is kinematic viscosity
37
(v). For the lOW oil used in Part A, kinematic viscosity is 0.1623 (,in. 2/sec). Performing
this calcuJation resuJts in a Reynolds Number of 36.92. In the general ca e, th refor the
following equation may be used to create a control point at the origin:
0= 2C/DA L1P
 pv(36.92)
(3.30)
The position of the second end point must also be related to Reynolds Number.
Reynolds Number is defined by the following equation:
(3.31)
From Part A, the second end point was located at:
P = 85 psi and Q = .10277 inJ/sec
Using this value for Qand the known values for A, Dand v, the Reynolds Number at the
end point is easily calculated as 67.19. With this Reynolds Number known, the flow rat
at the second end point for any orifice may be calculated by rearranging Equation 3.15 as
follows:
Q= 67.19Av
D
(3.32)
The corresponding pressure is calculated by rearranging the turbulent flow equation
(3.20):
(3.33)
38
Equations 3.32 and 3.33 represent a generalized method to calculate the location ofth
second end point for any orifice/fluid combination. The turbulent flow equation (3.20)
may now be used to create the control point as outlined in Part A.
To further generalize, it is desirable to allow the user to select the Reynold
Number at which the turbulent flow equation applies (NRr). This feature is made possible
by using a ratio of the two Reynolds Numbers calculated above. These values are:
Reynolds Number Used in Laminar Flow Equation (NRI) = 36.92
Reynolds Number Used in Turbulent Flow Equation (NRr) = 67.19
Rau. o= 36.92=0.55
67.19
(3.34 )
Using this infonnation, Equations 3.30 and 3.32, may be further generalized as follows:
Q= 2C/DA t1P
pv(0.55NRr)
NR,Av Q='"'
D
(3.35)
(3.36)
The value ofNRr may now be defined by the modeler without effecting the basic shape of
the cubic spline model.
Example 3. I was developed in the MATLAB computing environment. The
MATLAB script files used to complete this example are contained in Appendix C.
Example 3.2  Static Relief Valve Model
A common direct acting pressure relief valve is depicted in Figure 3.13. When
the force due to pressure at port P becomes large enough to overcome the spring preload
against the poppet, fluid begins to flow to port T. This pressure is known as the cracking
pressure. As the combined forces due to pressure and flow increase the spring continu
to compress until the mechanical stop is reached.
Port T
Port P
Figure 3.13: Typical Pressure Relief Valve
In many hydraulic circuits, it is acceptable to ignore the dynamics of the poppet.
The simplified mathematical model is called a sialic reliefvalve model. A typical static
relief valve model is presented in Figure 3.14. The flow through the valve remains at
zero until the cracking pressure (Per) is attained. At this point, a linear relationship
o
between pressure and flow is assumed as the pressure at port P works against th prt ng.
When the upper mechanical stop is reached (PmlCf) the poppet is static valve begins and
the turbulent flow equation applies.
Figure 3.14: General Static Relief Valve Model
o 100 200 300 400 500 600 700 800 900 1000 1100
Pressure (psi)
60
50
E 40
CL
.5!
~ 30 CG a:::
~
.2 20
u.
10
o
PmlLX ....
',.......H
Strnig l Tur ulcnl
Line ~..... Flo ·Eq.
Per ...... r....
The static relief valve model contains two discontinuities which could severely
Iimit the computational speed of a numerical ODE solver. The goal of this example i to
smooth these discontinuities using cubic splines. For this example, the following
infonnation is known:
Per 800 psi
Pmax = 850 psi
Qma'C = 45.05 in3/sec
Q@Per = 0 inJ/sec
P at Port T = 0' psi
41
The first step is to develop working mathematical relationships. From zero
pressure to Per the relationship is obviously:
Q=O for 0 ~ P ~ Per (3.37)
The relationship from Per to Pmax is simply a straight line equation:
Q= mP + b for Per ~ P ~ Pmax (3.38)
where, 45.05 0 90' 'I . 850800 =. In' seepsI (3.39)
b= Qmax  mPmax = 45.050.90(850) = 719.95 inJ/sec (3.40)
For pressures above Pmax the turbulent flow equation is utilized. The turbulent flow
equation (3.20) may be simplified by lumping the constants together as shown in
Equation 3.41. In this case, i1P equals P because port T is at zero pressure.
for P ~ Pmax (3.41 )
Given Pma"( and Qman kv is easily calculated:
(3.42)
With the mathematical relationships in place, it is now possible to generate cubic
splines using the fourpoint method. The discontinuity at Per will be considered tirst.
Reasonable end points must be selected to bridge the discontinuity. The discontinuity
occurs at Per (800 psi). End point pressure values of 720 psi and 825 psi will be selected
hecause this range spans the discontinuity with enough room to create a relatively smooth
curve between the functions. This choice is only one of many possibilities. Using a
computer algorithm allows the modeler to test any number of combinations. The flow
4
rates corresponding to the chosen endpoint pressures may be calculated using Equation
3.34 and 3.35. A computer algorithm based on Figure 3.7 wa used to compuli the end
point values and the control point values:
End Point 1: (7.200000000000000 x ]02,0)
Control Point I: (7.200000000000002 x 102,0)
End Point 2: (8.250000000000000 x 102, 2.252250000000004 x 10 1
)
Control Point 2: (8.500000000000002 x ]01,2.252250000000015 X ]01)
An inspection of End Point 1 and its control point reveals that no migration is necessary
for flow rate (dependent variable) because the original function (Equation 3.37) has a
slope of zero.
These points may now be used to generate the cubic spline according to the
procedure outlined in Appendix A. Again, the resulting cubics between the end point
and their control points may be ignored because the difference between these points is
physically negligible. The resulting cubic spline equation for the middle interval is:
Q = 6.440 X 106 p3 + 1.367 X 103 p2 + 1.943 X 1016 P (3.43 )
To complete the model, the entire procedure must be repeated using Equations 3.38 and
3.41. Using endpoints at P = 840 psi and P= 900 psi, the following cubic spline equation
is produced:
Q= 5.206 X 105 p3
 8.592 X 103 p2 + 5.000 X 101P + 3.604 (3.44)
The resulting static relief valve model is presented below (eqs. 3.453.49) and d picted
graphically in Figure 3.] 3. Although some small amount of accuracy i lost, th re ulting
model contains no discontinuities. This feature provides compatibility with sti fTODE
solvers and adaptive step size algorithms.
Q=O for 0 psi s P s 720 psi (3.45)
Q= 6.440 X 1O.Q (p_720)3 + 1.367 X 103 (p720i + 1.943 x 1016 (P720) (3.46)
Q= 5.206 X 105 (P840)3  8.592 X 103 (P840)2 + 5.000 X 10'\ (P840) + 3.604 (3.48)
for 720 psi < P s 825 psi
Q = 0.90P  719.95
for 825 psi < P s 840 psi
for 840 psi < P s 900 psi
(3.47)
Q = 1.545 JP (3.49)
for 900 psi < P
uhlc plmc'"r1,
I I UJ "U'''"I
FI w Eg.
• lTalQht
line .. ,
I
( ubic
.plil1~ ~,
50
60
10
o
o 100 200 300 400 500 600 700 800 900 1000 1100
Pressure (psi)
$ 30 l'lJ
0::
~
.2 20 u..
E 40
Q.
Cl
Figure 3.15: Smoothed Static ReJiefYalve Model

44
As with Example 3.1, Example 3.2 was developed in the MATLAB computing
environment. The MATLAB script files used to complete this example are contained in
Appendix D.

45
CHAPTER 4
THE EVENT SWlTCHING ALGORITHM
Background
The computational efficiency of a stiff numerical integration routine may be
severely impeded by the presence of state events. A state event is a discontinuity which
occurs when a variable reaches a physical limit during numerical integration. State
events signal a change to the original set of ODEs. A new integration problem is defined
at the exact time the critical value is achieved. The sudden switch to a new set ofODEs
produces an abrupt change to the Jacobian matrix of a stiff solver. This abrupt change
often causes the ODE solver to fail. In addition, adaptive step size algorithms tend to
"hunt" around these state events in an effort to traverse them by reducing the step size.
This hunting process requires an inordinate amount of processing time.
To combat these problems, researchers have recommended restarting the ODE
solver each time a state event is encountered (12,13]. This process divides the original
problem into continuous sections and solves them in a piecewise manner. Several event
location routines have been devised to pinpoint the precise time at which a state event
occurs (11,12,] 3,]4,15,16,17,18]. The goal of these routines is to control step size such
that the state event occurs at the end ofthe step. Commercially available integration
packages locate events by monitoring sign changes. Therefore, each state event must be
defined by a discontinuity[unction. The discontinuity function is a mathematical entity
involving the state variables of the original ODEs. The state event occurs when the
discontinuity function equals zero.
4
Commercial numerical integration packages are intended to be generic. As such
it is left to the user to define the discontinuity functions for a given application.
Hydraulic engineers would benefit greatly from a functionspecific "book keeping"
approach to state event handling. This type of approach may be developed b_
investigating the physical nature ofhydraulic system state events.
Hydraulic System State Events
fn the field of hydraulics, state events typically occur when actuators loads or
internal valve parts encounter mechanical travel limits. A frictionless springmass ·~· damper system may be used to investigate the effects of these mechanical stops (See
Figure 4.1). The motion ofthe mass (m) is constrained by the two mechanical stops. ff a
Mass
(m)
Dclrnp;lg C08!tIclent (bl
X=Xmax
X=Q
Figure 4. ] : SpringMassDamper System with Mechanical Stops
slowly increasing force (F) is applied to the mass, three distinct states are revealed.
These states are depicted in Figure 4.2. From Figure 4.2, State A shows a static condition
because the applied force (F) is not large enough to overcome the spring preload (f~reload)

State A:
Applied Force (F]
Mass .......of (m)
State B:
Applied Force (F)
Mass
(m)
State C:
Applied Force WI
Mass
(m]
Figure 4.2: Three States for SpringMassDamper System with Mechanical Stops
47
,
"·~c"
:'
•". '
and no motion occurs. The same is true of State C. In this case, the mass is pinned
against the left stop because the force due to spring compression (k·x) is limited by this
mechanical stop. By definition, therefore, the velocity ( i ) and acceleration (j )of the
mass are zero for States A and C.
If the magnitude ofthe applied force is between the spring preload and the
maximum spring force and the mass is not against a stop, the mass is either dynamic or
potentially dynamic. This condition is depicted as State B in Figure 4.2. The mass
would be potentially dynamic if the applied force (F) was constant. Any change in the
applied force would result in motion. The equations of motion for State B may be
48
developed by summing the forces acting on the mass.
mi = F  Fprdoad kx· bi
Solving for acceleration gives:
,t = .! (F  FpreloQlr Iexbi )
m
(4.] )
(4.2)
··•
~,·•
· '
For numerical integration, a set of first order differential equations is required. This set
ofequations is obtained by introducing XI and X2 as follows:
(4.3)
(4.4)
(4.5)
Therefore, the equations of motion for State Bare:
(4.6)
(4.7)
49
For States A and C. the velocity and acceleration are zero. The equations of motion for
these states are:
(4.8)
(4.9)
During the course of numerical integration, therefore, it is necessary to switch between
these two sets of differential equations. The appropriate set ofequations at any given
time is dictated by the magnitude of the applied force. Typically, the magnitude of the
applied force depends on timevarying hydraulic pressures. As a result, it is necessary to
track the magnitude of the applied force during integration and to switch equation sets
when a threshold value is achieved. This process assumes instant deceleration when a
mechanical stop is encountered. In reality, deceleration is not instantaneous but it is so
nearly instantaneous that the assumption is valid.
Switching between two sets of differential equations produces an abrupt change
to the Jacobian matrix of a stiffODE solver. In addition, adaptive step size algorithms
require excessive computational effort to traverse a switch discontinuity. As a result, it is
necessary to restart the solver when a state event is detected [12,13]. Before the solver is
restarted, however, it is necessary to integrate up to the precise time of the state event.
Commercial numerical integration packages often locate events by monitoring sigll
changes. In order to implement a commercial event location scheme, a discontinuity
function (t/J) must be created for each state event such that the discontinuity occurs when
'I
..
'I
'. J
o
the function equals zero. The discontinu"ty functions for the springmas damper yst m
must take the fonn:
More specifically, the discontinuities for each state may be written as follows:
(4.10)
State A:
State B:
State C:
¢a = X max  X
or
¢s= x
¢c = F  FpreloaJ k(xmax)
(4.11 )
(4.12)
(4.13)
(4.14) .
:'
Zero initial conditions will be assumed to demonstrate the purpose of these discontinuity
functions. Under these conditions, the set of equations for State A are applicable (4.8 and
4.9). It is also assumed that the applied force (F) increases with time. At each
successive time step, the discontinuity function tPA (4.11) is evaluated. If at any time ¢A
becomes negative, the event location routine is activated to pinpoint the time at which ¢A
equals zero. Physically, ¢A is zero when the applied force (F) is large enough to balance
the spring preload. The numerical integrator is stopped at this point. Using the state
variable values at the point oftennination as initial conditions, the ODE solver is
restarted using the State B differential equations (4.6 and 4.7). Coinciding with the
switch to a new set of ODEs, the discontinuity function must be changed. Equations 4.12
and 4.13 are appropriate for State B. One ofthese discontinuity functions will cross zero
if the displacement of the mass fans outside the limits imposed by the mechanical stops.
ff a sign change occurs, the event locator will again pinpoint the time at which the state
I
'1
change occurs and the solver may be restarted with the static equations of motion (4.8
and 4.9). In the case of an increasing applied force, the left stop will be encountered. A
such, the discontinuity function r/Jc (4.14) must be invoked. If the applied force later
begins to decrease, ¢c will define the point at which the applied force is no longer large
enough to keep the mass pinned against the stop. The event location and solver restarting
procedure may then be repeated as the mass enters State B. Similarly, the entire
procedure must be repeated if the mass were to contact the right stop and enter State A.
Development of the Event Switching Algorithm
This ongoing event location and integrator restarting process may be efficiently
handled by introducing an F;vent Switching Algorithm (ESA). The ESA takes advantage
of physical properties to simplify the "book keeping" required to numerically integrate
systems containing state discontinuities. The ESA is developed by comparing the
discontinuity functions with the second dynamic equation of motion (eq.4.7). For State
A these equations were:
::
¢A = Fpreload  F (4.15)
(4.16)
The mass is pressed against the right mechanical stop while in State A. Therefore, the
displacement (XI) of the mass is zero. The velocity (X2) of the mass is also zero.
Inserting this information into Equation 4. l6 gives:
(4.17)
If equation 4.17 is set equal to zero, the following is true:

1
0=  (F  Fpre1oot/)
m
o= F  Fpreload
.18)
(4.19)
The state event occurs when the discontjnuity function (tPA) from equation 4.15 equals
zero:
o= Fpreload  F
Comparing equations 4.19 and 4.20 reveals the following correlation between the
discontinuity function (tPA) dynamic equation of motion ( ;(2 ):
(4.20)
(4.21 )
Consequently, the existing equation of motion (4.16) may be used as the discontinuity
function in State A. This fortunate circumstance eliminates the need for a separate
discontinuity function. rn addition, the negative sign or "1" may be used as a flag to
identitY State A within the ESA. This negative sib'll is important because it dictates the
direction of zero crossing to the event location routine (i.e. positive to negative).
A similar analysis may be performed for State C using rPc and x2 . These
equations are reprinted below for convenience:
t/Jc = F  f~rdoad  k (x max ) (4.22)
(4.23)
In this case, Xl equals Xmax because the mass is pinned against the left mechanical stop.
The velocity(x2} is again equal to zero. Substituting this information into equation 4.23
and setting both equations (4.22 and 4.23) equal to zero reveals:

3
(4.24
Therefore, State C also does not requjre a separate discontinuity function (tA:). In
keeping with State A, a flag of U+ 1"may be used by the ESA to identify State A.
The discontinuity functions for State B (eqs. 4.12 and 4.13) are related simply to
displacement (Xl)' This information is readily available during the course of numerical
integration. A zero ("0") flag may be used to identify State B within the ESA.
With the discontinuity functions in place, an event location routine may be used
to locate the precise time ofthe state event. The ODE solver will almost always "step
over" the state event. At the time step just prior to the state event, the discontinuity
function is positive. The discontinuity function then becomes negative after the next
time step. Because the state event occurs when the discontinuity function is zero, the
actual time of the event is between the two steps. The goal of an event location
algorithm is to iteratively shorten the step size until the discontinuity function is zero at
the end of the current time step. This task is most simply performed by uccessively
bisecting the time step until the discontinuity function is within some tolerance ofzero.
More efficient methods often involve interpolation or false position schemes.
The event location routine effectively integrates to the precise moment ofthe
state event. The time at which the state event occurs, as well as all of the corresponding
state values, are known. At this point, the ODE solver is stopped. The necessary changes
are made to the system of differential equations and integration is restarted using the time
and state values obtained by the event location algorithm as initial conditions. The only
54
change to the initial conditions involves the velocity ofthe mass. This vel.ocity must be
reset to zero when the mass is pinned against either of the mechanical stops.
The resulting Event Switching Algorithm is best represented in flow chart form.
This flow chart is displayed in Figure 4.3. The ESA contains a simple method to
seamlessly integrate numerically stiff systems containing state event discontinuities.
Although useful in many engineering disciplines, the ESA is especially valuable to
hydraulic system engineers because stiffness and state discontinuities are commonplace.
5
Reset
Velocity
x, = a
set Initial
State Flog
\Ulue
set Initial No
Conditions
to Values
of Last
Step
set State Use New VOkJes
to EvokJate
Flag =1 2nd EQJOliOn
of Mallon fOl
State 8: Reset
)(,  l(X .. X,. F) VelocNy
Save x, .. 0
set State Results Flag =+ 1 ~Result
By Stcte Flog
\otlue to get
DlscontnJlly
Ft.r)c11on 0B
Use Event
Save Lcx:oIion
Results SCheme to Set Inltiol
Integrate 10 Condttfons
Time'M1ere 0 A=X,New =0 to Values
IvoftNn 1OIetOnCe~ Use Event of Lost
Location Step
SCheme to
Use Event
Integrate to
Time 'M'lere
Location
0B =0 SChema to Exit ODE
Integrate to Solver ,"""*' Iololoncol
Exit ODE Time Where
Solver 0c=
XMo>  X'New Extt ODE save =0
(WlthIn lolet\Ylce) Solver Results
Figure 4.3: Flow Chart for Event Switching Algorithm as Applied to General
SpringMassDamper System
Application Example
The application ofthe Event Switching Algorithm will be demonstrated with an
example. This example involves a hydraulic actuator encountering a travel limit. In this
case, an arbitrary stroke limit of 2 inches is placed on the cylinder. Only cylinder
extension is considered in this example.
Example 4. J  Hydraulic Cylinder Circuit
Consider the hydrau.lic circuit shown in Figure 4.4. This control circu.it is
designed to impart translational motion to a load using a hydraulic cylinder. Critical
pressure nodes have been identified on the schematic. The mathematical
Hydraulic
Cylinder rTE=~::JMass 1ooIIW\rL (m)
Meterlng\ < Orifice /
L.:.:...:..J 0 psi
Reservoir
Figure 4.4: Hydraulic Circuit Schematic for Example 4.1

models for the circuit components and the resulting set ofdifferential equations are
developed as follows:
Mathematical Models
LOAD & HYDRAULIC ACTUATOR:
A springmassdamper system will be used to model the load. The load mass is rigidly
fixed to the piston rod ofthe hydrauhc cylinder. As a result, the mass of the rod and of
the load must be combined. This combination is called a lumped load model. In this
case, the following information is known:
Lumped Load Weight (w): 3860lb
Spring Constant (k): 700 Iblin
Damping Coefficient (b): 200Ib/(inlsec)
The force applied to the load is provided by hydrauhc pressure acting against the piston
area. Hydraulic pressure is present on both sides of the piston. The resulting applied
force is:
7
(4.25)
The piston dimensions are as follows:
Bore Diameter:
Rod Diameter:
Therefore,
4.0 in.
2.5 in.
Abore = (4.0)2.]t /4
Arod = Abore  (2.5f1t /4
(4.26)
(4.27)
58
DIRECTIONAL CONTROL VALVE:.
When shifted, the directional control valve directs pressurized fluid to either the rod or
cap end ofthe hydraulic cylinder. Depending on the direction the valve is shifted, the
cylinder will either extend or retract. In this case, the valve will be shifted to force the
cylinder to extend at time zero (t = 0). The flow restriction through the open valve rna
be modeled as an orifice. For this directional control valve, the effective orifice
diameter is:
Ddc = 0.27842 in.
Given the following, an orifice model may be created as outlined in Example 3.1 Part B:
Making the flow area,
(4.28)
Orifice Diameter (DoT) = 0.27842 in.
Discharge Coefficient (Cd) = 0.61
Fluid Density (p) = 7.95 X to5 (Ibosec2lin4
)
Kinematic Viscosity (v) = 0.02 in2/sec
Transition Reynolds Number (NRt ) = 1500
The resulting model is:
(4.30)
This model will apply to flow through the valve in either direction. The return £10\
equation is simply:
For 0:::; (P4) :::; 1.240 psid,
Qdc2 = 9.350 X 10 1 (P4)3  4.649 x 10° (P4)2 + 9.617 X 102 (P4)
For CP4) >1.240 psid,
(4.31 )
(4.32)
METERING ORIFICE:
The metering orifice creates back pressure against the rod end of the cylinder. This back
pressure is responsible for controlled cylinder extension. The following infonnation is
known:
Orifice Diameter (Dur) = 0.11277 in.
Discharge Coefficient ( d) = 0.61
Fluid Density (P) = 7.95 x ] 05 (Ibosec2/inl)
Kinematic Viscosity (v) = 0.02 in2/sec
Transition Reynolds Number (NRr ) = 1500
Given the above data, it is again possible to develop an orifice model using the general
procedure outlined in Example 3. t Part B. The resulting model is:
o
For (P3P4) >7.560 psid,
Q~ ~ C.,4~ ~~(p,  p.,) (4.34)
RELIEF VAtVE:
The relief valve simply limits pump discharge pressure. A model for this relief valve was
developed in Example 3.2. This model is as follows:
(4.37)
(4.35)
(4.36)
for 0 psi ~ PI ~ 720 psi
Qrv = 6.440 X 1O~1 (P I720)3 + 1.367 X 103 (P\720)2 + 1.943 X 1016 (P I720)
for 720 psi < PI ~ 825 psi
Qr.' = 0.90 P1  719.95
for 825 psi < PI ~ 840 psi
Qr. = 5.206 X 105 (P\840)3  8.592 X 103 (P,840i + 5.000 x 10\ (P,840) + 3.604
for 840 psi < PI ~ 900 psi (4.38)
Qrv = 1.545 JP: (4.39)
for 900 psi < PI
PUMP AND MOTOR:
The pump supplies the system with pressurized fluid. An ideal pump model will be used.
Ideal pumps experience no internal leakage. Therefore, the pump flow may be calculated
using the following equation:
Qin = Displacement x Rotational Speed (4.40 )
61
Given a motor speed of 1800 rpm and a pump displacement of 1.28 in3/re ,the pump
flow is:
Qm = 1.28 in
3
x 1800 rev x 1min = 38.5 in
J
rev min 60 sec sec
(4.41 )
This flow rate will be applied to the system as a step input because the direction control
valve is shifted at time zero (t=O). Prior to time zero, the flow is simply bypassed to tank
with a negligible pressure drop through the directional control valve.
Differential Equations
The governing differential equation for pressure is:
p= ~(LQ) (4.42)
p= Pressure at a Given Point ~
13= Bulk Modulus of Elasticity
v= Fluid Volume
Q= Flow Rate
Bulk modulus is a fluid property. In this case, a fluid with a bulk modulus of 150,000 psi
will he used. By convention, a positive flow rate enters a pressure point while a negative
flow exits. The differential equations for each pressure point may be developed as
follows:
NODE I:
The pump discharge flow (Qit,) enters Node I while the flow through the relief valve and
the directional control valve exit. The volume of fluid trapped between the pump outlet
po
and valve inlet is known to be 10 in3
. As a result, the differential equation for this point
15:
NODE 2:
p. _ 150000(0 _Q _Q )
I  10 _ill de rv (4.43)
Node 2 has one incoming flow rate and one outgoing flow rate. The incoming flow
originates from the directional control valve (Qdc)' The outgoing flow fills the extending
cylinder. The volume of fluid in the extending cylinder at any point in time is simply:
The conduit connecting the directional control valve and the cylinder is known to have a
v = Abore oX
Therefore, the volumetric flow rate into the cylinder is
Qcyl in = Abort' 0 i
(4.44)
(4.45)
..•)
3
..
volume of 10 in3
. Given this information, the following differential equation may be
developed:
(4.46)
NODE 3:
The differential equation for Node 3 may be developed in the same manner as Node 2.
Again, a conduit volume of] 0 in) will be used. The resulting equation is:
jJ = 150000 (A 0 i _ Q )
.\ 10 +A. rod nr rod X
(4.47)
63
NODE 4:
The flow from the metering orifice enters Node 4 and the flow through the directional
control valve exits. Given a transmission line volume of lO in3
, the differential equation
for Node 4 is simply:
p = 150000(0 _ Q )
4 10 _or de (4.48)
The standard differential equation of motion for the springmassdamper system
(eq.4.2) is appropriate for the load:
.t = ~ (F  Fpreload  k.xbi )
m
(4.49)
The applied force (F) is provided by the cylinder. This force was modeled as (eq. 4.25):
(4.50)
In this case, the preload force (Fpreload) is zero. As a result, the differential equation of
motion is:
.x.;  ] (l):2Ahore  })yAt rod  k.x. hX. )
m
System ofDifferential Equations
(4.51 )
In order to create a vectorized set ofdifferential equations for the computer
algorithm, the following substitutions were made:
.\'5 = x
64
The resulting set of linked ODEs is:
Results
x = 150000 (A .x Q )
3 ]0 A. rod Ii flr + rod X j
•\ = ]50000 (0 _Q ,)
4 10 _or de ..
xIi = ~ (P._,AboreP.1Arod  k.x5bx(, ) m
This system of equations was solved using Gear's method for stiff ODEs in
(4.52)
(4.53)
(4.54)
(4.55)
(4.56 )
(4.57)
:,
'.,
~ .
.)
..l..,
•,
•
it
conjunction with the ESA as outlined in Figure 4.3. Adaptive step control was also
utilized. The cylinder stroke limit was assumed to be 2.0 inches. Timebased plots of
cylinder rod displacement and velocity are displayed in Figures 4.5 and 4.6 respectively.
The displacement becomes constant when the stroke limit is reached. In addition, the
velocity becomes zero when the stroke limit is encountered. In the absence of the ESA,
these sharp changes would cause the stiff numerical integration and adaptive step size
algorithms to become inefficient or to fail completely. The ESA, however, produces the
...
65
expected results in a reasonable amount of time. In fact, these results were obtained in
under 12 seconds using Pentium 200 computer with 32 megabytes of RAM. A detailed
study of computation time is left for Chapter 5 ofthis document. However, 12 second
is, by no means, unreasonable for a system of this complexity.
Figure 4.5: Cylinder Rod Displacement Response fOT Example 4.1
•
V V
./
,/
V V
/
/
V
_2.5
1Il
QI .r::.
u
:§. 2 c:
CP
E~
1.5
III
0.,
is 1
"8 a:
~ 0.5
~
.5:
>.
(.) 0
o 0.1 0.2 0.3 0.4 0.5
Time (seconds)
0.6 0,7
,•
.t
5...,,..r,.,,
U
QI
.!! 4 #!++'f4fl i
.~ 3 IH+~I~~~ ===+===F===t===:t====t:==rl
o IV
~
"8 2 +i+t1r+++;
a.:.:.
QI
~
.5: 1 +~++1f+++!
>.
u
0.2 0.3 0.4 0.5 0.6 0.7
Time (seconds)
0.1
O++!t+1++If+4<I++4
o
Figure 4.6: Cylinder Rod Velocity Response for Example 4.1
66
Pressure transients at Nodes 1 2 and 3 are depicted in Figure 4.7. Abrupt change
in pressure are also handled by the ESA When the cylinder rod reaches it stroke limit
the pressures at Nodes 1 and 2 quickly climb until the relief valve opens. Similarl
pressure at Node 3 drops to zero because the cylinder is no longer forcing fluid through
the metering orifice.
":1
0.2 0.3 0.4 0.5 0.6 0.7 ,
Time (seconds)
0.1
r
JlAA P3
Y
PI 
~~
P2
900
800
700
:= 600
til
Q.
 500
I!!
; 400
til
Ql 0: 300
200
100
o
a
Figure 4.7: Pressure Responses for Example 4.1
Example 4.1 was developed and solved in the MATLAB computing environment.
All of the relevant MATLAB script files are contained in Appendix E.
7
CHAPTER 5
CASE STUDY: PILOTOPERATEDRELIEF VALVE
Introduction
The fourpoint method for generating cubic splines was developed in this effort a
a modeling tool. Its primary function is to eliminate model discontinuities which cause
stiffODE solvers and adaptive step size algorithms to fail or become inefficient.
Similarly, an Event Switching Algorithm was developed in this work to effectively
eliminate failures caused by abrupt changes to the Jacobian matrix ofa stiffODE solver,
When combined~ these two techniques are used to create compatibility between
mathematically stiff hydraulic systems containing discontinuities and stiffODE solvers
with adaptive step size control. These stiffODE solvers offer a considerable advantage
over conventional solvers in tenns of processing time.
Overview of StiffODE Solvers
Two stiff numerical integration methods are predominant. These methods are
Gear's method for stiff ODEs and the Rosenbrock method, Ofthese, Gear's method is the
most common. Gear's method is a predictorcorrector method. The predictor and
corrector are based on Backward Differentiation Fonnulas (BDFs) or Numerical
Differentiation Formulas (NDFs) [21]. Though closely related, the NDFs are generally
more efficient than the BDFs [22]. The order of the BDFs or NDFs affects stability.
During the course of integration, the order is varied by the solver [23]. If a Iimit is
applied to the maximum order, greater stability is achieved. The degree of stability
decreases and efficiency increases as the order limit is raised.
II,
",
..,
II
"
""
To explore the effects oforder Example 4.1 was solved using a varie of order
limits. The results are compiled in Table 5.1. A stability problem occurred when th
Table 5.1: Simulation Results for Example 4. 1 Using Different Variation ofG ar'
Method.
Gear's Method: Integration Results: Processing Time•
Solver Parameters Success or Failure (seconds)
NDFs Success 39.08
Maximum Order = 1
NDFs Success 13.98
Maximum Order = 2
NDFs Success 11.72
Maximum Order = 3
NDFs Success 11.54
Maximum Order = 4
NDFs Failure Not Applicable
Maximum Order = 5
* These results were obtained in the MATLAB computing environment using a Pentium 200 computer with
32 megabytes of RAM. The absolute and relative tolerances were set at 10'7 and 104 respectively.
maximum order was set at 5. This stability problem ultimately led to failure of the
method. A review of the processing times reveals a threefold increase in efficiency a
the order limit is increased.
Although it is not as numerically efficient as Gear's method, the Rosenbrock
method has become popular for two reasons. First, the Rosenbrock method is relatively
simple. This method is a singlestep solver based on the familiar RungeKutta scheme.
The need for complicated predictorcorrector fonnulas with varying order is eliminated.
As a result, the method is conceptually easy to understand. Secondly, the Rosenbrock
method has stability properties which surpass those of Gear's method. The Rosenbrock
method can often solve problems which cause Gear's method to fail.
"
'II
'. .:;
":~
Tbe Pilotoperated Relief Valve
In an effort to investlgate the advantages offered b the e stiffODE 01 ers, the
dynamics ofa pilotoperated relief valve were studied. In addition, the costs ofrestarting
the numerical integrator in the ESA were explored by varying the travel limits ofthe
poppets. A typical pilotoperated relief valve is depIcted in Figure 5.1. The valve
Orifice
Port T (0 psij
~
Port P
Figure 5.1: PilotOperated Relief Valve
functions by controlling the force balance of the main poppet. As the pressure
increases, it distributes equally on all surfaces as long as the pilot poppet remains seated.
With the aid of a light spring, the closing force is larger than the opening force and the
main poppet remains against its seat. The valve opens when the pressure becomes large
enough to lift the pilot poppet off of its seat. When the pressure becomes large enough to
overcome the pilot spring, flow is established through the orifice. The orifice creates a
pressure drop in the spring chamber ofthe main poppet. As a result, the opening force
becomes larger than the closing force and the main poppet lifts offof its seat. The bulk
of the flow then passes to tank by flowing around the unseated main poppet.
".:
; I•~
70
Dynamically, this process is described by six first~rder differential equations.
These equations are listed below. Flow forces actjng on the poppets have been included
in this system of equations.
Where,
f3 = Fluid Bulk Modulus of Elasticity
VI = Fluid Volume at Inlet Port
V2 = Fluid Volume of Main Poppet Spring Chamber
Qil7 = Inlet Flow Rate
QI = Flow Rate Over Main Poppet Seat
Q2 = Flow Rate Over Pilot Poppet Seat
QaT = Flow Rate Through Orifice
D) = Diameter of Main Poppet Seat
D2 = Diameter of Pilot Poppet Seat
(5. I)
(5.1)
(5.3)
(5.4)
.~
(5.5)
I,
(5.6) " ,
'h

Ap =
Fpre/oadl =
Fpre/oatll. =
Fflowl =
Area ofMain Poppet Seat
Area ofPilot Poppet Seat
Area of Main Poppet on the Spring Chamber Side
Mass ofMain Poppet
Mass of Pilot Poppet
Force of Main Spring Preload
Force of Pilot Spring Preload
Flow Force Acting on Main Poppet
Flow Force Acting on Pilot Poppet
Damping Coefficient for Main Poppet
Damping Coefficient for Pilot Poppet
Spring Rate for Main Spring
Spring Rate for Pilot Spring
71
For the test problem, the following constants are known:
/3= 1.03 X 109 Pa VI = 3 X 104 m3
V2 = 1 x 107 mJ Q. = 1 X 103 m3/sec /II
Dj =0.017m D2 =0.005 m
AI = 2.27 X 104 m2 • J A2 = 1.96 x 10 m
Ap = 2.35 X 104 m2
ml = 0.045 kg
m2 = 0.020 kg f~reloadl = 100 N
bl =1000 N/(m/sec) b2 = 50 N/(m/sec)
k) = 5000 N/m k2 = 50000 N/m
....
7
The desired cracking pressure (Per) for the relief valve is 1x 107 Pascal. Therefore the
spring preload on the pilot poppet must be:
(5.7)
The flow rate through the fixed orifice (Qor) requires a mathematical model. A suitable
model may be developed by using the procedure outlined in Part B ofExample 3.1. In
this case, the following information is given:
Orifice Diameter (Dor) = 0.001 m
Discharge Coefficient (Cd) = 0.6]
Fluid Density (p) = 845 (kg/m3
)
Kinematic Viscosity (v) = 14.3 X 106 m2/sec
Transition Reynolds Number (NR1) = 1500
The following orifice flow model may be developed from this information:
For 0 ~ (PIP2) ~ 5.224 x 10 5 Pa.
Qor = 3.759 X 1023 (PI_P2)3  7.0]4 X 1017 (P 1P2i + 5.863 x 1011 (1"',1'2) (5.8)
For (PIP2) >5.224 x 10 5 Pa.
Q" ~ c~o, ~~(Ii  p,) (5.9)
The flow rate over the main poppet seat (QI) is essentially the flow through a variable
orifice created by the moving poppet. The flow area of this variable orifice is a function
of poppet displacement (Eq. 5.10).
(5.10)

T'
=
Flow Area ofMain Poppet/Seat
Diameter of Main Poppet Seat
Half Angle of Main Poppet
Main Poppet Displacement
The mathematical model for this variable orifice may be developed by first assuming a
fixed orifice with a diameter equal to the diameter ofthe main poppet seat (D!). Using
the procedure outlined in Part 8 of Example 3.1, a cubic spline may be generated. The
variability of the orifice is handled by simply multiplying this cubic spline by the ratio of
the flow area (Ajlowd to the seat area (AI). In this case, a cubic spline may be developed
using the following information:
Seat Diameter (DI ) = 0.017 m
Discharge Coefficient (Cd) = 0.61
Fluid Density (p) = 845 (kg/m3
)
Kinematic Viscosity (v) = 14.3 X 106 m2/sec
Transition Reynolds Number (NR1) = 1500
The resulting cubic spline is:
(5.11 )
"
:1
":1
;1
"
I,I'
Ii '."
Given a main poppet half angle of 30° (aJ = Tt/6 rad.), the variable orifice model is:
For 0 ~ (Pd ~ 1.808 x 10 3 Pa,
mJ) sin(a l Jx 14 3 II' 7 Ql= {1.543 x 10 (Pd 9.959xlO (Pd~+2.88)xlO (Pd}(5.12)
AI
74
For (PI) >1.808 x 10 3 Pa,
(5.13)
The mathematical model for the variable orifice created by the pilot poppet rna be
generated in the same fashion as that ofthe main poppet. The data for the pilot poppet is
as follows:
Seat Diameter (D2) = 0.0005 m
Discharge Coefficient (Cd) = 0.6]
Fluid Density (P) = 845 (kglmJ
)
Kinematic Viscosity (v) = 14.3 x 1O~ m2/sec
Transition Reynolds Number (NRt) = ]500 II
Q2 = ;rD2 s:(a2 )y {2.937 x 1018 (P2)J  2.192 X 1013 (P2)~ + 7.329 X 10<) (J12)1
2
Pilot Poppet Half Angle (a2) =
The resulting mathematical model is:
For 0 S; (P2) S; 2.090 x 10 4 Pa.
1C/6 radians
(5.14 )
For (P2) >2.090 x 10 4 Pa,
(5.15)
75
With the mathematical models in place, the discontinuity functions can b
identified. The pilotoperated relief valve has two springmassdamper stems. Both of
the masses may encounter state events in the fonn of mechanical stops. As a result the
ESA must be implemented twice. After each time step, the appropriate discontinuity
function for both masses is evaluated. If either discontinuity function becomes negative,
the integrator restarting procedure is initiated. For reasons discussed in Chapter 4 of this
document, Equations 5.4 and 5.6 will be used by the ESA as the discontinuity functions
for their respective poppets when the poppets are pinned against mechanical stops. If
either of the poppets is between its mechanical travel limits, the displacement (X3 or xs) is
monitored until a mechanical stop is reached as dictated by the ESA.
Numerical Integration Evaluation
The pilotoperated relief valve was inserted in a simple test circuit as depicted in
Figure 5.2. Shifting the directional control valve at a predetermined time produced a
~
Reservoir
Figure 5.2: PilotOperated Relief Valve Test Circuit
7
duty cycle consisting of an "on" response and an "off' response. A series oft :ts with
varying travel limits on both poppets was perfonned to evaJuate the effectiveness of the
stiffODE solvers and the ESA. In addition to Gear's Method for StiffODEs and
Rosenbrock's Method, three conventional (nonstift) solvers were evaluated. An attempt
to solve each test problem was made using all ofthe following ODE solvers:
1) RungeKutta (4,5): A singlestep nonstiff soJver credited to
Dormand and Prince [22]. This ODE solver uses 4th and 5th order
RungeKutta methods.
2) RungeKutta (2,3): A singlestep nonstiffsolver credited to
Bogacki and Shampine [22]. This ODE solver uses 2nd and 3rd
order RungeKutta methods.
3) AdamsBashforthMoulton [22]: A nonstiffmultistep solver. This
ODE solver uses a variable order routine to iteratively predict and
correct at each time step.
4) Gear's Stiff Method (22): A variableorder multistep solver based
on BDFs or NDFs. For comparison, the following combinations
were tested:
a) BDFs with order limit of2.
b) BDFs with order Iimit of 5.
c) NDFs with order limit of2.
d) NDFs with order limit of 5.
7
5) Rosenbrock's Method [22): A second order, singlestep stitTODE
solver based on semiimplicit RungeKutta Fonnulas.
TEST CONDITIONS:
Each ofthese ODE solvers is available in the MATLAB computing package.
MATLAB provided a convenient environment to implement the ESA and test a variety of
ODE solvers. The following test conditions prevailed for each test problem:
Computer Specifications: Pentium 200 with 32 megabytes ofRAM
Relative Error Tolerance: 1 x 103
Absolute Error Tolerance: 1 x 1O{i
Adaptive Step Size Control:
Event Switching Algorithm:
Initial Conditions:
Initial Time:
Final Time:
Active for all ODE solvers
Utilized for all tests
Zero
t = 0 seconds
t = 0.08 seconds
Initial Direction Control Valve Position:
Directional Control Valve Shift Time:
Open to Rehef Valve
t = 0.04 seconds
The goal of these tests was to explore the potential advantages of using stitTODE solvers.
For fixedstep integration, the nonstiff solvers may not require mathematically smooth
models or an event switching routine because no Jacobian matrix is present. Adaptive
step sizing, however, requires both. Because adaptive step size routines are
commonplace, adaptive step sizing was applied to all of the tests. This condition
provided commonality of mathematical models and event switching. For a given
7
problem, therefore, any variations in performance among the various ODE solvers were
solely due to the solver routines.
TEST PROBLEMS:
Travels limits for each test problem and the corresponding integration results ar
presented on the following pages. The pressures and poppet displacements were of
primary interest As such, only the transient responses of these parameters were
displayed graphically. In each case, the simulation results were identical (within the
specified error tolerance) for all of the integration methods. Table 5.2 contains
performance information for all ofthe tests. The number ofODE solver restarts required
by the ESA and the computation times for the various ODE solvers are also contained in
this table. Computation times presented in Table 5.2 are an average of three runs. As
expected, the difference in elapsed time between the three runs was negligible.
7
Test #1:
Main Poppet Maximum Travel Limit (x"rar) = 3 X 104 meters
Pilot Poppet Maximum Travel Limit (Ymar) = 0.5 X 104 meters
Test #1 Solution:
2.0E+07
1.5E+07
....... 11.!
1.0E+07 ::J
III
III
f
11.
5.0E+06
f"' Pl
f1 I"..
:'
P2
O.OE+OO
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
Figure 5.3: Pressure Response for Test #1
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
TIme (seconds)
M ~n Po pet
I \
I \
I Pi pt POJ pet , \ OE+OO
o
5ED4
eS
4E04
GI' gi
3E04
E•
U
III a. 2E04 •C 11E04
o
11.
Figure 5.4: Displacement Response for Test #1
80
Test #2:
Main Poppet Maximum Travel Limit (xmm) = 3 x 104 meters
Pilot Poppet Maximum Travel Limit (Ymar) = 0.8 X 104 meters
Test #2 Solution:
2.0E+07
1.5E+07 1lI e:..
.f::,:.J 1.0E+07
III
f!
a..
5.0E+06
PI
~ P2
O.OE+OO
o 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
Figure 5.5: Pressure Response for Test #2
it
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
M in Po pet
1\
...... ~ \ >,. • .. UP} '.... \ OE+OO t
o
5EQ4
!
.! 4EQ4
CII
.§.
~ 3E04
E
II
U
1lI
Q. 2EQ4 •o
..!. 1E04
Q. o
a..
Figure 5.6: Displacement Response for Test #2
81
Test #3:
Main Poppet Maximum Travel Limit (xmar) = 4.5 X 104 meters
Pilot Poppet Maximum Travel Limit (Ymar) = 0.8 x 10"" meters
Test #3 Solution:
2.0E+07
1.5E+07 III e::.
$ 1.0E+07
III
III e~
5.0E+06
A
~ I
PI V'' I
P2
O.OE+OO
a 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
Figure 5.7: Pressure Response for Test #3
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
r.
t\ lain ~
J'v \/'
'\
\
\ ;\) POP) et j .\ ~
OE+OO
o
5E04
f
~ 4E04
CII gi
3E04
E•u.. a 2E04
III
Ci
1E04
Q. o~
Figure 5.8: Displacement Response for Test #3
8
Test #4:
Main Poppet Maximum Travel Limit (x1Tun') = 4.75 x 104 meters
Pilot Poppet Maximum Travel Limit (Ymar) = 0.5 X 104 meters
Test #4 Solution:
2.0E+07
1.5E+07 llJ
~e
1.0E+07 ~
II)
II)
.I.».
Q..
5.0E+06
f\,
'r\
l~ PI rv
P2
O.OE+OO
o 0.01 0.02 0.03 0.04 0,05 0.06 0.07 0.08
Time (seconds)
Figure 5.9: Pressure Response for Test #4
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Time (seconds)
n l'fainF oppet
"
I / \f\('....,r.
I
V \
I \
I Pile ~ POPI et g \r~ \ OE+OO
o
5E04
f
.! 4E04
I»
§.
~ 3E04
E
II
Co) •a 2E04 •is
i 1E04
Q.. o
Q..
Figure 5. 10: Displacement Response for Test #4
Table 5.2: Compiled Test Results for Variations ofthe PilotOperated Relief
Valve Problem
Computation Times (seconds)
Test Travel Number NonstiffODE StiffODE Solvers
Number Limits Of Solvers
(meters) Restarts Runge RWliIe Adams Gear's Gear's Gear's Gear's Rosen
Kutta Kutta Bashforth BDFs BDFs NDFs NDFs brock
Xma..x (4,5) (2.3) Moulton Ma~ Ma~ Max Max
Ymax
Order =2 Order =5 Order =2 Order "'5
1 3.0x 10'" 8 239.9 373.8 346.5 10.9 10.4 10.1 10.0 14.9
0.5 x 104
2 3:0x 10'" 8 240.4 335.0 346.8 16.9 13.2 15.2 12.9 24.2
0.8 x 10'"
.... ' 4.5 x ]0'" 8 264.0 335.0 369.8 21.2 13.7 18.6 IJ8 32.6 0.8 x 10'·
4 4.75 x 10" 12 265.0 334.0 360.8 21.8 17.3 19.8 17.4 32.6
0.5 x 10'"
Discussion of Computation Time Results
The effectiveness ofthe various ODE solvers can be evaluated by studying Table
5.2. In the presence of mathematical stiffness, RungeKutta (4,5) was the most effective
of the conventional nonstiffODE solvers. The stiffODE solvers were significantly more
efficient. At worst, an eightfold improvement was achieved when using a stiffODE
solver. Processing speeds were over 20 times faster than that of RungeKutta (4,5) in
some cases. With the aid ofcubic splines generated using control points and the Event
Switching Algorithm, these computational efficiency gains may be realized. The stitT
solvers would have failed due to abrupt changes in their Jacobian matrices if these tools
were absent. As such, the effort required to implement these tools is well justified. This
is especially true of very large systems which could take days to inteblTate.
'4
Among the stiff solvers, all ofthe variations ofGear's method proved to be
slightly more efficient than Rosenbrock's method. This condition is problem specific.
There exists a set of problems for which Gear's method is ineffective [22]. However th
inherent stability ofRosenbrock's method will lead to a successful solution. The ability
of Rosenbrock's method to solve a wider variety of problems oft,en compensates for its
slightly longer processing time.
Table 5.2 also reveals some information about the variations of Gear's method.
Numerical differentiation fonnulas of a given order are generally more efficient than
their corresponding backward differentiation formulas. In the two cases where the NDFs
did not perfonn better than the BDFs, the computation times are nearly equal. For both
BDFs and NDFs, an increase from 2nd order to 5lh improved the computation time by as
little as I% or by as much as 35% depending on the problem. Again, stability is the key
issue. Differentiation formulas of 2nd order maintain Astability [22]. However, this
stability property is lost if the differentiation fonnula is above order two. For the pilotoperated
relief valve problem, stability is maintained through 5th order.
The Event Switching Algorithm performed eight ODE solver restarts each during
test numbers I, 2 and 3. Despite the equal number of restarts, computation times varied
significantly between these three tests for a given stiffODE solver. These results suggest
that restart costs are not the dominant factor in determining processing time. A careful
review ofthe solutions reveals the true cause of the varying computation times.
Inspection of Figure 5,4 reveals a very simple dynamic response. Both poppets spend
most of the response time against their mechanical stops during Test #1. During Test #2,
8
the pilot poppet spends only a short amount oftime against the upper tra ellimit (Figure
5.6). The dynamics ofTest #2 are slightly more complex than tho e ofTest #1. As a
result, slightly more computation time is required The solution for Test #3 (Figure 5.8)
is clearly more oscillatory than either Test #1 or #2. Again, the increased dynamic
complexity is accompanied by longer processing times.
Test numbers 3 and 4 may be compared to examine the effects of restarting the
various stiffODE solvers. Although dynamically similar, Test #4 requires four more
restarts than Test #3. The 2nd order BDFs for Gear's method required only a 2.8%
increase in processing time to accommodate the four additional restarts. Similarly, 2nd
order NDFs required 6.5%. However, the 5th order BDFs and NDFs required over 26%
more processing time to accomplish the same task. The costs of restarting the low order
differentiation formulas for Gear's method are clearly much lower than those ofthe high
order formulas. By comparison, Rosenbrock's method required no significant
computation time to perfonn the additional restarts. This fortunate result i inher nt to
the method itself Rosenbrock's method is a singlestep solver which requires no
historical information. As a result, the first step is virtually the same as any other step.
Gear's method is a multistep solver which requires historical information to complete a
step. When the solver is started, no historical information is available and special
measures must be introduced to start the solver.
Error Analysis
In an effort to investigate the reliability of the stiff solvers in terms of accuracy,
one final test was performed. Test #4 was repeated using RungeKutta (4,5) and a step
86
size limit of )07 seconds. This step size limitation is several orders of magnitude smaller
than the minimum step size produced by the adaptive step size control algorithm with no
size restrictions. Because accuracy is directly related to step size, the simulation
performed at 107 was considered to be "exact" for all practical purposes and s rved as a
baseline for comparison. The results from the inlet pressure (PI) response of Test #4
using the various ODE solvers were compared to the results obtained using RungeKutta
(4,5) with a step size restriction of 107 second. The results were refined to provide state
values at increments of 104 second. At each increment, the pressure difference between
the results from the solver of interest and the "exact" solution was calculated. This
difference was used to determine error percentage as follows:
IPsc/v",.  Pt!:U>C/ I Percent Error = x 100%
PUQC/
(5.16)
Vvllere, Inlet Pressure (PI) as computed by the ODE solver of
interest with no step size restriction.
Pe:rac/= Inlet Pressure (PI) as computed by RungeKutta (4,5) with
a step size restriction of 10.7 second.
The results of this error analysis, in terms of maximum error and average error are
compiled in Table 5.3.
Ta ble 53. E AnI'R IfiIIP R fT #4 . . rror alYSlS esu ts or net ressure es lonse 0 est
Numerical Integration Maximum Average
Method Error Error
RungeKutta (4,5) 0.0294 % 0.0024 %
AdamsBashforthMoulton 0.1930% 0.0077 %
Gear's:BDFs Max. Order =2 0.3932 % 0.0696 %
Gear's:BDFs Max. Order =5 0.3539 % _..__. 0.0522 %
Gear's:NDFs Max. Order =2 1.0272 % 0.1124 %
Gear's:NDFs Max. Order =5 0.1709 % 0.0434 %
Rosenbrock's Method 0.0756 % 0.0140 %
7
A graphical depiction ofthe error analysis results is contained in Figure 5. ] 1.
Adams Gear's Gear's Gear's Gear' Rosenbrock
Bash/orth BDFs BDFs NUFs NDFs
Moulton Max. Max. MELx. Max.
Order =2 Order =5 Order =2 Order =5
ID Maximum Error • Average Error I
,.....

~ ~ ....... It • .. IlLn a RlUlge
KutUl
(4,5)
0.2
~
.e.O.8 w
~O.6
Cb u.. ~O.4
1.2
1
Figure 5. J]: Comparison of Error Analysis Results
Discussion of Error Analysis Results
The worst case error was slightly over] %. Accuracy of this magnitude is more
than adequate for hydraulic pressure calculation. RungeKutta (4,5) wa the most
accurate of all the solvers but this accuracy is achieved at the expense of processing time.
The Rosenbrock method provides accuracy comparable to RungeKutta (4,5) without
excessive processing time. Among the various fonns ofGear's Method, Numerical
Differentiation Formul.as (NDFs) with a 5th order limit were the most accurate. This fonn
of Gear's method is also the most efficient in tenns of processing time making it a
powerful integration method for stiff hydraulic systems. Ultimately, stability will
determine the most suitable integration method because adequate accuracy may be
achieved using any of the stiff ODE solvers proposed in this study.
Summary
The optimum ODE solver for hydraulic systems is problem specific. However
Rosenbrock's method would be the most useful to the average hydraulic system designer.
Although not as computationally expeditious as Gear's Method, Rosenbrock's method is
considerably more efficient than the nonstiff solvers. The Rosenbrock method is al 0
conceptually easier to understand than Gear's method and its robust stability properties
facilitate the solution ofa wider variety of problems. In addition, the restart costs of
Rosenbrock's method are insignificant. This characteristics makes it compatible with the
Event Switching Algorithm. For systems with a large number of components, these
restart cost could become excessive if Gear's method were used.
89
CHAPTER 6
EXPERIMENTAL VERIFICATION
Introduction
The ultimate goal of computer modeling and simulation of hydraul ic systems is to
determine the dynamic response of"realworld" hardware. This response information is
invaluable to hydraulic system designers. In the computer environment, the effect of
design changes on the dynamic response can be detennined without the financial
consequences of complicated "breadboard" testing. However, computer simulation can
not completely replace laboratory tests. Computer analysis results are only as accurate as
the methods used to simulate the physical system. As a result, verification of computer
simulation results must be performed at some point in the design process. In this spirit,
the dynamic response of an actual directacting relief valve was simulated using the
modeling and integration techniques present in this work and then compared to available
test data.
Test Component
The hydraulic component under consideration is a directacting relief valve. A
typical directacting relief valve is represented in Figure 6.1. The flow through the valve
remains at zero until the force acting on the poppet is sufficient to lift the poppet from its
seat. Increasing pressure acting on the poppet area at Port P creates this force. As the
poppet lifts from the seat, a flow path from Port P to Port T is created and the pressure
relieved. A relief valve is typically used to limit hydraulic system pressure. Relief valves
o
protect other hydraulic system components from accidental overpre surization and
prevent injuries related to overpressurization.
Port T
Port P
Figure 6.1: Typical DirectActing Pressure Relief Valve
The dynamic performance of a relief valve is of paramount importance for
obvious safety reasons. As a result, performance testing of relief valves is common and
data from such tests is readily available. The availability of such test data makes the
directacting relief valve a perfect candidate for experimental verification of computer
simulation results.
In this case, data from a relief valve test performed at FES Incorporated in
Stillwater, Oklahoma was used. The test unit was a commercially available directacting
91
relief valve. Testing was performed with MILH5606 hydraulic fluid at 38 Celsiu . As a
result of laboratory work, the following valve characteristics were identified:
Test System
Poppet Mass (m):
Spring Rate (k):
Damping Coefficient (b):
Poppet HalfAngle (a):
Seat Diameter (d):
Cracking Pressure (Perl
Fluid Volume at Port P (v):
Discharge Coefficient of
Poppet/Seat (Cd):
Bulk Modulus of Fluid (/1):
Density of Fluid (p):
0.0031 kg
105076 N/m
35 N/(m/sec)
7.75 X 103 m
0.61
For testing purposes, the test unit was installed in the hydraulic system depicted
schematically in Figure 6.2. Shifting the directional control valve with an electrical
Plug
1
~
Reservoir
Figure 6.2: DirectActing Relief Valve Test Circuit
step input sends a flow rate (Qin) of6.305 x 104 m3/sec to the test valve. Ho\
input flow rate does not occur when the directional control val e is shifted. Full.
developed flow occurs over a finite period of time. In an effort to account for thi
condition, a 15 millisecond ramp to full flow was utilized for Qin.
Mathematical Models and Differential Equations
r a step
The mathematical model for variableorifice flow and the differential equations
for the directacting relief valve were developed in the same manner as those for the
pilotoperated relief valve in Chapter 5 of this study. The resulting differential equations
are as [0110\\'5:
(6.1 )
(6.2)
(6.3)
Where,
X I = Pressure at Port P
X2 = Poppet Displacement
X3 = Poppet Velocity
Q2 = Flow Rate Across Poppet Seat Arrangement
Results of Computer Analysis
A cubic spline based orifice model was used for the poppet/seat arrangement and
the Event Switching Algorithm was employed for boundary handling. Numerical
3
integration ofthe differential equations was perfonned using Gear s Method for Stiff
ODEs. Numerical Differentiation Fonnulas with a maximum order of 5 were utilized.
For a relief valve, the parameter of interest is the pressure at Port P. The dynamic
response for this pressure, as predicted by employing the methods outlined in this work,
is displayed in Figure 6.3. The following infonnation can be gleaned from this plot:
Peak Pressure:
Peak Time:
Steady State Pressure:
1.513 x 107 Pa
0.0068 seconds
8.23 X 106 Pa
~ 
/
/
i /
V ,
1.SE+07
1.4E+07
Ci'1.2E+07
f!:.1.0E+07
CI)
:5S.0E+06
Cf)
mS.OE+OS
r..
D..4.0E+06
2.0E+06
O.OE+OO
o 0.025 0.05 0075
Time (sec)
0.1
Figure 6.3: Simulation Results for DirectActing Relief Valve
Test Results
The dynamic response from laboratory testing is displayed in Figure 6.4. Again,
the perfonnance characteristics may be obtained from the plot:
Peak Pressure:
Peak Time:
Steady State Pressure (Average):
1.552 X 107 Pa
0.008 seconds
94
f\
I
I ~ l ......... ~I~ ~ r...
~ .  IV ~
\
,
V" ,
1
1.SE+07
1.4E+07
C;;1.2E+07
e:.1.0E+07
(I)
S8.0E+06
:'."C.S.OE+OS
Q. 4.OE+OS
2.0E+06
O.OE+OO
o 0.025 0.05 0.075
Time(sec)
0.1
Figure 6.4: Test Results for DirectActing Relief Valve
The computer simulation result is superimposed over the measured dynamic
response in Figure 6.5.
~
I
Test Data
v
JlI
~ rcr'\. ~  .
r/ . rv 'V
""I
I
~ ~;Simulation
V Data
I
1.SE+07
1.4E+07
_1.2E+07
e":'.1.0E+07
(I)
S8.0E+OS
'"f6.OE+OS
Q. 4.0E+06
2.0E+06
O.OE+OO
o 0.025 0.05 0.075
Time (sec)
0.1
Figure 6.5: Simulation and Test Results Superimposed
9
Discussion of Results
A comparison ofthe test results to the computer simulation results reveals
excellent correlation in tenns of peak pressure. The difference between the two value IS
less than 3%. Peak pressure is the most important characteristic of relief valve
performance. Pressure spikes are often responsible for hydraulic system damage. The
predicted steadystate pressure value also closely matches the experimental data.
Although steadystate analysis may be performed without numerical integration, the
steadystate value is a good measure of correlation between computergenerated results
and empirical results. Comparing the peak times reveals some difference in terms of
response time. The predicted response is about 0.00 I7 seconds faster. This slight
disparity is due the way the system was modeled. The dynamics of the pressure
transducer were not considered in the system model. As a result, the measured response
lags the predicted response. A more complex model could be developed to account for
instrumentation dynamics if greater response time accuracy is desired.
CHAPTER 7
CONCLUSION AND RECOMMENDAnONS
Conclusion
Hydraulic systems are generally considered to be mathematically stiff. This
stiffness is the direct result of widely varying time constants within a single system.
During numerical integration of a stiff set of ordinary differential equations,
mathematical stiffness forces conventional ODE solvers to use very smaJi time steps to
maintain stability. These small time steps lead to a considerable amount of
computational effort and processing time.
In an effort to reduce processing time, ODE solvers designed specifically for stiff
problems have been devised. In addition, adaptive step size routines have been
developed to improve computational efficiency. Hydraulic systems often contain
discontinuities which cause both ofthese advancements to fail or become inefficient.
Stiff ODE solvers use a Jacobian matrix to perpetuate a solution. Abrupt changes to the
Jacobian matrix due to discontinuities will cause the stiff solver to fail. Similarly, abrupt
changes will cause adaptive step size algorithms to "hunt" around a discontinuity until the
step size is reduced enough to traverse the problem area. This hunting process is
excessively time consuming.
In order to solve the problems caused by discontinuities, it is necessary to identify
the types of discontinuities found in hydraulic systems. Discontinuities may exist in
mathematical models or they may be encountered during the integration process. Model
discontinuities typically result from slope discontinuities in the mathematical model.
97
The lack of mathematical smoothness causes abrupt changes in the Jacobian matri and
results in solver failure. Discontinuities encountered during numerical integration ha e
the same disastrous effects. This type of discontinuity is termed a state nt. In
hydraulics, state events occur when mechanical devices encounter travel limits. For
example, the differential equations describing a system will change if an actuator reaches
a mechanical stop.
In this study, two powerful techniques have been developed to overcome the
problems associated with these discontinuities. Cubic splines generated with a fourpoint
method were created to smooth mathematical models. In addition, an Evenl Switching
Algorithm (ESA) was developed to provide seamless integration over state events.
The fourpoint cubic spline technique provides a means to join two discontinuous
functions. This technique involves positioning a control point in such close proximity to
the true end point of a curve that the difference is physically insignificant. However, the
relative position of the control point and the true end point determines the shape of the
resulting curve. Control points may be used to force the slope continuity of a cubic
spline bridging a point of discontinuity. This new technique has proven to be both
versatile and effective. In addition to bridging discontinuities, cubic splines may be
implemented to eliminate the infinite stiffness assoc1ated with infinite slopes.
The ESA provides a means to transcend state event discontinuities. Successful
integration is accomplished by restarting the ODE solver when a state event is
encountered. This restarting procedure eliminates abrupt changes to the Jacobian matrix.
State events are defined by discontinuity functions. The exact time ofa state event
8
occurs when a discontinuity function is equal to zero. For hydraulic systems these
discontinuity functions are based on the forces applied to a mass and on its displac ment.
After each step, the appropriate discontinuity functions are evaluated. If a stat event is
encountered (as signaled by a sign change), an event location routine is invo ed to
determine the precise time ofthe state event. Subsequently, the ODE olver is restarted
with the new set of differential equations. This technique has proven effective in the
successful integration of systems containing actuators or valve components which
t:ncounter physical travel limits.
With the aid ofthe cubic spline modeling tools and the ESA, hydraulic systems
are compatible with ODE solvers specifically designed for stiff differential equations. Of
these solvers, Gear's method and Rosenbrock's method are predominant. Gear's method
has greater potential for computational efficiency. Tests performed in this study reveal
processing times over 20 times faster than conventional (nonstifl) ODE solvers. While
not as efficient as Gear's method, Rosenbrock's method has superior stabil1ty properties.
As a result, Rosenbrock's method can solve a wider variety of problems.
Suggestions for Further Study
The value of curves developed by the fourpoint method and ofthe ESA have
been demonstrated in this study. Significant gains in processing speed were realized by
implementing these techniques. However, further study would refine and improve these
techniques. A list of topics for further study is presented below
1) The shape of the cubic spline used to bridge a model discontinuity may effect
processing speed. Curves approaching sharp corners would likely slow the ODE solver
significantly. However, these sharper curves would more accuratel represent the
original discontinuous model. A study of processing speed versus model accurac. could
give insight into this relationship.
2) The ESA uses on event location scheme to pinpoint the time at which a tate
event occurs. Several methods for event location are available. Among these methods
are bisection, false position, secant linear interpolation, inverse quadratic interpolation,
and the Illinois method [18]. A comparison study of these methods may reveal
advantages in terms of computational efficiency or stability. In this effort, the Illinois
method was used exclusively.
3) If the ESA is applied to multiple objects in a single system it is possible for more
than one state variable to cross its threshold value during a single time step. An event
location scheme designed to handle this condition would increase the versatility of the
ESA.
10
Works Cited
[IJ Ellman, Asko, Robert Piche', and MattI Vilenius. "Integration of Numerically
StifT Fluid Power Circuit Models Using an LStable RungeKutta Method." Paper
presented at 7th Bath Inti. Fluid Power Workshop, September 1994.
[2J Ellman, A.u. "Proposals for Utilizing Theoretical and Experimental Methods in
Modelling TwoWay Cartridge Valve Circuits." Ac/a Polylechnica MEl 0 I. 1992
pp.257.
[3J Bowns, D. E., S. K. Dugdale, and S. P. Tomlinson. "Progress Towards a General
Purpose Hydraulic System Simulation Language." Proe. 6th IntI. Fluid Power
Symposium, 1981, pp.115131.
[4] Bowns, D. E., and L. M. Wang. "The Digital Computation of Pressures in
Hydraulic Pipes with Small Volume Using an Iterative Technique." Froe. Insln.
Meeh. Engrs., 1990, Part C, Vol. 204, pp. 2936.
[5] Ellman, A., and R. Piche'. "Numerical Integration of Fluid Power CirCUIt Models
Using TwoStage SemiImplicit RungeKutta Methods." Proe. Ins/n. Meeh.
Engrs., 1994, PartC, Vol 208, pp. 167175.
[6] Krus, P. "The Simulation of Fluid Power Systems with Complex Load Dynamics."
International Journal ofModelling and Simulation, 1986, Vol. 6, No.2, pp. 5257.
Pl Krus, P.,and J. O. Palmberg. "Simulation of Fluid Power System in Time and
Frequency Domains." Proe. 7th Inti. Fluid Power Symposium, 1986, pp. 7379.
[8] Bowns, D. E., R. E. Dorey, and S. P. Tomlinson. "Computer Simulation
Techniques for the Dynamic Performance Assessment of Fluid Power Systems."
Proc. 71h Intf. Fluid Power S~vmposium. 1986, pp.8188.
[9] Ellman, A. u., and M. J. Vilenius. "Methods for Simulating the SteadyState and
Dynamic Behavior of TwoWay Cartridge Valve Ci.rcuits." SAE.!. ofCommercial
Vehkles. 1990, Vol. 99, No.2, pp. 384393.
[10] Carver, M. 8., and S. R. MacEwen. "Numerical Analysis ofa System Described
by implicitlyDefined Ordinary Differential Equations Containing Numerous
Discontinuities." Appl. Math. Modelling, December 1978, Vol. 2, pp. 280286.
[IIJ Ellison, D. "Efficient Automati.c Integration of Ordinary Differential Equations
with Discontinuities." Mathemalies and Computers in SimulaJion, 1981, Vol.
XXIII, pp. 1220.
101
[12] Caney, K. "integration Across Discontinuities in Ordinary Differential Equations
Using Gear's Method." University ofBath. School ofEngineering R port I ,
1979.
[13] Berzins, M., and A. 1. Preston. "Algorithms for Location of Discontinuities in
Dynamic Simulation Problems." Computers Chem. Engng, 1991, Vol. 15. No. 10,
pp.701713.
[14] Crosbie, R. E., and J. L. Hay. "Digital Techniques for the Simulation of
Discontinuities." Proc. ofthe 1974 Summer Computer Simulation Conference,
1974, pp. 879].
[15J Carver, M. B. "Efficient Integration Over Discontinuities in Ordinary Differential
Equation Simulations." Mathematics and Computers in Simulalion, 1978, Vol.
XX, pp. 190] 96.
[]6] O'Regan, P. G. "Step Size Adjustment at Discontinuities for Fourth Order RungeKutta
Methods." The Computer Journal, 1970, Vol. 13, No.4, pp. 401 404.
[17] Chaplin, RI., R.E. Crosbie, and 1. L. Hay. "Integration Routines for Systems with
Discontinuities." The Computer Journal, ]974, Vol. 17, No.3, pp. 275278.
[18] Moler, C. "Are We There Yet." Matlah News and NOles, Simulink 2 ,",pecial
Edilion, 1997, pp. 1617.
[19] Gerald, C. F., and P. O. Wheatley. Applied NumerIcal Analysis. 5th ed. Reading,
MA: AddisonWesley 1994, pp. 233244.
[20] Fitch, E. c., and l. T. Hong. Hydraultc Component Design and Selection.
Stillwater, OK: Bardyne 1997. pp. 3742.
[21] Lapidus, L., and W. E. Schiesser. Numerical Method\Ior Differential ..<"ystems:
Recent Developments in Algorithms. Software, and ApplicQlions, New York:
Academic Press 1976, pp. 97124.
[22] Using Mat/ab. Natick, MA: The Mathworks, Inc. ]997. pp. 8] to 852.
[23] Reichelt, M. W., and L. F. Shampine. "The Matlab ODE Suite" SIAM./. Sci.
('amput., January 1997, Vol. 18, No. ], pp. 122.
Appendix A
Procedure for Generating Cubic Splines
to Interpolate Between Four Points
10.
10
In matrix fonn, three cubic spline polynomials connecting four points rna be gen rated
as follows [19].
]) Define the four points in the points in Cartesian coordinates:
(A. I)
2) Create X and Ymatrices:
XI YI
x=
x 2 Y= Y2 (A.2)
x3 y,
xl Yl
3) Generate H matrix:
4) Create the RHS (right hand side) matrix:
(A.3)
RHS= (A4)
5) Build the A matrix:
0 0 0
h) 2(111 +112 ) h2 0
A= (A.5) a h2 2(~ + hJ ) h)
a a 0
6) Solve set of equations to get S matrix:
104
7) The three cubic spline equations are:
Between point 1 and point 2:
(A.6)
(A.7)
Beh.veen point 2 and point 3: (A.S)
Between point 3 and point 4; (A.9)
Appendix B
Experimental Determination of the Pressure
and Flow Relationship of an Orifice
10
Test System and Equipment
The hydraulic system used for this test is depicted if Figure B1. In addition to
the hydraulic system, a stop watch and graduated cylinder were required.
106
Pump
Heat
Exchanger
Orifice
V
Flow
Cantlal
Valve
GrOdUOled
CYlinder Water
Figure B1: Hydraulic System for Orifice Testing
Test Conditions
Test Fluid: SAE lOW Oil
Fluid Temperature: 73 of
Orifice Dimensions: Diameter = 0.012 inches Thickness = 0.010 inches
107
Test Procedure
1) The test system was filled with the specified test fluid.
2) The bypass valve was opened.
3) The pump drive motor was started.
4) Water flow rate through heat exchanger was adjusted to achieve the
desired fluid temperature.
5) The bypass valve and flow control valve were adjusted to obtain
incrementally larger pressure differences across the orifice.
6) At each increment, the flow rate through the orifice was measured and
recorded using a graduated cylinder and a stop watch.
Test Results
Table B1: Test Results for 0.012 Diameter Orifice
Pressure
I
Flow Rate
Difference (inJ/sec.)
(psid)
0.0 0.0000
2.5 0.0069
5.0 0.0135
7.5
10.0
12.5 0.0288
15.0 0.0335
17.5 0.0374
10
20.0
22.5
25.0
30.0
0.0412
0.0439
0.0474
0.0539
35.0 0.0593
40.0 0.0643
45.0 0.0685
50.0 0.0743
55.0 0.0801
60.0 0.0839
65.0 0.0893
70.0 0.0936
75.0 0.0970
80.0 0.0997
85.0 0.1024
90.0 0.1055
95.0 0.1084
100.0 0.1114
105.0
110.0
115.0
120.0
0.1141
0.1168
0.1196
0.1223
Appendix C
MATLAB Script File for Example 3.1
10
%File Name: sp8ex.m
%File Description:
% MATLAB Script file used to generate cubi spline for
~ orifice model based on experimental data (Example 3.1).
% 0.012 orifice; SAE lOW @73 F
%%DECLARE EMPTY MATRICES FOR PLOTTING PURPOSES
xv= [] ;
xp= [J;
xp1=[] ;
xp2=( ) ;
gx= [] ;
gxl=[ J;
gx2=[] ;
*%INITIALIZE VARIABLE CONSTANTS
n=4; %nurnber of points for cubic spline
dia=O.Ol2; %orifice diameter
area=((dia~2)*pi)/4i %orifice flow area
ro=8.l7le5; %fluid density (lbfsec~2/in~4)
v=.1623; %kinematic viscosity of fluid (in~2/sec)
cd=.63i %orifice discharge coefficient
cf=3.85i%conversion factor  in~3/sec to GPM
%% DEPARTING SLOPE AT ORIGIN
slope=.00220
%% GENERATE CONTROL POINT AT ORIGIN
%% (xl,yll True end point
~% (x2,y2) = Control point
xl=O;
yl=O;
epsx=l i
x2=xl+epsx*eps;
y2=slope*x2i
while( (y2yl)<eps)
x2=xl+epsx*eps:
y2=slope*x2i
epsx=epsx*lO;
end
~~ GENERATE CONTROL POINT AT 85psid
• ~ (x3, y3) True end point
%% (x4,y4) = Control point
x3=85; t Pressure difference at second end point
y3=cd*area*sqrt(2/rol*sqrt(x3) i tflow rate at 85 psid
epsx=li
x4=x3+epsx*epsi
y4=cd*area*sqrt(2/rol*sqrt(x4) :
whOle((y4y3)<eps)
x4=x3+epsx*eps;
y4=cd*area*sqrt(2/ro)*sqrt(x4) ;
epsx=epsx*lO;
end
%% RESULTING 4 POINTS FOR CUBIC SPLINES
X=[XliX2;x3;x4j
y=[yl:y2;y3;y4]
110
III
%%Generate H Matrix
for i=I:I:(nl);
k(i)=x(i+l)x(i) ;
end
h=k' ;
~~ Generate RHS Matrix
rhs (1,l) =0;
rhs(n,I)=O;
for i=2:1:(nl);
rhs (i, 1) =6* ( (y (i +1, 1) y (i, 1) ) Ih (i, 1)  (y (i, 1) y (i 1, 1) ) 1h (i 1, 1) ) ;
end
~% Generate A matrix
a(I,l)=I;
a(n,n)=I;
for i=2:1: (nl)
a (i, iI) =h (i 1, 1) ;
a (i, i +1) =h (i, 1) ;
a (i, i) =2 * (h (i 1, 1) +h (i, 1) ) ;
end
%% Solve System to get S Matrix
s=inv(a)*rhs;
!~ Use S values to ca culate coefficients of cubic splines
for i=I:1: (nl)
interval=i;
cubeco (i )= (s ( i +I, 1 )  s (i, 1) ) 1 (6*h ( i , 1 ) ) ;
squareco(i)=s(i,1)/2;
'nco(i)=((y(i+l, )y(i,1))/h(i,I))
(2 *h (i, 1) *s (i, 1) +h (i, 1) *s (i+1, 1) ) /6;
end
;~ Generate points for plotting
t~ Only the middle interval is important
incr=2.5; %pressure increment for plot points
i=l; ~~ generate plot points for cubic (xp,gx)
for xv=x(2,1) :ncr:x(3,1)
xp(i,I)=xv;
gx(i,I)=(( (xvx(2,1) )A3)*cubeco(2)+((xvx(2,1))A2)*squareco(2)+(xvy.
(2,1)) *linco (2) +y(2, 1));
i=i+l;
end
~~ Generate plot points for Turbulent flow Equation (xpl,gxl)
xv=eps:incr:90;
xpl=xv;
gxl=area* (cd*sqrt(2/ro)*(xv) .AO.5);
~ Plot the original test data
express=[O 2.5 5 7.5 10 12.5 15 17.5 20 22.5 25 30 35 40 45 50 55 60
65 70 75 80 85];
exq=cf*[O 1.8e3 3.5e3 4.98e3 6.35e3 7.48e3 8.7e3 9.71e3 1.07e2
1.14e2 1.23e2 1.40e2 1.54e2 1.67e2 1.78e2 1.93e2
2.08e2 2.18e2 2.32e2 2.43e2 2.52e2 2.5ge2 2.66e2];
plot (express,exq, 'b');
hold
plot (express, exq, '+') i
plot (xp,gx, 'r') %% plot cubic spline model
plot(xpl,gxl, 'g') %%plot turbulent flow eq
plot(x3,y3, '0') %% plot transition point
xlabel('Pressure (psid) ')
ylabel('Flow Rate (cu.in/sec) ')
%% print resulting cubic spline coefficients to screen
coeffs=[cubeco(2)isquareco(2)ilinco(2)1
112
Appendix D
MATLAB Scn.pt Files for Example 3.2
11
%File Name: staticrv3.m
%File Description:
% MATLAB Script file used to generate
~ the first of two cubic splines for
% the static relief valve model of
% Example 3.2
%%DECLARE EMPTY MATRICES FOR PLOTTING PURPOSES
xv= (]:
xp= [];
xpl=[];
xp2= [] ;
gx= (];
gxl=[] :
gx2=[] ;
%%INITIALIZE VARIABLE CONSTANTS
n=4; %nurnber of points for cubic spline
cd=.61; %orifice discharge coefficient
Pcr=800i %cracking pressure (psid)
Pmax=850; %Maximum pressure (psid)
Qmax=45.045; %Maximum flow rate (cu. in. !sec)
Kv=Qmax!sqrt(Pmax); %lumped constant for turbo flow EQ.
slope=Qmax!(PmaxPcr) i %slope (m) of linear (spring) portion of model
offset=Qmaxslope*Pmaxi %offset (b) of straight line equation
%% GENERATE CONTROL POINT AT 1st End point
%% (xl,y1) True end point
%% (x2,y2) = Control point
x1=720; %Desired endpoint pressure
yl=O; %flow rate corresponding to xl
epsx=l;
x2=x1+epsx*eps;
y2=O;
while ( (x2xl)<eps)
x2=xl+epsx*eps;
epsx=epsx*IO:
end
~% GENERATE CONTROL POINT AT 2nd End point
~~ (xl,yl) True end point
%% (x2,y2) = Control point
x3=825; %Desired endpoint pressure
y3=slope*x3+offset: ~flow rate corresponding to x3
epsx=l;
x4=x3+epsx*epsi
y4=slope*x4+offset;
while( (y4y3)<eps)
x4=x3+epsx*epsi
y4=slope*x4+offset;
epsx=epsx*lO;
end
~% RESULTING 4 POINTS FOR CUBIC SPLINES
x=[xI:x2;x3;x4j
y=[yl;y2;y3;y4]
114
%%Generate H Matrix
for i=l: 1: (n1) ;
k ( i )=x (i +1 )  x (i) ;
end
h=k' ;
%% Generate RHS Matrix
rhs(l,I)=O;
rhs(n,ll=O;
for i=2:1:(n1);
rhs (i, 1) =6* ( (y (i +1, 1)  Y ( i, 1) ) Ih (i, 1)  (y (i, 1) y (i 1, 1) ) /h (i 1, 1) ) ;
end
%% Generate A matrix
a(l,l)=l;
a(n,n)=l;
for i=2:1: (nl)
a(i,i1)=h(i1,1) ;
a (i, i +1) =h (i, 1) ;
a(i,i)=2*(h(i1,1)+h(i,1));
end
%% Solve System to get S Matrix
s=inv(a)*rhs;
%% Use S values to calculate coefficients of cubic splines
for i=l:l: (nl)
interval=i;
cubeco (i) = (s (i +1, 1)  s (i, 1) ) I (6*h (i , 1) ) ;
squareco(i)=s(i,1)/2i
1i nco (i )= ( (y (i +1, 1 )  Y (i, 1) ) I h (i, 1) ) (
2*h (i, 1) *s (i, 1) +h (i, 1) *s (i +1, 1) ) /6;
end
II
~% Generate points for plotting
.~:~ Only the middle interval is important
incr=(x(n,l)x(l,l) )/100; %pressure increment for plot points
i=l; %% generate plot points for cubic (xp,gx)
for xv=x(2,1) :incr:x(3, 1)
xp (i, 1) =x v ;
gx(i,1)=(((xvx(2,1) )A3)*cubeco(2)+((xvx(2,1))~2)*squareco(2)+(xvx(
2,1) )*linco(2)+y(2,1));
i=i+1;
end
~% Generate plot points for straight line
~% portion of original model
xv=O:incr:Pmax+incri
xp1=xv;
gxl=slope*xv+offset;
%% Generate plot points for turb flow Eq.
%% portion of original model
xv=O:incr:1.4*Pmax;
xp2=xv;
gx2=Kv*xp2.~.5;
plot (xp,gx, 'r') %% plot cubic spline
axis([O 1.4*Pmax 0 1.4*Qmax))
hold
plot(xp2,gx2, 'b')%% plot