Analytical Modeling and Simulation of Machinery Containing Hydraulic Lines with Fluid Transients

1. Introduction

Hydraulic power transmission lines are common in industrial equipment. Consider the schematic of a fluid line shown in Figure 1.

The boundary conditions for the line in Figure 1 are defined by the pressure variables $\left(P_a,P_b\right)$ and flow rate variables $\left(Q_a,Q_b\right)$ shown on the schematic. Analyses of the fluid transients in a line often focus on variables such as shear stress and pressure at different locations in the line. In this paper, however, the focus is on the effects of the fluid transients interacting with components such as pumps, valves, and actuators at the ends of the line. Such components define the boundary conditions at the ends of a line. By analytically solving the partial differential equations in terms of yet to be specified boundary conditions, it is possible to obtain transfer function equations applicable to any set of boundary conditions that meet the requirements of causality [1].

A survey paper [1] on modeling techniques for fluid transients in smooth circular lines (Figure 1) presents frequency domain analytical solutions for a variety of assumed complexities in the model for the line and fluid properties. Equation (1) is a matrix representation of the equations from this reference for the analytical solution of the ‘dissipative model’ which includes viscosity and heat transfer effects in the model for fluid transients in lines with Newtonian fluids.

$$\left[\begin{array}{c}{P}_a(s)\\ Q_a(s)\end{array}\right]=\left[\begin{array}{cc}cosh\mathsf{\Gamma }(s)& Zsinh\mathsf{\Gamma }(s)\\ {\displaystyle \frac{sinh\mathsf{\Gamma }(s)}{Z}}& cosh\mathsf{\Gamma }(s)\end{array}\right]\left[\begin{array}{c}{P}_b(s)\\ Q_b(s)\end{array}\right]$$

(1)

In Eq. (1), ‘s’ is the Laplace transform variable, $\mathsf{\Gamma }(s)$ is the propagation operator, and Z(s) is the characteristic impedance. Equation (1) has been proven to be experimentally accurate for laminar mean flow conditions and “over a limited turbulent mean flow range”. The hyperbolic Bessel functions in Eq. (1) can be expressed in terms of infinite order product series for identifying dominant modes, resonant frequencies, and damping ratios [2].

Previous publications document reformulations of Eq. (1) for turbulent flow [3,4] and for a certain class of non-Newtonian fluids [5,6,7,8,9]. A review of the derivation of these reformulations is presented in Appendix A. In all cases, the reformulations fit the normalized small perturbation format shown in Eq. (2).

$$\left[\begin{array}{c}{\displaystyle \frac{∆P_b(\overline{s})}{P_e}}\\ {\displaystyle \frac{∆Q_b(\overline{s})}{Q_e}}\end{array}\right]={\left[\begin{array}{cc}cosh\Gamma +{\displaystyle \frac{R_T}{nZ}}sinh\Gamma & -{\displaystyle \frac{sinh\Gamma +\frac{R_T}{nZ}cosh\Gamma }{(8D_n\sqrt{1-B}+\frac{R_T}{Z})}}\\ -(8D_n\sqrt{1-B}+{\displaystyle \frac{R_T}{Z}})sinh\Gamma & cosh\Gamma \end{array}\right]}^n\left[\begin{array}{c}{\displaystyle \frac{∆P_a(\overline{s})}{P_e}}\\ {\displaystyle \frac{∆Q_a(\overline{s})}{Q_e}}\end{array}\right]$$

(2)

where

$$Z(\overline{s})={\displaystyle \frac{{\rho }_{e}a}{\pi r_o^2\sqrt{1-B}}}$$

(3)

$$\Gamma (\overline{s})={\displaystyle \frac{D_n\overline{s}}{\sqrt{1-B}}}$$

(4)

$$B={\displaystyle \frac{2J_1(jy)}{j\mathrm{y}{J}_0(j\mathrm{y})}}$$

(5)

$$D_n={\displaystyle \frac{{\mu }_{e}L}{{\rho }_{e}a{r}_o^2}}$$

(6)

$$y=\sqrt{{\displaystyle \frac{\overline{s}}{h\left(\overline{s}\right)}}}$$

(7)

$$h\left(\overline{s}\right)=\left[{\displaystyle \frac{1+{\overline{\lambda }}_2\overline{s}}{1+{\overline{\lambda }}_1\overline{s}}}\right]$$

(8)

where ${\rho }_e$ is the fluid density, $a$ is the wave speed, $r_o$ is the internal line radius, and ${\mu }_e$ is the pre-transient absolute viscosity. The perturbation variables in Eq. (2) are associated with changes (transients) relative to pre-transient steady-flow conditions. In Eq. (2), the pressure perturbations have been normalized by dividing by the pre-transient steady-flow pressure differential between the ends of the line, $P_e$. The flow rate perturbations have been normalized by dividing by the pre-transient steady-flow rate through the line, $Q_e$. The normalized Laplace variable, $\overline{s}=r_o^{2}s/v_e$, in Eq. (10) indicates that time has been normalized; $v_e$ is the pre-transient kinematic viscosity.

In Eq. (2), for pre-transient steady-flow Reynolds numbers < 1187.6, the flow is assumed laminar and n=1 and $R_T=0.$For Reynolds numbers > 1187.6, the pre-transient flow is assumed turbulent; n is an integer between 1 and 10 (based on convergence at relatively high Reynolds numbers) and $R_T$ is based on an empirical friction factor formula for smooth pipe (see Appendix A), i.e.

$${\displaystyle \frac{R_T}{Z}}=D_n\sqrt{1-B}\left(0.039544R_n^{0.75}-8\right)R_n>1187.6$$

(9)

It is important to note that the transition Reynolds number 1187.6 has been obtained by equating the laminar flow Darcy friction factor equation to the Blasius empirical formula for turbulent flow. Equation (8) is the Oldroyd-B shear stress transfer function for a class of non-Newtonian fluids. For ${\overline{\lambda }}_1>{\overline{\lambda }}_2$, viscosity decreases (shear thinning) as transient frequencies increase. For ${\overline{\lambda }}_1<{\overline{\lambda }}_2$, an increase in viscosity occurs as the frequency increases. For ${\overline{\lambda }}_1{=\overline{\lambda }}_2$, the viscosity is constant corresponding to a Newtonian fluid.

Section 2 of this paper contains examples using Eq. (2) to obtain total system transfer function models when components such as pumps, valves, and actuators are attached to the ends of the line. Since Eq. (2) is applicable to laminar or turbulent pre-transient flows with Newtonian or a certain class of non-Newtonian fluids, it is convenient to reformulate and express Eq. (2) in an all-encompassing simplified format to use for each of the examples, i.e.

