Dynamic stiffness method for free vibration analysis of variable diameter pipe conveying fluid

1. Introduction

Fluid-conveying pipe is widely used in nuclear industry, petroleum industry, and many other industries. The failure caused by the interaction between fluid and pipe wall occurs very frequently in practice. So flow-induced vibration of pipe conveying fluid attracts more and more interests of researchers. Ashley and Haviland [1] begun the research of bending vibrations of pipe line containing flowing fluid. Paidoussis [2] has studied the dynamic characteristics of pipe conveying fluid for many years, the great achievement he accomplished in nonlinear vibration of pipes set the foundation of pipe dynamics. Recently, flow-induced vibration of pipe conveying fluid is researched more exhaustively. Doare et al. [3] studied the dissipation effect on the local and global stability of pipe conveying fluid and proposed a numerical method to analyze the stability of finite-length system. Jung, Chung [4] investigated the stability of semi-circular pipe conveying harmonically oscillating fluid. They considered the extensibility and nonlinearity of the fluid-conveying pipe with different boundary conditions through Hamilton principle. Rinaldi, Prabhakar [5] studied micro-scale resonators containing internal flow, modeled as micro-fabricated pipe conveying fluid, and investigated the effects of flow velocity on damping, stability, and frequency shift. Huang, Liu [6] studied the natural frequencies of fluid conveying pipe with different boundary conditions through Garlerkin approach. Ragulskis, Fedaravicius [7] have developed a harmonic balance method in the FEM analysis of the fluid in pipe with taking the interactions between the vibrating pipe and fluid flow into consideration.

Uniform pipe was always adopted in the flow-induced vibration analysis of pipe conveying fluid. However, fluid-conveying pipes with variable section are widespread in nuclear engineering, such as pipe-expanding and nozzle etc. Hannoyer and Paidoussis [8, 9] have studied the dynamics of a slender tapered beam with internal and external flow in theory and experiment. Based on the researches of Hannoyer and Paidoussis [8, 9], we studied the dynamic characteristics of variable diameter pipes conveying fluid through stiffness matrix method. Natural frequencies and modal shapes of the variable diameter pipe were calculated in this paper and contrasts were made between the results in this paper and previous researches, which prove that the method employed in this paper is reasonable.

2. Motion equation of variable diameter pipe conveying fluid

Here, we use Hamilton’s principle to build the dynamic model of variable diameter pipe conveying fluid. Because that the fluid-conveying pipe is an open system. The Hamilton’s principle for open system is used here, i.e.:

where, $\mathcal{L}$ is the system’s Lagrange function, that is $\mathcal{L}=T_k-T_p$, where $T_k$, $T_p$ stand for the kinetic energy and potential energy of the system respectively.$\delta W$ is the work done by non-conservative force, $m_f$ is the fluid mass per unit length, ${\overrightarrow{r}}_e$ is the position vector of pipe exit. ${\overrightarrow{e}}_t$ is the unit vector in tangent direction of the deformed pipeline, as shown in Fig. 1.

In the coordinate shown in Fig. 1, the absolute fluid velocity can be written as follow:

