An Analytical Solution for Unsteady Laminar Flow in Tubes with a Tapered Wall Thickness

Abstract

Transient fluid flows through tubes are critical in such topics as water hammer, ram pumps and pipeline dynamics. While analytical solutions exist in the literature for simple geometries such as tapered and non-tapered tube diameters, one area that is lacking is the case where the wave speed changes along the length. An example of this is a flexible pipe with a tapered wall thickness. In order to calculate the transient pressure response of such a system, this previously required a computationally expensive gridded method of characteristics (MOC) solution. This paper describes an analytical solution to the dynamic laminar flow of liquid in a tube where the wave speed varies along its length. This frequency-domain solution includes frequency-dependent friction effects. A comparison to a method of characteristics (MOC) solution is used to verify the solution. The paper also discusses some numerical issues and provides an approximate method that can be used for high-frequency calculations where limited numerical precision can cause errors. Finally, a preliminary comparison of the computational performance is presented, in which the new method is an order of magnitude faster to calculate than an MOC solution.

1. Introduction

This paper describes an analytical solution for pressure and flow dynamics for laminar liquid flow within a circular tube with constant inner diameter and varying wave propagation speed along its length, as shown in Figure 1. Fluid wave speed inside a non-rigid tube is dependent on the mechanical properties and the geometry of the tube as well as the fluid properties. Considering the pressure wave propagation along a tube, mechanical defects such as erosion, corrosion, and leakage have significant effects on the dynamic response of the pressure and flow. The motivation for this research is the use of such dynamic effects to allow for estimation of the remaining wall thickness in an eroded pipe at a lower cost than existing inspection techniques such as ultrasonic inspection, fiber-optic monitoring or “smart pigs” [1,2,3]. Other applications of this model would include research into water hammer [4], ram pumps [5,6,7], switched inertance hydraulic converters [8,9,10,11] and transient flow measurement [12].

The frequency-domain solution method presented here is computationally faster than existing gridded methods (such as the method of characteristics [13]), as it only calculates pressure and/or flow at the end-points of the tube. It can also be expected to be more accurate than time-domain solutions, such as the transmission line method (TLM), that require an approximation for the frequency-dependent friction terms [13,14,15]. While frequency-domain analytical solutions exist for tubes with tapered and non-tapered geometries [16,17,18], the authors are not aware of any such solution for the case where the wave speed varies.

This paper is organized with the preceding introduction, a development of the frequency-domain model and some sample results comparing the model to a conventional method, followed by conclusions.

2. Materials and Methods

Consider the laminar flow of a liquid through a tube, as shown in Figure 1. If the tube walls are not perfectly rigid, the wave speed is affected by the compliance of both the fluid and the containing tube walls. For example, if the tube is axially constrained, with expansion joints throughout, the elastic wave speed is given by

$$c=\sqrt{\frac{\left(\frac{K}{\rho }\right)}{1+\frac{2Kr}{Ee}}}$$

(1)

where $K$ is the fluid bulk modulus (usually the adiabatic tangent bulk modulus), $\rho$ is the fluid density, $r$ is the inner radius of the tube, $e$ is the wall thickness and $E$ is the wall material’s elastic modulus [4]. (Similar relationships are given in [4] for other anchoring conditions.) If the wall thickness changes along the length of the tube, then so too will the wave speed.

Zielke [19] gives the equations of continuity and motion for 1D laminar flow through a tube as follows (as restated using nomenclature from [20])

$$\frac{\partial q}{\partial t}+\frac{q}{A}\frac{\partial q}{\partial x}+\frac{A}{\rho }\frac{\partial p}{\partial x}+h\left(q\right)=0$$

(2)

$$\frac{\partial p}{\partial t}+\frac{q}{A}\frac{\partial p}{\partial x}+\frac{\rho c^2}{A}\frac{\partial q}{\partial x}=0$$

(3)

where $x$ is the position in the streamwise direction along the pipe, $p\left(t,x\right)$ is the pressure, $q\left(x,t\right)$ is the average volumetric flow rate, $A$ is the cross-sectional area and $h\left(q\right)$ is a friction function.

Following Johnston [20], if the fluid velocity is much less than the sonic velocity, then the $\frac{q}{A}$ terms can be neglected. Taking the Fourier transform then gives

$$j\omega Q+\frac{A}{ \rho }\frac{\partial P}{\partial x}+\frac{\nu }{r^2} W\left(j\omega \right)Q=0$$

(4)

$$j\omega P+\frac{\rho c^2}{A}\frac{\partial Q}{\partial x}=0$$

(5)

where $W\left(j\omega \right)$ is a friction function [20]. A common friction function [20,21] is