$$\left[\begin{array}{c}\overline{\mathsf{\Delta }{P}_b}(\overline{s})\\ \overline{\mathsf{\Delta }{Q}_b}(\overline{s})\end{array}\right]=\left[\begin{array}{cc}{a}_{11}(\overline{s})& a_{12}(\overline{s})\\ a_{21}(\overline{s})& a_{22}(\overline{s})\end{array}\right]\left[\begin{array}{c}\overline{\mathsf{\Delta }{P}_a}(\overline{s})\\ \overline{\mathsf{\Delta }{Q}_a}(\overline{s})\end{array}\right]$$

(10)

It is important to note that for values of n > 1, the $a_{ij}$ terms in Eq. (10) are defined only after raising the matrix in Eq. (2) to the nth power. Obtaining transfer functions for a system with components such as pumps, valves, and actuators interconnected by fluid lines is the focus of Section 2.

Using these infinite order transfer functions in Eq. (10) to generate time domain simulations and important parameters such as resonant frequencies, time constants, and damping ratios is the focus of Section 3. Section 3 of the paper explains how to accurately convert total system transfer functions containing the $a_{ij}$ terms in Eq. (10) to accurate linear finite order equivalent transfer functions containing the dominant modes and the correct DC gain. Having these equivalent transfer functions enables one to perform time domain simulations and analyses without having to decide which terms in infinite series formulations to include, without having to figure out how to revise the series to achieve the correct DC gain [2], without having to deal with the complexities of an inverse Laplace transform [10], and without having to redo each method-of-characteristics numerical solution for each change in boundary conditions [11] associated with components at the ends of the line. A MATLAB mfile is included in Appendix B for implementing the inverse frequency algorithm and generating time domain simulations used in examples in Section 2.

2. Application Examples

There are 4-variables in the 2-equations in Eq. (10): $\overline{\mathsf{\Delta }{P}_a},$ $\overline{\mathsf{\Delta }{P}_b},$ $\overline{\mathsf{\Delta }{Q}_a},$ and $\overline{\mathsf{\Delta }{Q}_b}$. Thus, additional equations, considering the requirements of causality, associated with the boundary conditions must be specified to complete the analytical solution for a specified total system. The process of adding additional equations associated with the boundary conditions associated with component hardware is demonstrated in the examples that follow.

2.1. Water Hammer

This example is likely a case where the fluid transients have a detrimental effect on system components. Suppose the pressure at the left end of the line is constant, then $\overline{\mathsf{\Delta }{P}_a}(\overline{s})=0$. In this case, solving the two equations in Eq. (10) gives

$$\overline{\mathsf{\Delta }{P}_b}\left(\overline{s}\right)={\displaystyle \frac{a_{12}(\overline{s})}{a_{22}(\overline{s})}}\overline{\mathsf{\Delta }{Q}_b}(\overline{s})$$

(11)

The transfer function ${\displaystyle \frac{a_{12}(\overline{s})}{a_{22}(\overline{s})}}$ is the frequency domain infinite order analytical solution typically used to model the fluid transients (water hammer) due to sudden changes in the flow $\overline{\mathsf{\Delta }{Q}_b}\left(t\right).$ The transfer function in Eq. (11) can be used to generate a time domain solution for $\overline{\mathsf{\Delta }{P}_b}(\overline{t})$ given the time history of $\overline{\mathsf{\Delta }{Q}_b}(\overline{t})$. In water hammer studies, a typical assumption is that $Q_b\left(t\right)$ instantly decreases from the pre-transient flow rate $Q_e$ to 0; that is, $\overline{\mathsf{\Delta }{Q}_b}\left(\overline{t}\right)$ is a negative step function from 1 to 0 [11].

For laminar or relatively low levels of turbulence, n = 1 can be assumed [3,4]; in this case

$${\displaystyle \frac{a_{12}(\overline{s})}{a_{22}(\overline{s})}}=-{\displaystyle \frac{sinh\Gamma +\frac{R_T}{nZ}cosh\Gamma }{(8D_n\sqrt{1-B}+\frac{R_T}{Z})cosh\Gamma }}$$

(12)

When working with an infinite order transfer function, it is convenient to first get a finite order approximation for the transfer function. This can be accomplished using product series for hyperbolic and Bessel functions [1,2]; this approach has been proven to be most difficult when it comes to satisfying the low frequency gain associated with steady-state conditions. A more simple and direct approach for getting finite order approximations is by approximating the frequency response of the transfer function over a specified frequency range using an inverse frequency algorithm [3]. A frequency response pertains to the output if the input to the transfer function is a sinewave. Assuming a sinusoidal input to a transfer function, the output will become a sinusoidal function with the same input frequency; however, the output amplitude will either be greater or smaller than the input amplitude depending on the magnitude (dynamic gain) of the transfer function at the input frequency. The frequency response of a transfer function consists of two plots: (1) The magnitude plot representing the dynamic gain as a function of the input frequency, and (2) The phase angle plot resenting the time shift between the input sinewave and the output sinewave as a function of the input frequency. Transfer function approximations must match both parts of the frequency response over a specified range of frequencies. As will be demonstrated, an inverse frequency algorithm achieves this goal.

To generate the water hammer time domain solution for $\overline{\mathsf{\Delta }{P}_b}(\overline{t})$, the recommended first step is to get a frequency response plot of ${\left[{\displaystyle \frac{a_{12}(\overline{s})}{a_{22}(\overline{s})}}\right]}_{\overline{s}=j\overline{\omega }}$. The second step will be to use an inverse frequency algorithm to get a finite order approximation for the frequency response over a specified range of frequencies (see Section 3). An example of the results using the inverse frequency algorithm is shown in Figure 2, Figure 3 and Figure 4 for a water line 37.23 m in length, with a 22.1 mm diameter, and with a pre-transient Reynolds number of 5600. These line dimensions and pre-transient flow conditions are associated with experimental data [12] used in the verification of the analytical solutions for turbulent flow [3]. The MATLAB mfile for generating Figure 2, Figure 3 and Figure 4, that is also applicable to other transfer functions, is listed in Appendix B.

Referring to the magnitude plots in Figure 2, it is of interest to note that the input specification to the inverse frequency algorithm was to get a match over the normalized frequency range from 10^-11 to 8x10⁴; the minimum value 10^-11 was chosen to get the DC gain near the correct value of -1 and the maximum value 8x10⁴ chosen to accurately include at least seven 2nd order dominant modes. For purposes of resolution, the minimum normalized frequency shown on the graph is 500 even though 10^-11 was used in the algorithm. Also, note that the first three normalized resonant frequencies and peak values have been accurately identified with data tips. For example, at the normalized resonant frequency of 5689, the sinusoidal pressure amplitude amplification is ${10}^{{\displaystyle \frac{79.7}{20}}}=3.99.$ The order of the transfer function approximation needed to be 19 to get a good fit over the specified frequency range. The mfile generates eigenvalues, damping ratios, and time constants.

An accurate match of the phase plots over the frequency range is confirmed in Figure 3.