In Eq. (2), $V\left(x\right)$ stands for the velocity of the fluid relative to the pipe. ${\overrightarrow{e}}_1$ and ${\overrightarrow{e}}_2$ are unit vectors in $x$- and $y$-directions.$\left(\right)\text{'}$ and $\dot{\left(\right)}$ stand for the partial derivatives with respect to $x$ (coordinate) and $t$ (time). According to the law of conservation of fluid mass, as shown in Eq. (3), the velocity of the fluid is change along the length of the conical pipe. So we use $V\left(x\right)$ and $V_a\left(x\right)$ here:

where $C$ is a constant.

Fig. 1The displacement sketch of simply supported pipe

Fig. 2The sketch of variable diameter pipe

The sketch of a variable diameter pipe is shown in Fig. 2. The fluid flows from left end to right end. The outer and inner diameters of left end of the pipe are $D_i$ and $d_{i}s$ respectively. The non-uniformity of the pipe is characterized by the diameter variation coefficient $\gamma$. The outer and inner diameters of the section at coordinate $x$ can be expressed as $D\left(x\right)=\left(1-\gamma \frac{x}{L}\right)D_i$ and $d\left(x\right)=\left(1-\gamma \frac{x}{L}\right)d_i$ respectively, where $L$ is length of the pipe. Because that the area of the fluid is $A_f\left(x\right)=\frac{1}{4}\pi d^2\left(x\right)$, the velocity of the fluid along the length of the pipe can be expressed as:

Given that the flow velocity at the entrance of the pipe is $V_0$, then Eq. (4) can be transformed into following form:

Now, we use Hamilton’s principle for open system to deduce the motion equation of variable diameter pipe conveying fluid. The kinetic energy of fluid in pipe can be written in the following form:

Here, the axial inextensible assumption was used, i.e. $\mathrm{}{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial w}{\partial x}\right)}^2=1$. And high order infinitesimals have been ignored in Eq. (6). ${\rho }_f$ is the density of fluid in pipe, $A_f\left(x\right)$ is area of the fluid section located at $x$. $m_f\left(x\right)dx={\rho }_f{A}_f\left(x\right)dx$ is mass of the fluid element analyzed.

The kinetic energy of the pipe is as follow:

where, ${\rho }_p$ is density of the pipe material, $A_p\left(x\right)$ is area of the pipe cross section at $x$. $m_p\left(x\right)dx={\rho }_p{A}_p\left(x\right)dx$ is mass of the pipe element.

The potential energy of the pipe is as follow:

In Eq. (8), $I_x$ is the inertial moment of the pipe section. Obviously, $I_x$ also changes along $x$axis. Only considering the work done by fluid pressure, the work done by non-conservative force is zero [2], that is:

Substituting Eq. (6)-(9) into Eq. (1), the following equation can be obtained:

Here, we take the simply supported pipe (as shown in Fig. 1) for example to deduce the lateral vibration equation of the variable diameter pipe conveying fluid.

The boundary conditions of simply supported pipe are:

At time $t_1$ and $t_2\text{:}$

Combine Eq. (11) with Eq. (12), the following equation is obtained:

The detail manipulation of the third term in Eq. (10) is as follow:

Considering the assumption of axial inextensible, we obtain that $u_e=-\frac{1}{2}\underset{0}{\overset{L}{\int }}{w^{\text{'}}}^{2}dx$, where the subscript ‘$e$’ means “exit”, and the following equation can be gained:

Now, we can get the motion equation of variable diameter pipe conveying fluid as follow:

3. Free vibration analysis of variable diameter pipe

3.1. Frobenius method

To harmonic vibration, we have:

17

$w\left(x,t\right)=W\left(x\right)e^{i\omega t}.$

Substituting Eq. (17) into Eq. (16), we can get the following ordinary differential equation:

18

$E{I}_x\frac{d^{4}W}{d{x}^4}+2E\frac{d{I}_x}{dx}\frac{d^{3}W}{d{x}^3}+\left(E\frac{d^2{I}_x}{d{x}^2}+m_f{V}^2\right)\frac{d^{2}W}{d{x}^2}+\left(m_{f}Vi\omega +m_f^{\text{'}{V}^2}+2m_{f}V{V}^{\text{'}}\right)\frac{dW}{dx}+\left[\left(m_f^{\text{'}}V+m_f{V}^{\text{'}}\right)i\omega -\left(m_p+m_f\right){\omega }^2\right]W=0,$

where:

19

$m_f\left(x\right)={\rho }_f{A}_f\left(x\right)={\rho }_f\pi \frac{{\left(1-\frac{\gamma x}{L}\right)}^2{d}_i^2}{4}=\left(1+{\alpha }_1\frac{x}{L}+{\alpha }_2\frac{x^2}{L^2}\right){\rho }_f{A}_f\left(0\right),$
$m_p\left(x\right)={\rho }_p{A}_p\left(x\right)={\rho }_f\pi \frac{{\left(1-\frac{\gamma x}{L}\right)}^2\left(D_i^2-d_i^2\right)}{4}=\left(1+{\alpha }_1\frac{x}{L}+{\alpha }_2\frac{x^2}{L^2}\right){\rho }_p{A}_p\left(0\right),$

and

20

$I_x=\pi \frac{D_x^4-d_x^4}{64}=\pi \frac{{\left(1-\frac{\gamma x}{L}\right)}^4\left(D_i^4-d_i^4\right)}{64}=\left(1+{\beta }_1\frac{x}{L}+{\beta }_2\frac{x^2}{L^2}+{\beta }_3\frac{x^3}{L^3}+{\beta }_4\frac{x^4}{L^4}\right)I_0.$

For the purpose of analyzing Eq. (18) conveniently, the following dimensionless parameters are introduced here:

21

$\xi =\frac{x}{L},\eta =\frac{W}{L},\mathrm{\Omega }=\sqrt{\frac{{\rho }_f{A}_f\left(0\right)+{\rho }_p{A}_p\left(0\right)}{EI\left(0\right)}}\omega L^2,$
$v_0=\sqrt{\frac{{\rho }_f{A}_f\left(0\right)}{EI\left(0\right)}}L{V}_0,\beta =\frac{{\rho }_f{A}_f\left(0\right)}{{\rho }_f{A}_f\left(0\right)+{\rho }_p{A}_p\left(0\right)}.$

The dimensionless form of Eq. (18) is as follow:

22

$\left(1+{\beta }_1\xi +{\beta }_2{\xi }^2+{\beta }_3{\xi }^3+{\beta }_4{\xi }^4\right)\frac{d^4\eta }{d{\xi }^4}+2\left({\beta }_1+2{\beta }_2\xi +3{\beta }_3{\xi }^2+4{\beta }_4{\xi }^3\right)\frac{d^3\eta }{d{\xi }^3}+\left[\left(2{\beta }_2+6{\beta }_3\xi +12{\beta }_4{\xi }^2\right)+\left(1+{\alpha }_1\xi +{\alpha }_2{\xi }^2\right)\frac{v_0^2}{{\left(1-\gamma \xi \right)}^4}\right]\frac{d^2\eta }{d{\xi }^2}+\left[\left({\alpha }_1+2{\alpha }_2\xi \right)\frac{v_0^2}{{\left(1-\gamma \xi \right)}^4}+\left(1+{\alpha }_1\xi +{\alpha }_2{\xi }^2\right)\frac{\gamma v_0^2}{{\left(1-\gamma \xi \right)}^5}+i\mathrm{\Omega }{v}_0\frac{1}{4{\left(1-\gamma \xi \right)}^2}\sqrt{\beta }\right]\frac{d\eta }{d\xi }+\left\{i\mathrm{\Omega }\sqrt{\beta }\left[\left({\alpha }_1+2{\alpha }_2\xi \right)\frac{1}{{\left(1-\gamma \xi \right)}^2}+\left(1+{\alpha }_1\xi +{\alpha }_2{\xi }^2\right)\frac{\gamma v_0}{{\left(1-\gamma \xi \right)}^3}\right]-\left(1+{\alpha }_1\xi +{\alpha }_2{\xi }^2\right){\mathrm{\Omega }}^2\right\}\eta =0.$

Given that the contraction angle of the conical pipe is very small, that is, ($\gamma \xi$) is sufficiently small so that we can use Tylor series here as:

23

$\frac{1}{1-\gamma \xi }=\sum _{n=0}^{\infty }{\left(\gamma \xi \right)}^n\approx 1+\gamma \xi +{\left(\gamma \xi \right)}^2+{\left(\gamma \xi \right)}^3+{\left(\gamma \xi \right)}^4+{\left(\gamma \xi \right)}^5+{\left(\gamma \xi \right)}^6,$
$\frac{1}{{(1-\gamma \xi )}^2}\approx 1+2\gamma \xi +3{\left(\gamma \xi \right)}^2+{4\left(\gamma \xi \right)}^3+{5\left(\gamma \xi \right)}^4+{6\left(\gamma \xi \right)}^5+7{\left(\gamma \xi \right)}^6,$
$\frac{1}{{(1-\gamma \xi )}^3}\approx 1+3\gamma \xi +6{\left(\gamma \xi \right)}^2+{10\left(\gamma \xi \right)}^3+{15\left(\gamma \xi \right)}^4+{21\left(\gamma \xi \right)}^5+{28\left(\gamma \xi \right)}^6,$
$\frac{1}{{(1-\gamma \xi )}^4}\approx 1+4\gamma \xi +10{\left(\gamma \xi \right)}^2+{20\left(\gamma \xi \right)}^3+35{\left(\gamma \xi \right)}^4+56{\left(\gamma \xi \right)}^5+{84\left(\gamma \xi \right)}^6,$
$\frac{1}{{(1-\gamma \xi )}^5}\approx 1+5\gamma \xi +15{\left(\gamma \xi \right)}^2+{35\left(\gamma \xi \right)}^3+{70\left(\gamma \xi \right)}^4+{126\left(\gamma \xi \right)}^5+{210\left(\gamma \xi \right)}^6.$

Here, we apply Frobenius method to analyze Eq. (22). Given that the solution of Eq. (22) takes the following form:

24

$\eta \left(\xi \right)=\sum _{n=0}^{\infty }{a}_n{\xi }^{n+k},a_0\ne 0.$

In which $a_n$ is the coefficient of the polynomial. The derivatives of $\eta \left(\xi \right)$ are:

25

$\frac{d\eta }{d\xi }=\sum _{n=0}^{\infty }{a}_n\left(n+k\right){\xi }^{n+k-1},\frac{d^2\eta }{d{\xi }^2}=\sum _{n=0}^{\infty }{a}_n\left(n+k\right)\left(n+k-1\right){\xi }^{n+k-2},$
$\frac{d^3\eta }{d{\xi }^3}=\sum _{n=0}^{\infty }{a}_n\left(n+k\right)\left(n+k-1\right)\left(n+k-2\right){\xi }^{n+k-3},$
$\frac{d^4\eta }{d{\xi }^4}=\sum _{n=0}^{\infty }{a}_n(n+k)(n+k-1)(n+k-2)(n+k-3){\xi }^{n+k-4}.$

When $n=0$, substituting Eq. (24) and Eq. (25) into Eq. (22), because that the lowest power of $\xi$ equals to zero, so the following equation can be obtained:

26

$a_{0}k(k-1)(k-2)(k-3)=0.$

Because that $a_0\ne 0$, so we obtain:

27

$k(k-1)(k-2)(k-3)=0.$

Eq. (27) is namely indicial equation and it plays an important role in the analysis [10]. Obviously, Eq. (27) have four roots, i.e. $k_j=$0, 1, 2, 3 ($j=$1, 2, 3, 4) respectively. To each $k$, only one coefficient $a_n\left(k\right)$ is related. Substituting Eq. (24) and Eq. (25) into Eq. (22), then all the coefficients can be obtained:

28a

$V_0=0,$

28b

$a_1\left(k\right)=-\frac{{\beta }_1\left(k-1\right)}{k+1}{a}_0\left(k\right),$

28c

$a_2\left(k\right)=-\frac{{\beta }_{1}k}{k+2}{a}_1\left(k\right)-\frac{{\beta }_{2}k\left(k-1\right)+{\zeta }_1}{\left(k+2\right)\left(k+1\right)}{a}_0\left(k\right),$

28d

$a_3\left(k\right)=-\frac{{\beta }_1\left(k+1\right)}{k+3}{a}_2\left(k\right)-\frac{{\beta }_{2}k\left(k+1\right)+{\zeta }_1}{\left(k+3\right)\left(k+2\right)}{a}_1\left(k\right)$
$-\frac{{\beta }_3\left(k+1\right)k\left(k-1\right)+{\alpha }_1{\zeta }_{1}k+4{\zeta }_1\gamma \left(k-1\right)+\left({\zeta }_1+{\zeta }_2\right)}{\left(k+3\right)\left(k+2\right)\left(k+1\right)}{a}_0\left(k\right),$

28e

$a_4\left(k\right)=-\frac{{\beta }_1\left(k+2\right)}{k+4}{a}_3\left(k\right)-\frac{{\beta }_2\left(k+2\right)\left(k+1\right)+{\zeta }_1}{\left(k+4\right)\left(k+3\right)}{a}_2\left(k\right)$
$-\frac{{\beta }_3\left(k+2\right)\left(k+1\right)k+{\alpha }_1{\zeta }_1\left(k+1\right)+4\gamma {\zeta }_{1}k+\left({\zeta }_1+{\zeta }_2\right)}{\left(k+4\right)\left(k+3\right)\left(k+2\right)}{a}_1\left(k\right)-\frac{\left[\begin{array}{c}{\beta }_4\left(k+2\right)\left(k+1\right)k\left(k-1\right)+\left(4{\alpha }_1\gamma +10{\gamma }^2+{\alpha }_2\right){\zeta }_{1}k\left(k-1\right)\\ +\left[\left(4{\alpha }_1\gamma +2{\alpha }_2+5\gamma +{\alpha }_1\right){\zeta }_1+2\gamma {\zeta }_2\right]k+\left({\alpha }_1{\zeta }_3+{\zeta }_4-{\zeta }_5\right)\end{array}\right]}{(k+4)(k+3)(k+2)(k+1)}{a}_0\left(k\right),$

28f

$a_{n+5}\left(k\right)=-\frac{{\beta }_1\left(k+n+3\right)}{k+n+5}{a}_{n+4}\left(k\right)-\frac{{\beta }_2\left(k+n+3\right)\left(k+n+2\right)+{\zeta }_1}{\left(k+n+5\right)\left(k+n+4\right)}{a}_{n+3}\left(k\right)$
$-\frac{\left[\begin{array}{c}{\beta }_3\left(k+n+3\right)\left(k+n+2\right)\left(k+n+1\right)\\ +{\alpha }_1{\zeta }_1\left(k+n+2\right)+4\gamma {\zeta }_1\left(k+n+1\right)+\left({\zeta }_1+{\zeta }_2\right)\end{array}\right]}{\left(k+n+5\right)\left(k+n+4\right)\left(k+n+3\right)}{a}_{n+2}\left(k\right)-\frac{\left\{\begin{array}{c}{\beta }_4\left(k+n+3\right)\left(k+n+2\right)\left(k+n+1\right)\left(k+n\right)\\ +\left(4{\alpha }_1\gamma {\zeta }_1+10{\gamma }^2{\zeta }_1+{\alpha }_2\right)\left(k+n+1\right)\left(k+n\right)\\ +\left[\left(4{\alpha }_1\gamma +2{\alpha }_2+5\gamma +{\alpha }_1\right){\zeta }_1+2\gamma {\zeta }_2\right]\left(k+n+1\right)+\left({\alpha }_1{\zeta }_3+{\zeta }_4-{\zeta }_5\right)\end{array}\right\}}{\left(k+n+5\right)\left(k+n+4\right)\left(k+n+3\right)\left(k+n+2\right)}{a}_{n+1}\left(k\right)-\frac{{\zeta }_3\left(2{\alpha }_2+2{\alpha }_1\gamma \right)+{\zeta }_4\left(3\gamma +{\alpha }_1\right)-{\zeta }_5{\alpha }_1}{\left(k+n+5\right)\left(k+n+4\right)\left(k+n+3\right)\left(k+n+2\right)}{a}_n\left(k\right),$

where:

29

${\zeta }_1=v_0^2,{\zeta }_2=i\mathrm{\Omega }{v}_0\sqrt{\beta },{\zeta }_3=i\mathrm{\Omega }\sqrt{\beta },{\zeta }_4=i\mathrm{\Omega }\sqrt{\beta }\gamma ,{\zeta }_5={\mathrm{\Omega }}^2.$

Here, the high orders of $\gamma$ have been neglected. After getting these coefficients, the solution of Eq. (22) can be written in the following form:

30

$\eta \left(\xi \right)=A_1{f}_1\left(\xi \right)+A_2{f}_2\left(\xi \right)+A_3{f}_3\left(\xi \right)+A_4{f}_4\left(\xi \right),$

where $A_j$ is constant. The expression of $f_j\left(\xi \right)$ is as follow:

31

$f_1{\left(\xi \right)=a}_0{(k}_1{)+a}_1{(k}_1{)\xi +a}_2{(k}_1{)\xi }^2{+a}_3{(k}_1{)\xi }^3+\dots ,$
$f_2{\left(\xi \right)=a}_0{(k}_2{)\xi +a}_1{(k}_2{)\xi }^2{+a}_2{(k}_2{)\xi }^3{+a}_3{(k}_2{)\xi }^4+\dots ,$
$f_3{\left(\xi \right)=a}_0{(k}_3{)\xi }^2{+a}_1{(k}_3{)\xi }^3{+a}_2{(k}_3{)\xi }^4{+a}_3{(k}_3{)\xi }^5+\dots ,$
$f_4{\left(\xi \right)=a}_0{(k}_4{)\xi }^3{+a}_1{(k}_4{)\xi }^4{+a}_2{(k}_4{)\xi }^5{+a}_3{(k}_4{)\xi }^6+\dots .$

Fig. 3The sketch of node force

4. Examples

4.1. Uniform pipe conveying fluid

In order to validate the proposed method, we calculated the natural frequencies of a uniform pipe conveying fluid. The uniform pipe is actually the case of variable diameter pipe with section variation coefficient $\gamma =\text{0}$. The pipe and fluid parameters are chosen the same as previous reference, that is, Yang’s modulus is $E=\text{210}\text{GPa}$, the outer and inner diameters at origin of the pipe are $D_i=\text{324}\text{mm}$ and $d_i=\text{292}\text{mm}$ respectively, and the pipe span is $L=\text{32}\text{m}$, the density of pipe material is ${\rho }_p=\text{8200}{\text{kg/m}}^{\text{3}}$, the density of fluid is ${\rho }_f=\text{908.2}{\text{kg/m}}^{\text{3}}$ [11]. The first five orders of natural frequencies are computed under three different velocities, i.e. $V=\text{0}$, $V=\text{15}\text{m/s}$, $V=\text{25}\text{m/s}$ respectively. And comparisons are presented between the results obtained by the proposed method and that published in the research of Housner [11]. The computation results of natural frequencies and error are listed in Table 1.

Table 1Natural frequencies of uniform simply supported pipe conveying fluid under different fluid velocities (rad/s)

Fluid velocity	Natural frequency
		${\omega }_1$	${\omega }_2$	${\omega }_3$	${\omega }_4$	${\omega }_5$
$V=\text{0}$	This paper	4.373	17.493	39.359	69.971	109.330
	Reference	4.3732	17.4928	39.3587	–	–
	error (%)	0.004	0.001	0.0007	–	–
$V=\text{15}\text{m/s}$	This paper	4.287	17.417	39.286	69.900	109.259
	Reference	4.2870	17.4171	39.2858	–	–
	error (%)	0	0.0006	0.0005	–	–
$V=\text{25}\text{m/s}$	This paper	4.131	17.282	39.156	69.772	109.133
	References	4.1293	17.2816	39.1559	–	–
	error (%)	0.04	0.002	0.0003	–	–

It can be found from Table 1 that a well agreement is obtained between the results in this paper and that in reference.

Fig. 4Curves of lg(h)-frequency under different fluid velocities: a) V=0, b) V=15 m/s, c) V=25 m/s

a)

b)

c)

Fig. 4 shows the curve of $\text{lg}h\left(\omega \right)$ versus frequency $\omega$ under different fluid velocities, and the natural frequencies have been marked in the figure. It should be noted here that when the slope of the curve is negative infinite, the real part and the imaginary part of $h\left(\omega \right)$ are both zeros. In another word, the frequencies related to the negative cuspidal points are the natural frequencies of the fluid-conveying pipe.