$$W\left(j\omega \right)=\frac{j\alpha }{\frac{z{J}_0\left(z\right)}{2J_1\left(z\right)}-1}$$

(6)

$$\alpha =\frac{\omega r^2}{\nu }$$

(7)

$$z=j\sqrt{j\alpha }$$

(8)

where $J$ is the Bessel function of the first kind. This function takes into account unsteady frequency-dependent friction, but not thermal effects.

We now diverge from Johnston [20] to assume that $A$ is constant and that $K$ (and therefore $c$) is a function of $x$. The fluid density, $\rho$, is still assumed to be constant. Solving Equation (5) for $P$ and differentiating gives

$$\frac{\partial P}{\partial x}=-\frac{\rho }{j\omega A}\left[c^2\frac{{\partial }^{2}Q}{\partial x^2}+\frac{\partial Q}{\partial x}\frac{\partial c^2}{\partial x}\right]$$

(9)

which can be substituted into Equation (4) to give the second order ODE

$$\frac{{\partial }^{2}Q}{\partial x^2}+\frac{1}{c^2}\frac{\partial c^2}{\partial x}\frac{\partial Q}{\partial x}-j\omega \left[\frac{\nu }{r^2} W\left(j\omega \right)+j\omega \right]\frac{Q}{c^2}=0$$

(10)

We now assume that the wave speed’s variation along the tube takes the form

$$c^2=c_A^2{e}^{Bx}$$

(11)

where

$$B=\frac{1}{L}\ln \left(\frac{c_B^2}{c_A^2}\right)$$

(12)

and where $c_A$ and $c_B$ are the wave speeds at the inlet ($x=0$) and outlet ($x=L$), respectively. As demonstrated in Figure 2, this nonlinear taper approximates a linear taper as long as the differences in $c$ are small relative to its absolute value. If this assumption is broken, the tube can be divided into a number of shorter segments.

For the nonlinear taper given above, the average velocity (as derived in Appendix A) is

$$\overline{c}=c_A\left[\frac{\ln \left(\frac{c_B}{c_A}\right)}{1-\frac{c_A}{c_B}}\right].$$

(13)

This taper has the nice property of

$$\frac{1}{c^2}\frac{\partial c^2}{\partial x}=B$$

(14)

which allows Equations (11) and (14) to be substituted into Equation (10) to give

$$\frac{{\partial }^{2}Q}{\partial x^2}+B\frac{\partial Q}{\partial x}-j\omega \left[\frac{\nu }{r^2} W\left(j\omega \right)+j\omega \right]\frac{Q}{c_A^2{e}^{Bx}}=0 .$$

(15)

If we define

$$C=\frac{j\omega }{c_A^2}\left[\frac{\nu }{r^2} W\left(j\omega \right)+j\omega \right]$$

(16)

then the ODE takes the form of

$$\frac{{\partial }^{2}Q}{\partial x^2}+B\frac{\partial Q}{\partial x}-C\frac{Q}{e^{Bx}}=0$$

(17)

which has an analytical solution:

$$Q\left(j\omega ,x\right)=e^{-\left(\frac{Bx}{2}\right)}[C_2{Y}_1\left(\gamma \right)+C_1{J}_1\left(\gamma \right)]$$

(18)

$$\gamma \left(x\right)=-\frac{2}{B}\sqrt{-C}{e}^{-\frac{Bx}{2}}$$

(19)

where $J$ and $Y$ are the Bessel functions of the first and second kind, and $C_1$ and $C_2$ are constants of integration.

This can be substituted into Equation (5) to obtain the pressure solution

$$P\left(\omega ,x\right)=-\frac{c_A^2\rho }{Aj\omega }\sqrt{-C}\left[C_2{Y}_0\left(\gamma \right)+C_1{J}_0\left(\gamma \right)\right]$$

(20)

We can then apply open or closed boundary conditions to obtain the transmission matrix in the form

$$\left[\begin{array}{c}{P}_A\left(j\omega \right)\\ Q_A\left(j\omega \right)Z_{cA}\end{array}\right]=\left[\begin{array}{cc}{t}_{11}\left(j\omega \right)& t_{12}\left(j\omega \right)\\ t_{21}\left(j\omega \right)& t_{22}\left(j\omega \right)\end{array}\right]\left[\begin{array}{c}{P}_B\left(j\omega \right)\\ Q_B\left(j\omega \right)Z_{cB}\end{array}\right]$$

(21)

where $P_A$ and $P_B$ are the pressures at the inlet and outlet, $Q_A$ and $Q_B$ are the corresponding flow rates, and $Z_{cA}$ and $Z_{cB}$ are the characteristic impedances at the inlet and outlet:

$$Z_{cA}=\frac{\rho c_A}{\pi r^2}$$

(22)