As shown in detail in Appendix B, the normalized water hammer pressure time history shown in Figure 4 was obtained by assuming a nearly instant valve closure input to the transfer function approximation for 5 different values of ‘n’. A comparison of the plots reveals that n = 1 would be sufficient for this low level of turbulence. As shown in Figure 4, the water hammer initial pressure surge following valve closure peaks at about 150$P_e$ = 4.155x10⁵ N/m².

2.2. Pump and Valve on the Ends of the Line

Consider the schematic in Figure 5 of a line with a pump and a valve on the ends of the line. For this example, the passage area through the valve may be a function of time. For the pump, the output flow or pressure may be a function of time.

Assume the flow through the valve is defined by the orifice equation with variable area $A_v$ and coefficient $C.$The orifice equation that models forward or reverse flow is

$${\mathrm{P}}_{\mathrm{b}}-{\mathrm{P}}_{\mathrm{d}}={\displaystyle \frac{{\mathrm{Q}}_{\mathrm{b}}\left|{\mathrm{Q}}_{\mathrm{b}}\right|}{{\mathrm{C}}^2{\mathrm{A}}_{\mathrm{v}}^2}}$$

(13)

Linearizing Eq. (13) about the pre-transient steady flow conditions gives

$${∆\mathrm{P}}_{\mathrm{b}}-{∆\mathrm{P}}_{\mathrm{d}}={\displaystyle \frac{2{\mathrm{Q}}_{\mathrm{e}}}{{\mathrm{C}}^2{\mathrm{A}}_{\mathrm{v}}^2}}{∆\mathrm{Q}}_{\mathrm{b}}-{\displaystyle \frac{2{\mathrm{Q}}_{\mathrm{e}}^2}{{\mathrm{C}}^2{\mathrm{A}}_{\mathrm{v}}^3}}∆{\mathrm{A}}_{\mathrm{v}}$$

(14)

To maintain consistency with the line transfer function equations, the pressure and flow variables are normalized by the equation for pre-transient steady flow through the line, i.e.

$$\overline{{∆\mathrm{P}}_{\mathrm{b}}}-\overline{{∆\mathrm{P}}_{\mathrm{d}}}=\mathrm{q}\left(\overline{{\mathsf{\Delta }\mathrm{Q}}_b}-\overline{{\mathsf{\Delta }\mathrm{A}}_v}\right)$$

(15)

The dimensionless parameter q represents the ratio of the pre-transient valve resistance to the line steady-flow resistance $q={\displaystyle \frac{2Q_e^2}{C^2{A}_{ve}^2{P}_e}}$.

The two equations in Eq. (10) and Eq. (15) contain 6-normalized perturbation variables: $\overline{{\mathsf{\Delta }\mathrm{P}}_b},$ $\overline{{\mathsf{\Delta }P}_a}$, $\overline{{\mathsf{\Delta }Q}_a}$, $\overline{{\mathsf{\Delta }Q}_b}$, $\overline{{\mathsf{\Delta }P}_d}$, and $\overline{{\mathsf{\Delta }A}_v}$, i.e.

$$\overline{{\mathsf{\Delta }\mathrm{P}}_b}\left(\overline{s}\right)=a_{11}\overline{{\mathsf{\Delta }P}_a}\left(\overline{s}\right)+a_{12}\overline{{\mathsf{\Delta }Q}_a}(\overline{s})$$

(16)

$$\overline{{\mathsf{\Delta }Q}_b}\left(\overline{s}\right)=a_{21}\overline{{\mathsf{\Delta }P}_a}\left(\overline{s}\right)+a_{22}\overline{{\mathsf{\Delta }Q}_a}\left(\overline{s}\right)$$

(17)

$$\overline{{\mathsf{\Delta }\mathrm{P}}_b}\left(\overline{s}\right)-\overline{{\mathsf{\Delta }P}_d}\left(\overline{s}\right)=q\left[\overline{{\mathsf{\Delta }Q}_b}\left(\overline{s}\right)-\overline{{\mathsf{\Delta }A}_v}\left(\overline{s}\right)\right]$$

(18)

Thus, to complete the solution, one or more of the variables must be designated as input(s) and/or additional equations added. Consider the examples in the following section.

2.2.1. Constant Flow from Pump with Changing Valve Area

For example, suppose $\overline{{\mathsf{\Delta }A}_v}(t)$ is designated to be an input, the pump has a fixed displacement with a constant output flow, $\overline{{\mathsf{\Delta }Q}_a}\left(t\right)=0,$ and the value output pressure is constant, $\overline{{\mathsf{\Delta }\mathrm{P}}_b}\left(t\right)=0.$ If the output of interest is $\overline{{\mathsf{\Delta }Q}_b}$, then Eq.’s (16-18) simplify to give the desired transfer function equation, i.e.

$$\overline{{\mathsf{\Delta }Q}_b}\left(\overline{s}\right)=\left[{\displaystyle \frac{1}{1-\frac{a_{11}}{q{a}_{21}}}}\right]\overline{{\mathsf{\Delta }A}_v}\left(\overline{s}\right)$$

(19)

Time domain plots of $\overline{{\mathsf{\Delta }Q}_b}\left(\overline{t}\right)$ for different valve area changes $\overline{{\mathsf{\Delta }A}_v}\left(\overline{t}\right)$ can be generated once the transfer function in Eq. (19) is approximated by a ratio of finite order rational polynomials.

2.2.2. Constant Valve Area with Changing Flow from a Variable Displacement Pump

In this example, fluid transients are desired and likely improve the performance of the system. Pumps are used to force fluids into downhole fissures in rock formations for oil and gas production. Pulsating flow rates from variable displacement pumps at resonant frequencies of the pressure transients in lines can enhance the performance of the system significantly. The transfer function for perturbations in the downhole flow rate can be obtained from Eq.’s (16-18). Assuming $\overline{{\mathsf{\Delta }Q}_a}(t)$ is specified to be an oscillating input and ${\overline{∆P}}_d(t)=0$ and $\overline{{\mathsf{\Delta }A}_v}\left(t\right)=0$, Eq.’s (16-18) simplify to

$$\overline{{\mathsf{\Delta }Q}_b}\left(\overline{s}\right)=\left[{\displaystyle \frac{a_{12}{a}_{21}-a_{11}{a}_{22}}{q{a}_{21}-a_{11}}}\right]\overline{{\mathsf{\Delta }Q}_a}\left(\overline{s}\right)$$

(20)

As with previous transfer functions derived in this section, the transfer function in Eq. (20) is infinite order. As shown in Figure 6, a frequency response of this transfer function with q = 10 reveals that the first mode resonant frequency occurs at a normalized frequency of 5689 with dynamic gain 21 dB. Thus, for a sinusoidal flow rate out of the pump at this resonant frequency superimposed on a mean flow, the amplitude of the sinusoidal output flow rate perturbations at the end of the line will be ${10}^{{\displaystyle \frac{21}{20}}}=11.22$ times larger than the input amplitude. The results and accuracy of the inverse frequency algorithm are also shown in Figure 6; having an approximation for the transfer function enables one to generate time domain plots of $\overline{{\mathsf{\Delta }Q}_b}\left(\overline{t}\right)$ for any other pump output flow rate specifications such as pulses, etc. It is important to note that the same mfile was used to generate Figure 6; the only changes in the mfile are the equation for the transfer function and the label on the y-axis of the graph.

2.3. Rotary Actuator Controlled by a Variable Displacement Pump

A variable displacement pump is used to power and control the speed of a hydraulic motor (rotary actuator) [13]. A schematic of the system is shown Figure 7.

The rotation speed, ${\dot{\theta }}_p,$of the pump input shaft is assumed constant. The rotation direction and speed of the motor is determined by the direction and magnitude of the flow in the lines which is indirectly determined by the pump displacement control lever. The significance of the fluid dynamics and flow losses in the lines are functions of the length and diameter of the lines and the frequencies of the input $\varnothing (t)$ and the load torque $T_L(t)$. The variables labeled on the schematic indicate the positive direction of the flows corresponding to counterclockwise rotation of the motor and clockwise torque representing the load on the motor shaft opposing the rotation of the motor. Most likely, it would be desirable for the input frequencies from the pump and/or load torque to be well below the resonant frequencies of the total system which are defined by the fluid transients coupled with the load dynamics. The following derivation formulates the transfer function for the total system which provides the eigenvalue frequencies that should be avoided.

The pump displacement $D_p$ is proportional to the swash plate angle $\varnothing$; thus, the flow from the pump in terms of the pump shaft angular velocity ${\dot{\theta }}_p$ is

$$Q_{ap}=D_p{\dot{\theta }}_p=C_p{\dot{\theta }}_p\varnothing$$

(21)

The flow through the motor is proportional to the motor displacement $D_m$, i.e.

$$Q_{am}=D_m{\dot{\theta }}_m$$

(22)

The torque equation on the motor shaft is

$$J{\ddot{\theta }}_m+b_f{\dot{\theta }}_m+C_m{D}_m\left(P_{am}-P_{bm}\right)+T_L=0$$

(23)

where $J$ is the rotational inertia of the load, $b_f$ is a viscous damping coefficient, and $C_m$ is the motor torque coefficient.

Laplace transforming Eq.’s (20-22) and then rewriting in small perturbation normalized format gives

$$\overline{∆Q_{ap}}(\overline{s})={\displaystyle \frac{C_p{\dot{\theta }}_p}{Q_e}}∆\varnothing (\overline{s})$$

(24)

$$\overline{∆Q_{am}}(\overline{s})={\displaystyle \frac{D_m{\nu }_e}{Q_e{r}_o^2}}\overline{s}∆{\theta }_m(\overline{s})$$

(25)

$$\left[\left({\displaystyle \frac{J{\nu }_e^2}{P_e{r}_o^4}}\right){\overline{s}}^2+\left({\displaystyle \frac{{\nu }_e{b}_f}{P_e{r}_o^2}}\right)\overline{s}\right]∆{\theta }_m\left(\overline{s}\right)+C_m{D}_m\left[\overline{∆P_{am}}(\overline{s})-\overline{∆P_{bm}}(\overline{s})\right]+{\displaystyle \frac{1}{P_e}}\mathsf{\Delta }{T}_L=0$$

(26)

Using Eq. (10) line ‘a’ and for line ‘b’ (see Figure 7), we get the normalized perturbation equations for the fluid transients in the lines, i.e.

$$\left[\begin{array}{c}\overline{∆P_{am}}(\overline{s})\\ \overline{∆Q_{am}}(\overline{s})\end{array}\right]=\left[\begin{array}{cc}{a}_{11}(\overline{s})& a_{12}(\overline{s})\\ a_{21}(\overline{s})& a_{22}(\overline{s})\end{array}\right]\left[\begin{array}{c}\overline{\mathsf{\Delta }{P}_{ap}}(\overline{s})\\ \overline{\mathsf{\Delta }{Q}_{ap}}(\overline{s})\end{array}\right]$$

(27)

$$\left[\begin{array}{c}\overline{∆P_{bp}}(\overline{s})\\ \overline{∆Q_{bp}}(\overline{s})\end{array}\right]=\left[\begin{array}{cc}{a}_{11}(\overline{s})& a_{12}(\overline{s})\\ a_{21}(\overline{s})& a_{22}(\overline{s})\end{array}\right]\left[\begin{array}{c}\overline{\mathsf{\Delta }{P}_{bm}}(\overline{s})\\ \overline{\mathsf{\Delta }{Q}_{bm}}(\overline{s})\end{array}\right]$$

(28)

Assuming no leakage, we can write

$${∆Q}_{ap}={∆Q}_{bp}$$

(29)

$${∆Q}_{am}={∆Q}_{bm}$$

(30)

In Eq.’s (24-30 ), there are two input variables, $∆\mathrm{\varnothing }$ and $∆T_L$ and there are nine unknowns: $\overline{∆Q_{ap}}(\overline{s})$, $\overline{∆Q_{am}}(\overline{s})$,$\overline{∆Q_{bp}}(\overline{s}),\overline{∆Q_{bm}}(\overline{s}),∆{\theta }_m(\overline{s})$, $\overline{∆P_{am}}(\overline{s})$, $\overline{∆P_{bm}}\left(\overline{s}\right),\overline{∆P_{ap}}(\overline{s})$, and $\overline{∆P_{bp}}\left(\overline{s}\right)$. Thus, it is possible to solve for the transfer functions for any of the unknown variables. For example, suppose that the output-of-interest is the perturbation in the motor speed $∆{\dot{\theta }}_m(t)$; this variable can be rewritten in normalized Laplace format, i.e.

$$L\left\{∆{\dot{\theta }}_m\left(t\right)\right\}=s∆{\theta }_m\left(s\right)={\displaystyle \frac{v_e}{r_o^2}}\overline{s}∆{\theta }_m\left(\overline{s}\right)$$

(31)

Eliminating all of the unknowns except $∆{\theta }_m\left(\overline{s}\right)$ from Eq.’s (24-30) gives

$${\displaystyle \frac{v_e}{r_o^2}}\overline{s}∆{\theta }_m\left(\overline{s}\right)={\displaystyle \frac{\left[\frac{C_p{C}_m{D}_m{\dot{\theta }}_p}{Q_e}\left(1+a_{11}{a}_{22}-a_{12}{a}_{21}\right)\right]\mathsf{\Delta }\varphi \left(\overline{s}\right)-\left[\frac{a_{21}}{P_e}\right]\mathsf{\Delta }{T}_L\left(\overline{s}\right)}{\frac{J{\nu }_e{a}_{11}}{P_e{r}_o^2}\overline{s}+\frac{a_{21}{b}_f}{P_e}+\frac{C_m{D}_m^2}{Q_e}\left(a_{11}+a_{22}\right)}}$$