$$Z_{cB}=\frac{\rho c_B}{\pi r^2}.$$

(23)

The transmission matrix terms are

$$t_{11}=\frac{J_0\left({\theta }_1\right)Y_1\left({\theta }_2\right)-J_1\left({\theta }_2\right)Y_0\left({\theta }_1\right)}{D}$$

(24)

$$t_{12}=\left(\frac{c_A\sqrt{-C}}{j\omega }\right)\frac{J_0\left({\theta }_1\right)Y_0\left({\theta }_2\right)-J_0\left({\theta }_2\right)Y_0\left({\theta }_1\right)}{D}$$

(25)

$$t_{21}=-\left(\frac{j\omega }{c_A\sqrt{-C}}\right)\frac{J_1\left({\theta }_1\right)Y_1\left({\theta }_2\right)-J_1\left({\theta }_2\right)Y_1\left({\theta }_1\right)}{D}$$

(26)

$$t_{22}=-\frac{J_1\left({\theta }_1\right)Y_0\left({\theta }_2\right)-J_0\left({\theta }_2\right)Y_1\left({\theta }_1\right)}{D}$$

(27)

where

$${\theta }_1=-\frac{2}{B}\sqrt{-C}$$

(28)

$${\theta }_2=\frac{c_A}{c_B}{\theta }_1.$$

(29)

By applying the Wronskian, the denominator can be simplified to

$$D=J_0\left({\theta }_2\right)Y_1\left({\theta }_2\right)-J_1\left({\theta }_2\right)Y_0\left({\theta }_2\right)=-\frac{2}{\pi {\theta }_2} .$$

(30)

The above equations are undefined at $\omega =0$, but steady state physical requirements can be used to find the values. For continuity in flow between inlet and outlet, and zero pressure drop for zero flow:

$$t_{11}\left(\omega =0\right)=1$$

(31)

$$t_{21}\left(\omega =0\right)=0$$

(32)

$$t_{22}\left(\omega =0\right)=\frac{c_A}{c_B}$$

(33)

For the proper steady state pressure drop

$$R=\frac{8\nu \rho L}{\pi r^4}$$

(34)

$$t_{12}\left(\omega =0\right)=\frac{R}{Z{c}_B}=\frac{8\nu L}{r^2{c}_B}$$

(35)

where $R$ is the Hagen–Poiseuille laminar resistance [22]. Note the relationship between $t_{12}$ and the dissipation number $\beta =\frac{R}{8{\overline{Z}}_c}$ where ${\overline{Z}}_c$ is the lumped characteristic impedance, defined for a lossless line by

$${\overline{Z}}_c=\sqrt{\frac{\mathbb{L}}{\mathbb{C}}}$$

(36)

where $\mathbb{L}$ is the lumped inductance and $\mathbb{C}$ is the lumped capacitance [14,23]. For the tapered-wall tube defined in this paper, the dissipation number is

$$\beta =\frac{\nu L}{c_{A}r}\sqrt{\frac{1-{\left(\frac{c_A}{c_B}\right)}^2}{\ln \left({\left[\frac{c_B}{c_A}\right]}^2\right)}}.$$

(37)

Numerical issues exist for combinations of high frequency and small tapers. This is due to the transmission matrix numerators being a small difference of large numbers. This problem can be mitigated by using quadruple precision floating point representations. An alternative is to use high-frequency approximations of the Bessel functions [24]. In the limit as $\theta \to \infty$

$$J_v\left(\theta \right)\approx \sqrt{\frac{2}{\pi \theta }}\cos \left(\theta -\frac{1}{2v\pi }-\frac{1}{4}\right)$$

(38)

$$Y_v\left(\theta \right)\approx \sqrt{\frac{2}{\pi \theta }}\sin \left(\theta -\frac{1}{2v\pi }-\frac{1}{4}\right)$$

(39)

This allows for the following simplifications of the transmission matrix

$$t_{11}\left(\omega \to \infty \right)=\sqrt{\frac{c_A}{c_B}}\cos \left({\theta }_2\left[\frac{c_B}{c_A}-1\right]\right)$$

(40)

$$t_{12}\left(\omega \to \infty \right)=\left(\frac{c_A\sqrt{-C}}{j\omega }\right)\sqrt{\frac{c_A}{c_B}}\sin \left({\theta }_2\left[\frac{c_B}{c_A}-1\right]\right)$$

(41)

$$t_{21}\left(\omega \to \infty \right)=-\left(\frac{j\omega }{c_A\sqrt{-C}}\right)\sqrt{\frac{c_A}{c_B}}\sin \left({\theta }_2\left[\frac{c_B}{c_A}-1\right]\right)$$

(42)