(32)

There are two transfer functions in Eq. (32) for the inputs $∆\varnothing$ and $∆T_L$. Thus, the time history of the rotational speed of the motor, considering the fluid transients, can be obtained for any input specifications once the transfer functions have been approximated by ratios of finite order polynomials.

3. Suggestions for Approximating Infinite Order Transfer Functions Using an Inverse Frequency Algorithm

The MATLAB inverse frequency command is as follows:

[NumCoeffs,DenCoeffs]=invfreqs(FreqResp,Freq,NumOrder,DenOrder,Wt)

where ‘FreqResp’ are the complex frequency response values at the frequencies ‘Freq’, ‘NumOrder’ and ‘DenOrder’ are the desired orders of the polynomials in the transfer function approximation, and ‘Wt’ are optional weighting coefficients. A maximum frequency and a minimum frequency are used to generate the values ‘Freq’. I use 1000 points per decade and all weighting terms of one. Assigning values to the inputs in the ‘invfreqs’ function is a trial-and-error process. The following suggestions are recommended:

1. The first step is to estimate a minimum value and a maximum value that define a frequency range for a frequency response of the total system model such as in Eq.’s (11), (19), (20), and (32). Examination may reveal the presence of resonant frequencies and the low frequency requirements for an accurate DC gain. Change the minimum and maximum values until an acceptable frequency range is achieved. Consider the following statements used to generate the frequencies ‘Freq’ in your specified frequency range for the curve fit approximation.

wmin=.00001;wmax=120000; % Defines the frequency range for the curve fit

nd=ceil(log10(wmax)-log10(wmin)); %Determines the number of decades in the frequency range and rounds up to the next integer

w=logspace(log10(wmin),log10(wmax),1000*nd); %Generates at least 1000 points per decade to be used when getting the curve fit approximation

2. If resonant frequencies exist, then it is possible to estimate the order of the transfer function that will be required to accurately curve fit the response out to an upper frequency; at least 2 orders are required for each resonant peak and possibly additional orders for first order modes that are not apparent from examining the frequency response. To determine the number of dominant modes to include, keep increasing the order until there is little change in the eigenvalues and the time response.

dorder=23; % Desired order for the denominator of the transfer function

3. For a specific desired upper frequency, it is better to start with a lower order trial value and work up to a higher order with each trial to obtain an accurate fit out to that upper frequency

4. The order of the numerator should always be one less than the order of the denominator.

norder=dorder-1; % Desired order for the numerator of the transfer function

5. With each attempt, the DC gain of the transfer function should be checked. There will always be a trade off in accuracy at the lower and upper frequencies. The lower frequency should be decreased until the desired DC gain is achieved to the desired accuracy.

6. To generate time domain plots, it is important to note that the coefficients of the transfer function approximation are stored in memory. Simply create a transfer function with these coefficients and then use whatever linear simulation function is appropriate for a specific input function.

G=tf(NumCoeffs,DenCoeffs); % Creates the transfer function

Dcgain=dcgain(G) % Calculates the low frequency gain of the transfer function

step(G) % Generates a plot for a unit step input to the transfer function

The complete MATLAB mfile used to generate the approximation and time domain plots in Figure 2, Figure 3 and Figure 4 is listed in Appendix B.

4. Discussion

Accurate models for fluid transients in lines are defined by a set of simultaneous partial differential equations. As demonstrated in Appendix A, the analytical solution to these partial differential equations is obtained without having to specify specific details about components at the ends of the line. This is very desirable and provides significant versatilities to the analysis of the interaction of fluid transients with different end-of-line components. The analytical solution is applicable to either laminar or turbulent pre-transient flow for Newtonian [3] and a specific class of non-Newtonian fluids [9]. Compared to numerical solution approaches, such as the method-of-characteristics, the analytical solution to the partial differential equations does not have to be repeated with each end-of-the-line component design change and/or with each change in the external inputs.

Components at the ends of the line are typically modeled with ordinary differential equations and/or algebraic equations. These equations define the boundary conditions for a line which must satisfy causality requirements. Coupling the transfer functions for a line with equations for the components on the ends of the line to get a total system model, enables the generation of (1) time domain plots of system outputs-of-interest using linear or nonlinear simulation algorithms, (2) frequency response analyses, (3) computation of eigenvalue time constants, resonant frequencies, and damping ratios, and (4) models for controller designs for optimizing total system performance. Examples of coupling the analytical model for the line transients to equations for various fluid power components have been presented in Section 2. In some cases, the fluid transients are undesirable and should be avoided by operating the system at frequencies well below the system resonant frequencies. A waterhammer example was formulated and MATLAB code used to generate transfer function approximations and time domain simulations demonstrating the detrimental effects of fluid transients in lines. In a different example of a downhole fracking system, performance was enhanced by oscillating the flow rate input to the line at the first mode resonant frequency of the total system transfer function.

As demonstrated with the examples, once the boundary condition equations have been coupled with the Laplace transform solution to the fluid transient partial differential equations, a total system transfer function is achieved. These solutions have been derived using the Laplace transform which typically has not been a favorable approach due to the complexities of getting inverse Laplace transforms of infinite series formulations. In this paper, however, the complexities associated with the inverse Laplace transform are avoided by using an inverse frequency algorithm. Tips on using the inverse frequency algorithm are provided in Section 3; these tips are based on years of experience of trial-and-error selection of the frequency range and order of the polynomials in the transfer function approximation.

A MATLAB mfile is provided in Appendix B. The specific parameters and transfer function in the mfile pertain to the examples in Section 2.1 and Section 2.2.2. This mfile can be used with any transfer function formulated with the fluid transient equations and component equations with proper adjustment of the parameters as suggested in Section 3. Comment statements have been added to the MATLAB code to simplify the use of the mfile for other system transfer functions. The mfile contains many lines of code; however, it is important to note that the majority of the code pertains to the generation of output plots with appropriate labeling some of which may not be needed.

Appendix A.1. Newtonian Fluids

The closed form Newtonian analytical model for laminar flow transient perturbations in circular lines has been documented, experimentally verified, and used successfully for more than 50 years [1]. The primary objective of this Appendix is to extend this model for transient perturbations in non-Newtonian fluids and then extend it further for turbulent flow. The partial differential equations for the momentum, continuity, energy, and state for the dissipative model for laminar flow in a horizontal smooth line are listed in Eq.’s (A1-A3) [1]. The simplified isothermal version of these equations for compressible liquids with only $x$-direction axial flow are as follows

$$Momentum{\rho }_e{\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial p}{\partial x}}={\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial (r\sigma )}{\partial r}}$$

(A1)

$$Continuity{\displaystyle \frac{\partial \rho }{\partial t}}+{\rho }_e{\displaystyle \frac{\partial u}{\partial x}}=0$$

(A2)

$$State{\displaystyle \frac{d\rho }{{\rho }_e}}={\displaystyle \frac{dp}{{\beta }_e}}$$

(A3)

The axial velocity $u$ is a function of $r$ and $x$; $p$ is the pressure which is assumed to be uniform at any specific $x$; $\sigma$ is the stress tensor as a function of $r$ and $x$; ${\rho }_e$ is the pre-transient fluid density, and ${\beta }_e$ is the equivalent bulk modulus of the fluid which is a function of the elasticity of the line walls and can be adapted to include the effects of entrained gas.

For Newtonian fluids and laminar flow, the well-known frequency domain solution to Eq.’s (A1-A3) in terms of the normalized Laplace variable $\overline{\mathrm{s}}$is shown in Eq.’s (A4-A9) [1]. The details of this solution will be repeated in this appendix but with the addition of the non-Newtonian terms [4,5,6,7] and the impedance effects of turbulence [3]. This solution can be expressed in terms of normalized frequency using $\overline{\mathrm{s}}=j\overline{\omega }.$

$$\left[\begin{array}{c}{P}_b(\overline{s})\\ Q_a(\overline{s})\end{array}\right]=\left[\begin{array}{cc}{\displaystyle \frac{1}{cosh\Gamma (\overline{s})}}& -{\displaystyle \frac{Z(\overline{s})sinh\Gamma (\overline{s})}{cosh\Gamma (\overline{s})}}\\ {\displaystyle \frac{sinh\Gamma (\overline{s})}{Z(\overline{s})cosh\Gamma (\overline{s})}}& {\displaystyle \frac{1}{cosh\Gamma (\overline{s})}}\end{array}\right]\left[\begin{array}{c}{P}_a(\overline{s})\\ Q_b(\overline{s})\end{array}\right]$$

(A4)

$$D_n={\displaystyle \frac{{\mu }_{e}L}{{\rho }_{e}a{r}_o^2}}$$

(A5)

where

$$\Gamma ={\displaystyle \frac{D_n\overline{s}}{\sqrt{1-B}}}$$

(A6)

$$Z={\displaystyle \frac{{\rho }_{e}a}{\pi r_o^2\sqrt{1-B}}}$$

(A7)

$$B={\displaystyle \frac{2J_1(j\sqrt{\overline{s}})}{j\sqrt{\overline{s}}{J}_0(j\sqrt{\overline{s}})}}$$

(A8)

$$a={\displaystyle \frac{\sqrt{{\beta }_L/{\rho }_e}}{\sqrt{1+\frac{2r_o{\beta }_L}{\delta E}}}}$$

(A9)

The dimensionless parameter $D_n$is indicative of the attenuation and distortion along a line. It is the ratio of the wave travel time along a line to a term indicative of transient decay [1,14].

Appendix A.2. Extension of the Newtonian Analytical Model to Non-Newtonian Fluids

The well-known and widely accepted derivation for the Newtonian analytical model in Eq. (A4) is repeated here but, as mentioned above, expanded to include the Oldroyd-B stress tensor model [5,6,7,8,9]. Combining Eq.’s (A2) and (A3) gives

$${\displaystyle \frac{1}{{\beta }_e}}{\displaystyle \frac{\partial p}{\partial t}}+{\displaystyle \frac{\partial u}{\partial x}}=0$$

(A10)

Laplace transforming Eq.’s (A1) and (A10) assuming zero initial conditions and normalizing $s$ gives

$${\displaystyle \frac{{\rho }_e{v}_e\overline{s}}{r_o^2}}U\left(\overline{s}\right)+{\displaystyle \frac{dP(\overline{s})}{dx}}={\displaystyle \frac{1}{r}}{\displaystyle \frac{d(r\sigma \left(\overline{s}\right))}{dr}}$$

(A11)

$${\displaystyle \frac{v_e\overline{s}}{{r_o^2\beta }_e}}P\left(\overline{s}\right)+{\displaystyle \frac{dU(\overline{s})}{dx}}=0$$

(A12)

Substituting Eq. (3) for the Oldroyd-B shear stress model into Eq. (A11) gives

$${\displaystyle \frac{{\rho }_e{v}_e\overline{s}}{r_o^2}}U\left(\overline{s}\right)+{\displaystyle \frac{dP(\overline{s})}{dx}}={\mu }_{e}h(\overline{s})\left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{dU(\overline{s})}{dr}}+{\displaystyle \frac{d^{2}U(\overline{s})}{d{r}^2}}\right]$$

(A13)

Rearranging and collecting the terms in Eq. (A13) gives