$$t_{22}\left(\omega \to \infty \right)=\sqrt{\frac{c_A}{c_B}}\cos \left({\theta }_2\left[\frac{c_B}{c_A}-1\right]\right)$$

(43)

We found that the approximation error of these equations is exceeded by the precision error of the original equations (for double precision floating point numbers) when the imaginary part of ${\theta }_2$ is greater than about 15.

3. Results

In order to verify the equations, we compared the calculated responses to a method of characteristics (MOC) solution. The MOC solution divides the tube into N = 100 segments, each with a constant wave velocity, and uses an approximate frequency-dependent friction model from [13]. Since the analytical solution presented in our paper does not include these approximations, one should not view the MOC solution as the “true” solution as much as a sanity check to ensure that no blunders were committed in our development.

Frequency-domain solutions were calculated using the transmission matrix terms for frequencies scaled so that 100 time points would be calculated in the time required for the wave to transit the length of the tube, and the solution would span 8000 such cycles. These 800,000 frequency points are not strictly required, but the calculation is fast and it was desirable to ensure the entire response was captured for low dissipation numbers. Considerably fewer frequency points would be required for more damped responses. The inverse Fourier transform was used to calculate the time-domain response from the frequency-domain response.

Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 show solutions for four boundary condition situations, involving all four of the transmission matrix terms. Figure 3, Figure 4, Figure 5 and Figure 6 show the effect of varying the outlet wave speed while holding other parameters constant (which would result in $\beta =0.0127$ for $c_A=c_B)$. The solutions agree well, within the expected accuracy of the MOC solution. Figure 7, Figure 8, Figure 9 and Figure 10 show the same boundary conditions, but now holding $\frac{c_B}{c_A}=0.5$ constant and varying the viscosity to show a range of dissipation numbers (refer to Equation (37)). Again, the solutions agree well across the range of responses. At this point, we have not identified bounds for the applicability of our solution method, but we expect it should be suitable for all practical situations where the assumptions are satisfied.

In order to provide a sense of the relative computational cost for the new frequency domain method, responses were calculated for the case of $\beta ={10}^{-3}$ and $\frac{c_A}{c_B}=1$ while varying the number of frequency points in the solution. The method of characteristics solution was then calculated, with the grid selected to give results at the same time resolution. The computation time was recorded for each calculation, as implemented in MATLAB 2019b and run on a desktop computer with a i7-8700 processor with 32 GB RAM. In each case, the solution was calculated to an end time of four times the response’s settling time (defined as the point where transients reduce to 1% of the maximum value). The relative root mean squared error between the calculated response, $p\left(t\right),$ and the highest resolution frequency-domain method response calculated, $\widehat{p}\left(t\right)$, normalized by the standard deviation, was calculated as

$$E=\frac{\sqrt{\overline{{\left(p-\widehat{p}\right)}^2}}}{\sqrt{\overline{{\left(\overline{\widehat{p}}-\widehat{p}\right)}^2}}}$$

(44)

where an overbar designates the mean [25].

Figure 11 and Figure 12 shows the resultant calculation time and error as more terms are used in the solution. As expected, as the resolution is increased, the computation time and accuracy increase for both methods. However, for a given resolution, the new method is both faster (by around an order of magnitude) and more accurate. While the presented results are for a specific problem, we observed similar results across a wide range of dissipation numbers and taper ratios. It should be noted that neither the MOC solution nor the new frequency domain solution code was extensively optimized for speed, so these results should be viewed as a preliminary indication of relative computational speed.

One downside of the new method is that sufficient frequency points must be selected so that the transients have died out by the end of the corresponding time-domain response. Otherwise, the periodic nature of the Fourier transform will cause the response beyond the calculated end time to be added to the start of the response. If only the very beginning of a long response is of interest, the MOC method can be faster, as the computation can be halted early, while the new frequency-domain method must be fully calculated.

4. Summary and Conclusions

This paper presents an analytical solution for pressure and flow dynamics for laminar liquid flow within a tube, where the wave speed varies along the tube length in response to a change in compliance. This describes the situation of a tube with a flexible pipe wall where the wall thickness varies along the tube length but may also be applicable in other situations where the wave speed varies.

We have identified a numerical issue with standard double-precision floating point calculations and recommend two methods of mitigating these effects.

This paper also presents a comparison of the new method to a method of characteristics solution which agree well over a range of taper ratios and dissipation numbers, leading to confidence in the solution. A preliminary comparison of the calculation time and accuracy as the number of terms in the solution are varied is presented, which indicates that the new method is at least as accurate and much faster than the MOC solution.

Future work includes a more rigorous evaluation of the range of applicability of this new solution method, in terms of dissipation number and taper ratio, as well as an analysis of the number of frequency terms required for acceptable accuracy.