$${\displaystyle \frac{d^{2}U(\overline{s})}{d{r}^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{dU(\overline{s})}{dr}}-{\displaystyle \frac{\overline{s}}{r_o^{2}h\left(\overline{s}\right)}}\left[U\left(\overline{s}\right)+{\displaystyle \frac{r_o^2}{\overline{s}{\rho }_e{v}_e}}{\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}\right]=0$$

(A14)

or

$${\displaystyle \frac{d^2\widehat{U}(\overline{s})}{d{r}^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d\widehat{U}(\overline{s})}{dr}}-c^2\widehat{U}(\overline{s})=0$$

(A15)

where

$$\widehat{U}\left(\overline{s}\right)=U\left(\overline{s}\right)+{\displaystyle \frac{r_o^2}{\overline{s}{\rho }_e{v}_e}}{\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}$$

(A16)

$$c^2={\displaystyle \frac{\overline{s}}{r_o^{2}h\left(\overline{s}\right)}}$$

(A17)

Eq. (A15) is a Bessel equation with solution

$$\widehat{U}\left(\overline{s}\right)=f\left(x,\overline{s}\right)J_0\left(jcr\right)=f(x,\overline{s})J_0(jy{\displaystyle \frac{r}{r_o}})$$

(A18)

where

$$y=\sqrt{{\displaystyle \frac{\overline{s}}{h\left(\overline{s}\right)}}}$$

(A19)

The function $f\left(x,\overline{s}\right)$ is not a function of $r$ and will be defined by the boundary conditions at the ends of the line. Using Eq.’s (A16) and (A18), the solution for $U\left(\overline{s}\right)$ can be written, i.e.

$$U\left(\overline{s}\right)=f\left(x,\overline{s}\right)J_0\left(jy{\displaystyle \frac{r}{r_o}}\right)-{\displaystyle \frac{r_o^2}{\overline{s}{\rho }_e{v}_e}}{\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}$$

(A20)

Considering that there is no slip at $r=r_o$, Eq. (A20) gives

$$0=f\left(x,\overline{s}\right)J_0\left(jy\right)-{\displaystyle \frac{r_o^2}{\overline{s}{\rho }_e{v}_e}}{\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}$$

(A21)

Thus,

$${\displaystyle \frac{r_o^2}{\overline{s}{\rho }_e{v}_e}}{\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}=f\left(x,\overline{s}\right)J_0\left(jy\right)$$

(A22)

Substituting Eq. (A22) into Eq. (A20) gives an equation for the velocity as a function of $x$, $r$, and $\overline{s}$, i.e.

$$U\left(\overline{s}\right)=f\left(x,\overline{s}\right)J_0\left(jy{\displaystyle \frac{r}{r_o}}\right)-f\left(x,\overline{s}\right)J_0\left(jy\right)$$

(A23)

At this point in the derivation, the assumption will be made that the flow passage is a circular pipe with radius $r_o$. Since the boundary conditions at the ends of the line will be in terms of the flow rates and pressures, it is convenient to compute the average velocity. Multiplying Eq. (A23) by $2\pi r$ and then integrating, and dividing by $\pi r_o^2$ gives

$$\overline{U}\left(\overline{s}\right)={\displaystyle \frac{1}{\pi r_o^2}}{\int }_0^{r_o}2\pi r\left[f\left(x,\overline{s}\right)J_0\left(jy{\displaystyle \frac{r}{r_o}}\right)-f\left(x,\overline{s}\right)J_0\left(jy\right)\right]dr$$

(A24)

$$\overline{U}\left(\overline{s}\right)=f\left(x,\overline{s}\right)M$$

(A25)

where

$$M={\displaystyle \frac{2}{jy}}{J}_1\left(jy\right)-J_0(jy)$$

(A26)

Solving for $P(\overline{s})$ in Eq. (A12) and then taking the derivative with respect to $x$ gives

$${\displaystyle \frac{dP(\overline{s})}{dx}}=-{\displaystyle \frac{{r_o^2\beta }_e}{v_e\overline{s}}}{\displaystyle \frac{d^{2}U(\overline{s})}{d{x}^2}}$$

(A27)

Taking the 2nd derivative of Eq. (A25) with respect to $x$ gives

$${\displaystyle \frac{d^{2}U(\overline{s})}{d{x}^2}}=M{\displaystyle \frac{d^{2}f\left(x,\overline{s}\right)}{d{x}^2}}$$

(A28)

From Eq. (A22)

$${\displaystyle \frac{dP\left(\overline{s}\right)}{dx}}={\displaystyle \frac{\overline{s}{\rho }_e{v}_e}{r_o^2}}f\left(x,\overline{s}\right)J_0\left(jy\right)$$

(A29)

Substituting Eq. (A27) into Eq. (A29) gives

$$-{\displaystyle \frac{{r_o^2\beta }_e}{v_e\overline{s}}}{\displaystyle \frac{d^{2}U\left(\overline{s}\right)}{d{x}^2}}={\displaystyle \frac{\overline{s}{\rho }_e{v}_e}{r_o^2}}f\left(x,\overline{s}\right)J_0\left(jy\right)$$

(A30)

Substituting Eq. (A28) into Eq. (A30) gives

$$-{\displaystyle \frac{{r_o^2\beta }_e}{v_e\overline{s}}}\left[M{\displaystyle \frac{d^{2}f\left(x,\overline{s}\right)}{d{x}^2}}\right]={\displaystyle \frac{\overline{s}{\rho }_e{v}_e}{r_o^2}}f\left(x,\overline{s}\right)J_0\left(jy\right)$$

(A31)

Rearranging Eq. (A31) gives the following equation

$${\displaystyle \frac{d^{2}f\left(x,\overline{s}\right)}{d{x}^2}}-A^{2}f\left(x,\overline{s}\right)=0$$

(A32)

where

$$A^2=-{\displaystyle \frac{v_e^2{\overline{s}}^2{\rho }_e{J}_0\left(jy\right)}{M{\beta }_e{r}_o^4}}={\displaystyle \frac{v_e^2{\overline{s}}^2{J}_0\left(jy\right)}{a^2{r}_o^4\left[J_0(jy)-\frac{2}{jy}{J}_1\left(jy\right)\right]}}$$

(A33)

The solution to Eq. (A32) is

$$f\left(x,\overline{s}\right)=C_1\cosh \left(Ax\right)+C_2\sinh \left(Ax\right)$$

(A34)

where $C_1$ and $C_2$ are determined by the boundary conditions at the ends of the line defined by the flow rates and pressures. Using Eq. (A25), we get

$$Q_a(\overline{s})=\pi r_o^2{\overline{U}}_{x=0}=\pi r_o^{2}Mf{\left(x,\overline{s}\right)}_{x=0}=\pi r_o^{2}M\left[C_1\cosh \left(0\right)+C_2\sinh \left(0\right)\right]$$

(A35)

$$Q_a(\overline{s})=\pi r_o^{2}M{C}_1$$

(A36)

$$Q_b(\overline{s})=\pi r_o^2{\overline{U}}_{x=L}=\pi r_o^{2}Mf{\left(x,\overline{s}\right)}_{x=L}=\pi r_o^{2}M\left[C_1\cosh \left(AL\right)+C_2\sinh \left(AL\right)\right]$$

(A37)

The equations for the pressures require the use of Eq. (A12) and the derivative of Eq. (A25), i.e.

$$P_a\left(\overline{s}\right)=-{\displaystyle \frac{{r_o^2\beta }_e}{v_e\overline{s}}}{\left[{\displaystyle \frac{dU\left(\overline{s}\right)}{dx}}\right]}_{x=0}=-{\displaystyle \frac{{r_o^2\beta }_{e}M}{v_e\overline{s}}}{\left[{\displaystyle \frac{df\left(x,\overline{s}\right)}{dx}}\right]}_{x=0}$$

(A38)

$$=-{\displaystyle \frac{{r_o^2\beta }_{e}M}{v_e\overline{s}}}{\left[C_{1}Asinh(Ax)+C_2\mathrm{Acosh}\left(Ax\right)\right]}_{x=0}$$

(A39)

Thus,

$$P_a\left(\overline{s}\right)=-{\displaystyle \frac{{r_o^2\beta }_{e}M}{v_e\overline{s}}}\left[C_{1}Asinh\left(0\right)+C_2\mathrm{Acosh}\left(0\right)\right]=-{\displaystyle \frac{{r_o^2\beta }_{e}M}{v_e\overline{s}}}{C}_{2}A$$

(A40)

$$P_b\left(\overline{s}\right)=-{\displaystyle \frac{{r_o^2\beta }_{e}M}{v_e\overline{s}}}\left[C_{1}Asinh(AL)+C_2\mathrm{Acosh}\left(AL\right)\right]$$

(A41)

From Eq. (A36)

$$C_1={\displaystyle \frac{Q_a\left(\overline{s}\right)}{\pi r_o^{2}M}}$$

(A42)

From Eq. (A40)

$$C_2=-{\displaystyle \frac{{v_e\overline{s}P}_a\left(\overline{s}\right)}{{r_o^2\beta }_{e}MA}}$$

(A43)

Thus, substituting Eq. (A42) and Eq. (A43) into Eq. (A37) and Eq. (A41) gives

$$\left[\begin{array}{c}{P}_b\left(\overline{s}\right)\\ Q_b(\overline{s})\end{array}\right]=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(AL)& -Z(\overline{s})sinh(AL)\\ -{\displaystyle \frac{1}{Z(\overline{s})}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(AL)& \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(AL)\end{array}\right]\left[\begin{array}{c}{P}_a\left(\overline{s}\right)\\ Q_a(\overline{s})\end{array}\right]$$

(A44)

or

$$\left[\begin{array}{c}{P}_b(\overline{s})\\ Q_a(\overline{s})\end{array}\right]=\left[\begin{array}{cc}{\displaystyle \frac{1}{cosh\Gamma (\overline{s})}}& -{\displaystyle \frac{Z(\overline{s})sinh\Gamma (\overline{s})}{cosh\Gamma (\overline{s})}}\\ {\displaystyle \frac{sinh\Gamma (\overline{s})}{Z(\overline{s})cosh\Gamma (\overline{s})}}& {\displaystyle \frac{1}{cosh\Gamma (\overline{s})}}\end{array}\right]\left[\begin{array}{c}{P}_a(\overline{s})\\ Q_b(\overline{s})\end{array}\right]$$

(A45)

where

$$Z(\overline{s})={\displaystyle \frac{{\rho }_{e}a}{\pi r_o^2\sqrt{1-B}}}$$

(A46)

$$\Gamma (\overline{s})={\displaystyle \frac{D_n\overline{s}}{\sqrt{1-B}}}$$

(A47)

$$B={\displaystyle \frac{2J_1(jy)}{j\mathrm{y}{J}_0(j\mathrm{y})}}$$

(A48)

Equation (A45) represents the non-Newtonian version of Eq. (A4). It is of interest to note that Eq. (A45) is identical to Eq. (A4) except for the inclusion of the Oldroyd-B shear stress terms in $y$ and in the parameter $B$ in Eq. (A48). Also, Eq.’s (A46) and (A47) are identical in form to Eq.’s (A7) and (A6), respectively. Setting ${\lambda }_1={\lambda }_2=0$ makes $h\left(\overline{s}\right)=1$; so, Eq. (A45) for non-Newtonian becomes identical to Eq. (A4) for Newtonian.

Appendix A.3. Extension to Pre-Transient Turbulent Flow

For Newtonian fluid, Eq. (A4) was extended to a model for pre-transient turbulent flow [3]. Since Eq. (A45) is identical to Eq. (A4) except for the inclusion of the non-Newtonian shear stress parameters in the parameters $B$ and $y$, the extension of Eq. (A4) to turbulent flow is applicable to the non-Newtonian model in Eq. (A45). The extension to pre-transient turbulent flow is summarized here in order to clarify the derivations of equations for pre-transient turbulent normalized shear stress and normalized velocity which are included in this Appendix.

The pre-transient turbulent flow model is based on using the analytical model in Eq. (A45) in series with a lumped resistance $R_T$ sized so that the resulting model for the line has the same steady flow resistance of a smooth line with pre-transient turbulent flow [3]. To potentially improve the accuracy of the model for relatively large Reynolds numbers, the resistance $R_T$ is distributed along the line by dividing the line into n-sections. The schematic for such a multi-sectioned model is shown in Figure A1 where $P_{ai},P_{bi,}{Q}_{ai},and{Q}_{bi}$ are intermediate pressures and flows.

Figure A1. Line of length L with a lumped resistance at the end of each segment of length L/n.

Figure A1. Line of length L with a lumped resistance at the end of each segment of length L/n.

The Blasius empirical model for the Darcy friction factor for smooth lines is utilized to establish the turbulent flow resistance $R_T$. For Reynolds numbers up to 200,000 [3,4] in smooth lines with steady flow $Q_e$, the corresponding pressure differential between the ends of a horizontal line is

$${P_e=(P}_a-{\mathrm{P}}_{\mathrm{b}}{)}_{\mathrm{e}}=\mathrm{b}{\mathrm{Q}}_{\mathrm{e}}^{1.75}$$

(A49)

where

$$b=\left[{\displaystyle \frac{0.008972{\rho }_e^{0.75}{\mu }_e^{0.25}L}{r_o^{4.75}}}\right]$$

(A50)

Denoting the steady flow laminar resistance by $R_L$, it can be shown that the additional resistance $R_T$ needed to match the total resistance for turbulent flow is given by

$${\displaystyle \frac{R_T}{Z}}={\displaystyle \frac{b{Q}_e^{0.75}-R_L}{Z}}=D_n\sqrt{1-B}(0.039544R_n^{0.75}-8)$$

(A51)

Using normalized perturbation variables representing pressure and flow changes from the pre-transient turbulent level, the resulting model for the transients in the line is [3].

$$\left[\begin{array}{c}{\displaystyle \frac{∆P_b(\overline{s})}{P_e}}\\ {\displaystyle \frac{∆Q_b(\overline{s})}{Q_e}}\end{array}\right]={\left[\begin{array}{cc}cosh\Gamma +{\displaystyle \frac{R_T}{nZ}}sinh\Gamma & -{\displaystyle \frac{sinh\Gamma +\frac{R_T}{nZ}cosh\Gamma }{(8D_n\sqrt{1-B}+\frac{R_T}{Z})}}\\ -(8D_n\sqrt{1-B}+{\displaystyle \frac{R_T}{Z}})sinh\Gamma & cosh\Gamma \end{array}\right]}^n\left[\begin{array}{c}{\displaystyle \frac{∆P_a(\overline{s})}{P_e}}\\ {\displaystyle \frac{∆Q_a(\overline{s})}{Q_e}}\end{array}\right]$$

(A52)

where

$$\Gamma ={\displaystyle \frac{D_n\overline{s}}{n\sqrt{1-B}}}$$

(A53)

and

$$P_e=(8D_n\sqrt{1-B}+{\displaystyle \frac{R_T}{Z}})Q_e$$

(A54)

It has been demonstrated that Eq. (A52) for Newtonian fluids gives good agreement with water hammer experimental data and MOC solutions with values of n from 1 to 10 and Reynolds numbers up to 15,800 [3]. It has been demonstrated for water hammer simulations, a value of 10 is sufficient for decay rate convergence for Reynolds numbers out to 200,000 [4].

Referring to Eq. (A52), it can be shown that $R_T=0$ if $R_n=1187.6$; this is the transition Reynolds number for which the Blasius and laminar flow friction factors are identical. It is also important to note that for $n$ =1 and R_T = 0, Eq. (A52) gives the solution for pre-transient laminar flow in a horizontal smooth line.