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

## Liu Yongshou1, Zhang Zijun2, Li Baohui3, Gao Hangshan4

1, 2, 3, 4Department of Engineering Mechanics, Northwestern Polytechnical University, Xi’an, P. R. China

2Corresponding author

Journal of Vibroengineering, Vol. 16, Issue 2, 2014, p. 832-845.
Received 7 December 2013; received in revised form 10 January 2014; accepted 24 January 2014; published 31 March 2014

Copyright © 2014 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Views 188
Abstract.

The governing equation of flow-induced vibration is deduced in terms of Hamilton’s principle for a variable diameter pipe conveying axial steady flow. Frobenius method is applied to analyze the governing equation. After the recursion formulas of coefficients are obtained, dynamic stiffness method is proposed for free vibration analysis of the variable diameter pipe conveying fluid. In the example, the natural frequencies of uniform pipes conveying fluid are computed and comparisons are made to validate the dynamic stiffness method. Then, the natural frequencies and modal shapes are obtained for the variable diameter pipe conveying fluid with different section variation coefficients and fluid velocities.

Keywords: variable diameter pipe conveying fluid, dynamic stiffness method, natural frequencies, modal shapes.

#### 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.:

(1)

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. is the work done by non-conservative force, ${m}_{f}$ is the fluid mass per unit length, ${\stackrel{\to }{r}}_{e}$ is the position vector of pipe exit. ${\stackrel{\to }{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:

(2)
$\stackrel{\to }{{V}_{a}}\left(x\right)=\left(\stackrel{˙}{u}+V\left(x\right){u}^{\text{'}}\right){\stackrel{\to }{e}}_{1}+\left(\stackrel{˙}{w}+V\left(x\right){w}^{\text{'}}\right){\stackrel{\to }{e}}_{2}.$

In Eq. (2), $V\left(x\right)$ stands for the velocity of the fluid relative to the pipe. ${\stackrel{\to }{e}}_{1}$ and ${\stackrel{\to }{e}}_{2}$ are unit vectors in $x$- and $y$-directions. and 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:

(3)
$V\left(x\right){A}_{f}\left(x\right)=C,$

where $C$ is a constant.

Fig. 1. The displacement sketch of simply supported pipe

Fig. 2. The 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:

(4)
$V\left(x\right)=\frac{C}{{\left(1-\frac{\gamma x}{L}\right)}^{2}{d}_{i}^{2}}.$

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

(5)
$V\left(x\right)=\frac{1}{{\left(1-\frac{\gamma x}{L}\right)}^{2}}{V}_{0}.$

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:

(6)
${T}_{kf}=\frac{1}{2}\underset{0}{\overset{L}{\int }}{\rho }_{f}{A}_{f}\left(x\right){V}_{a}^{2}\left(x\right)dx=\frac{1}{2}\underset{0}{\overset{L}{\int }}{\rho }_{f}{A}_{f}\left(x\right)\left({\stackrel{˙}{w}}^{2}+{V}^{2}\left(x\right)+2\stackrel{˙}{w}{w}^{\text{'}}V\left(x\right)+2V\left(x\right)\stackrel{˙}{u}\right)dx.$

Here, the axial inextensible assumption was used, i.e. . 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:

(7)
${T}_{kp}=\frac{1}{2}\underset{0}{\overset{L}{\int }}{\rho }_{p}{A}_{p}\left(x\right){\stackrel{˙}{w}}^{2}dx=\frac{1}{2}\underset{0}{\overset{L}{\int }}{m}_{p}\left(x\right){\stackrel{˙}{w}}^{2}dx,$

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:

(8)
${T}_{pp}=\frac{1}{2}E\underset{0}{\overset{L}{\int }}{I}_{x}{{w}^{″}}^{2}dx.$

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:

(9)
$\delta W=0.$

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

(10)
$\delta \underset{{t}_{1}}{\overset{{t}_{2}}{\int }}\frac{1}{2}\underset{0}{\overset{L}{\int }}\left[{m}_{f}\left({\stackrel{˙}{w}}^{2}+{V}^{2}+2{\stackrel{˙}{w}w}^{\text{'}}V+2V\stackrel{˙}{u}\right)+{m}_{p}{\stackrel{˙}{w}}^{2}-E{I}_{x}{{w}^{″}}^{2}\right]dxdt-\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}{m}_{f}V\left({\stackrel{˙}{\stackrel{\to }{r}}}_{e}+V{\stackrel{\to }{e}}_{t}\right)\cdot \delta {\stackrel{\to }{r}}_{e}dt=0.$

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:

(11)

At time ${t}_{1}$ and ${t}_{2}\text{:}$

(12)
$\left\{\begin{array}{l}\delta {u}_{{t}_{1}}=\delta {u}_{{t}_{2}}=0,\\ \delta {w}_{{t}_{1}}=\delta {w}_{{t}_{2}}=0.\end{array}\right\$

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

(13)
$\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}\underset{0}{\overset{L}{\int }}\left[E\left({I}_{x}{w}^{″}{\right)}^{″}+\left({m}_{p}+{m}_{f}\right)\stackrel{¨}{w}+\left({m}_{f}V\stackrel{˙}{w}\text{'}+{m}_{f}\mathrm{\text{'}}V\stackrel{˙}{w}+{m}_{f}{V}^{\text{'}}\stackrel{˙}{w}\right)+\left({{m}^{\text{'}}}_{f}{w}^{\text{'}}{V}^{2}+{m}_{f}{w}^{″}{V}^{2}+2{m}_{f}V{V}^{\text{'}}{w}^{\text{'}}\right)\right]\delta wdxdt=0.$

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

(14)

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:

(15)
$\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}{m}_{f}V\left({\stackrel{˙}{\stackrel{\to }{r}}}_{e}+V{\stackrel{\to }{e}}_{t}\right)\cdot \delta {\stackrel{\to }{r}}_{e}dt=\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}{m}_{f}{V}^{2}\delta {u}_{e}dt=-\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}\underset{0}{\overset{L}{\int }}{m}_{f}{V}^{2}{w}^{\text{'}}\delta {w}^{\text{'}}dxdt=\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}\underset{0}{\overset{L}{\int }}\left({m}_{f}{w}^{\text{'}}{V}^{2}\right)\text{'}\delta wdxdt=\underset{{t}_{1}}{\overset{{t}_{2}}{\int }}\underset{0}{\overset{L}{\int }}\left({{m}^{\text{'}}}_{f}{w}^{\text{'}}{V}^{2}+{m}_{f}{w}^{″}{V}^{2}+2{m}_{f}V{V}^{\text{'}}{w}^{\text{'}}\right)\delta wdxdt.$

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

(16)
$E{\left({I}_{x}{w}^{″}\right)}^{\text{'}\text{'}}+\left({m}_{p}+{m}_{f}\right)\stackrel{¨}{w}+\left({m}_{f}V{\stackrel{˙}{w}}^{\text{'}}+{m}_{f}V\stackrel{˙}{w}+m{\text{'}}_{f}{V}^{\text{'}}\stackrel{˙}{w}\right)+\left({{m}^{\text{'}}}_{f}{w}^{\text{'}}{V}^{2}+{m}_{f}{w}^{″}{V}^{2}+2{m}_{f}V{V}^{\text{'}}{w}^{\text{'}}\right)=0.$

#### 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}}+2{m}_{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)

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}{{\left(1-\gamma \xi \right)}^{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}{{\left(1-\gamma \xi \right)}^{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}{{\left(1-\gamma \xi \right)}^{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}{{\left(1-\gamma \xi \right)}^{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)

In which ${a}_{n}$ is the coefficient of the polynomial. The derivatives of $\eta \left(\xi \right)$ are:

(25)

$\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}\left(n+k\right)\left(n+k-1\right)\left(n+k-2\right)\left(n+k-3\right){\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\left(k-1\right)\left(k-2\right)\left(k-3\right)=0.$

Because that ${a}_{0}\ne 0$, so we obtain:

(27)
$k\left(k-1\right)\left(k-2\right)\left(k-3\right)=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. 0, 1, 2, 3 (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)$
(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)$
(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)$

where:

(29)

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}{\left(k}_{1}{\right)+a}_{1}{\left(k}_{1}{\right)\xi +a}_{2}{\left(k}_{1}{\right)\xi }^{2}{+a}_{3}{\left(k}_{1}{\right)\xi }^{3}+\dots ,$
${f}_{2}{\left(\xi \right)=a}_{0}{\left(k}_{2}{\right)\xi +a}_{1}{\left(k}_{2}{\right)\xi }^{2}{+a}_{2}{\left(k}_{2}{\right)\xi }^{3}{+a}_{3}{\left(k}_{2}{\right)\xi }^{4}+\dots ,$
${f}_{3}{\left(\xi \right)=a}_{0}{\left(k}_{3}{\right)\xi }^{2}{+a}_{1}{\left(k}_{3}{\right)\xi }^{3}{+a}_{2}{\left(k}_{3}{\right)\xi }^{4}{+a}_{3}{\left(k}_{3}{\right)\xi }^{5}+\dots ,$
${f}_{4}{\left(\xi \right)=a}_{0}{\left(k}_{4}{\right)\xi }^{3}{+a}_{1}{\left(k}_{4}{\right)\xi }^{4}{+a}_{2}{\left(k}_{4}{\right)\xi }^{5}{+a}_{3}{\left(k}_{4}{\right)\xi }^{6}+\dots .$

Fig. 3. The sketch of node force

#### 3.2. Dynamic stiffness method

Here, we use dynamic stiffness method to analyze vibration of variable diameter pipe conveying fluid. As shown in Fig. 3, the node displacement of the pipe is:

(32)
${\delta }^{e}=PA\text{,}$

where:

(33)

$P=\left[\begin{array}{cccc}{a}_{0}\left({k}_{1}\right)& 0& 0& 0\\ {a}_{1}\left({k}_{1}\right)& {a}_{0}\left({k}_{2}\right)& 0& 0\\ {f}_{1}\left(1\right)& {f}_{2}\left(1\right)& {f}_{3}\left(1\right)& {f}_{4}\left(1\right)\\ {f}_{1}^{\mathrm{\text{'}}}\left(1\right)& {f}_{2}^{\mathrm{\text{'}}}\left(1\right)& {f}_{3}^{\mathrm{\text{'}}}\left(1\right)& {f}_{4}^{\text{'}}\left(1\right)\end{array}\right].$

As shown in Fig. 3, the shear force and bending moment can be expressed as:

(34)
$M=E{I}_{x}\frac{{\partial }^{2}w}{\partial {x}^{2}},$
$Q=-\frac{\partial M}{\partial x}=-E{I}_{x}\frac{{\partial }^{3}w}{\partial {x}^{3}}-E\frac{\partial {I}_{x}}{\partial x}\frac{{\partial }^{2}w}{\partial {x}^{2}}.$

Introducing the dimensionless form of shear force and bending moment as follows:

(35)
$\stackrel{-}{M}=\frac{ML}{E{I}_{i}}=\left(1+{\beta }_{1}\xi +{\beta }_{2}{\xi }^{2}+{\beta }_{3}{\xi }^{3}+{\beta }_{4}{\xi }^{4}\right)\frac{{d}^{2}\eta }{d{\xi }^{2}},$
$\stackrel{-}{Q}=\frac{Q{L}^{2}}{E{I}_{i}}=-\left(1+{\beta }_{1}\xi +{\beta }_{2}{\xi }^{2}+{\beta }_{3}{\xi }^{3}+{\beta }_{4}{\xi }^{4}\right)\frac{{d}^{3}\eta }{d{\xi }^{3}}$

The directions of $M$, $Q$ are shown in Fig. 3, then it is easy to get the vector of node force:

(36)
$F={\left[\begin{array}{cccc}{\stackrel{-}{Q}}_{1}& {\stackrel{-}{M}}_{1}& {\stackrel{-}{Q}}_{2}& {\stackrel{-}{M}}_{2}\end{array}\right]}^{T},$

where ${\stackrel{-}{Q}}_{1}=-\stackrel{-}{Q}\left(0\right)$, ${\stackrel{-}{M}}_{1}=-\stackrel{-}{M}\left(0\right)$, ${\stackrel{-}{Q}}_{2}=\stackrel{-}{Q}\left(1\right)$, ${\stackrel{-}{M}}_{2}=\stackrel{-}{M}\left(1\right)$.

The relation between vector of node force and vector $A$ can now be expressed as follow:

(37)
$F=BA,$

where:

(38)

In which ${\sigma }_{1}=1+{\beta }_{1}+{\beta }_{2}+{\beta }_{3}+{\beta }_{4}$, ${\sigma }_{2}={\beta }_{1}+2{\beta }_{2}+3{\beta }_{3}+4{\beta }_{4}$, and what should be noted here is that, stands for the derivative with respect to $\xi$, which is different from that in Section 2.

According to Eq. (32), we can get the coefficient vector as follow:

(39)
$A={P}^{-1}{\delta }^{e}.$

Now, substituting Eq. (39) into Eq. (37), the following result will be gained:

(40)
$F=B{P}^{-1}{\delta }^{e}=S{\delta }^{e}.$

Here, the matrix $S$ is named as the dynamic stiffness matrix of the system. Once the dynamic stiffness matrix is obtained, the natural frequencies and modal shapes can be easily obtained.

Here, we take the simply supported pipe for an example to illuminate how to compute natural frequencies of the pipe through dynamic stiffness method. The boundary conditions of simply supported pipe can be written as:

(41)

Rewriting Eq. (41) in the matrix form, we have:

(42)
$\left[\begin{array}{c}{\stackrel{-}{Q}}_{1}\\ 0\\ {\stackrel{-}{Q}}_{2}\\ 0\end{array}\right]=\left[\begin{array}{cccc}{S}_{11}& {S}_{12}& {S}_{13}& {S}_{14}\\ {S}_{21}& {S}_{22}& {S}_{23}& {S}_{24}\\ {S}_{31}& {S}_{32}& {S}_{33}& {S}_{34}\\ {S}_{41}& {S}_{42}& {S}_{43}& {S}_{44}\end{array}\right]\left[\begin{array}{c}0\\ {\eta }_{1}^{\text{'}}\\ 0\\ {\eta }_{2}^{\text{'}}\end{array}\right].$

The following characteristic equation can be obtained from Eq. (42):

(43)
$\left[\begin{array}{cc}{S}_{22}& {S}_{24}\\ {S}_{42}& {S}_{44}\end{array}\right]\left[\begin{array}{c}{\eta }_{1}^{\mathrm{\text{'}}}\\ {\eta }_{2}^{\mathrm{\text{'}}}\end{array}\right]=0.$

It is easy to find from Eq. (43) that the following equation must be satisfied:

(44)
$h\left(\omega \right)=\mathrm{d}\mathrm{e}\mathrm{t}\left|\begin{array}{cc}{S}_{22}& {S}_{24}\\ {S}_{42}& {S}_{44}\end{array}\right|=0.$

The natural frequencies can be obtained from Eq. (44). In the same way, we can get the characteristic equation for cantilevered pipe and clamped-pinned pipe conveying fluid as follows:

Cantilevered pipe (left clamped and right free):

(45)
$h\left(\omega \right)=\mathrm{d}\mathrm{e}\mathrm{t}\left|\begin{array}{cc}{S}_{33}& {S}_{34}\\ {S}_{43}& {S}_{44}\end{array}\right|=0.$

Clamped-pinned pipe (left clamped and right pinned):

(46)
$h\left(\omega \right)=\mathrm{d}\mathrm{e}\mathrm{t}\left|\begin{array}{ccc}{S}_{11}& {S}_{12}& {S}_{14}\\ {S}_{21}& {S}_{22}& {S}_{24}\\ {S}_{41}& {S}_{42}& {S}_{44}\end{array}\right|=0.$

#### 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 , the outer and inner diameters at origin of the pipe are and respectively, and the pipe span is , the density of pipe material is , the density of fluid is [11]. The first five orders of natural frequencies are computed under three different velocities, i.e. $V=\text{0}$, , 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 1. Natural 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 – – 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 – – 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. 4. Curves of $\text{lg}\left(h\right)$-frequency under different fluid velocities: a) $V=\text{0}$, b) , c)

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.

#### 4.2. Variable diameter pipe conveying fluid

To simply supported variable diameter pipes conveying fluid, we choose three different diameter variation coefficients, that is 0.1, 0.2, 0.3, respectively. Given that the Yang’s modulus of the pipe material is , the outer and inner diameters of pipe at the origin of the pipe are and respectively. The pipe length is , the density of pipe material is , and the density of fluid is . The natural frequencies are computed under the following three different fluid velocities: $V=\text{0,}$. The results are listed in Table 2.

Table 2. Natural frequencies of variable diameter pipe conveying fluid with different diameter variation coefficients under different fluid velocities

 $\gamma$ and fluid velocities Natural frequencies ${\omega }_{1}$ ${\omega }_{2}$ ${\omega }_{3}$ ${\omega }_{4}$ ${\omega }_{5}$ $\gamma =\text{0.1}$ ${V}_{0}=\text{0}$ 3.54 14.14 31.82 56.56 88.36 3.10 13.80 31.52 56.26 88.08 – 11.92 29.84 54.68 86.54 $\gamma =\text{0.}2$ ${V}_{0}=\text{0}$ 3.34 13.48 30.32 53.88 84.18 2.88 13.12 29.98 53.56 83.88 – 11.10 28.20 51.88 82.84 $\gamma =\text{0.3}$ ${V}_{0}=\text{0}$ 3.16 12.92 29.04 51.62 80.58 2.66 12.52 28.68 51.28 80.42 – 10.30 26.76 49.44 77.98

From Table 2, we can find that to the same diameter variation coefficient $\gamma$, the natural frequencies decrease as the fluid velocity increases, and that the larger $\gamma$ we choose, the smaller the natural frequencies we get. The results obtained above agree with that in the researches of Hannoyer and Paidoussis [8].

Additionally, we notice that when the fluid velocity reaches 25 m/s, the first order natural frequency has already vanished. So we can conclude that such velocity has already exceeded the critical velocity of the pipe. Actually, the critical velocities of variable diameter pipe conveying fluid with different diameter variation coefficients are 0.1, 0.2, ; 0.3, , respectively, as shown in Fig. 5.

Fig. 5. Curves of first order natural frequency versus fluid velocity

Fig. 6. Buckling of the variable diameter pipe conveying fluid

When the fluid velocity reaches the critical value, the simply supported variable diameter pipe would yield buckling, as a form of static instability, as shown in Fig. 6. Further, the instability would result in the failure of the pipe.

It is obvious that the larger the diameter variation coefficient is, the smaller the critical velocity we get, which agree with the result in the researches of Hannoyer and Paidoussis [8].

Fig. 7 shows the curves of $\text{lg}h\left(\omega \right)$ versus frequency $\omega$ for variable diameter pipe conveying fluid with different diameter variation coefficients under velocity ${V}_{0}=\text{0}$ and .

Fig. 7. Curves of $\text{log}\text{(}h\right)$-frequency with different diameter variation ratio under different fluid velocities: a) $\gamma =\text{0.1}$, ${V}_{0}=\text{0}$; b) $\gamma =\text{0.1}$, ; c) $\gamma =\text{0.2}$, ${V}_{0}=\text{0}$; d) $\gamma =\text{0.2}$, ; e) $\gamma =\text{0.3}$, ${V}_{0}=\text{0}$; f) $\gamma =\text{0.3}$,

a)

b)

c)

d)

e)

f)

After we obtain the natural frequencies of the variable diameter pipe, the modal shapes can be obtained correspondingly. The modal shapes of variable diameter pipe conveying fluid are shown in Fig. 8 and Table 3.

Table 3. Dimensionless amplitudes of the modal shape under different flow velocities

 Fluid velocity Amplitudes 1st order 2nd order 3rd order 4th order ${V}_{0}=\text{0}$ 0.22 0.29 0.21 0.12 0.25 0.31 0.23 0.15 0.30 0.35 0.26 0.18

Table 3 shows that the amplitude of the modal shape increases with increasing flow velocity. This illuminate that the increasing fluid velocity can weaken the pipe stiffness.

Fig. 8. The first four orders modal shapes of variable diameter pipe conveying fluid with different diameter variation coefficients under different fluid velocities

#### 4.3. Variable diameter pipe with damping

Damping is ignored in above calculations. But in practice, damping cannot be neglected. Here we just consider the viscous damping, then the motion equation governing the vibration of the pipe is:

(47)
$E\left({I}_{x}{w}^{″}\right)\text{'}\text{'}+\left({m}_{p}+{m}_{f}\right)\stackrel{¨}{w}+\left({m}_{f}V{\stackrel{˙}{w}}^{\text{'}}+{m}_{f}V\stackrel{˙}{w}+m{\text{'}}_{f}{V}^{\text{'}}\stackrel{˙}{w}\right)+c\stackrel{˙}{w}+\left({{m}^{\text{'}}}_{f}{w}^{\text{'}}{V}^{2}+{m}_{f}{w}^{″}{V}^{2}+2{m}_{f}V{V}^{\text{'}}{w}^{\text{'}}\right)=0,$

where the term ($c\stackrel{˙}{w}$) is the viscous damping force, which may be introduced by immersing the pipe in liquid. Suppose that 0.05, then we can calculate the natural frequencies of the variable diameter pipe. A contrast is illustrated in Fig. 9 between damping pipe and non-damping pipe (0).

Fig. 9. First order natural frequencies of damping pipe and non-damping pipe with variable fluid velocity

From Fig. 9, we conclude that damping would decrease the natural frequencies and critical flow velocity of the system. This has been illustrated by Paidoussis [2] in the research for a uniform fluid-conveying pipe. Moreover, damping delays the decreasing speed of natural frequency versus flow velocity. That is, as seen in Fig. 9, there is just a small difference between the two critical flow velocities, but a relatively larger difference between the natural frequencies when flow velocity is zero.

#### 5. Conclusion

In this paper, the dynamic model was built by open system’s Hamilton’s principle on the basis of Euler-Bernoulli beam model and axial inextensible assumption. The motion equation of variable diameter pipe conveying fluid is more complex than that of uniform pipe. We employed Frobenius method to analyze the motion equation and proposed the dynamic stiffness method for free vibration analysis of variable diameter pipe conveying fluid. By using the presented method, the natural frequencies of a uniform pipe conveying fluid under different fluid velocities were obtained and the comparisons between our results and that in previous researches were made. A well agreement validates our method. The natural frequencies and modal shapes were obtained for variable diameter pipe conveying fluid with different diameter variation coefficients under different fluid velocities. We find that the natural frequencies and critical velocities of variable diameter pipe both decreased with increasing diameter variation coefficient.

#### Acknowledgements

This work was supported by Northwestern Polytechnical University Basic Research Fund (Grant No. JC201114) and Aeronautical Science Foundation of China (2011ZA53014).

#### References

1. Ashley H., HavilandG. Bending vibrations of a pipe line containing flowing fluid. Journal of Applied Mechanics-Transactions of the ASME, Vol. 17, Issue 3, 1950, p. 229-232. [Search CrossRef]
2. Paidoussis M. P. Fluid structure interactions: slender structures and axial flow. Volume 1. London, Academic Press, 1998. [Search CrossRef]
3. Doare O. Dissipation effect on local and global stability of fluid-conveying pipes. Journal of Sound and Vibration, Vol. 329, Issue 1, 2010, p. 72-83. [Search CrossRef]
4. Jung D., Chung J., Mazzoleni A. Dynamic stability of a semi-circular pipe conveying harmonically oscillating fluid. Journal of Sound and Vibration, Vol. 315, Issue 1-2, 2008, p. 100-117. [Search CrossRef]
5. Rinaldi S., et al. Dynamics of microscale pipes containing internal fluid flow: damping, frequency shift, and stability. Journal of Sound and Vibration, Vol. 329, Issue 8, 2010, p. 1081-1088. [Search CrossRef]
6. Huang Y. M., et al. Natural frequency analysis of ﬂuid conveying pipeline with different boundary conditions. Nuclear Engineering and Design, Vol. 240, Issue 3, 2010, p. 461-467. [Search CrossRef]
7. Ragulskis M., Fedaravicius A., Ragulskis K. Harmonic balance method for FEM analysis of fluid flow in a vibrating pipe. Communications in Numerical Methods in Engineering, Vol. 22, Issue 4, 2006, p. 347-356. [Search CrossRef]
8. Hannoyer M. J., Paidoussis M. P. Dynamics of slender tapered beams with internal or external axial flow 1 theory. Journal of Applied Mechanics-Transactions of the ASME, Vol. 46, Issue 1, 1979, p. 45-51. [Search CrossRef]
9. Hannoyer M. J., Paidoussis M. P. Dynamics of slender tapered beams with internal or external axial flow 2 experiments. Journal of Applied Mechanics-Transactions of the ASME, Vol. 46, Issue 1, 1979, p. 52-57. [Search CrossRef]
10. Arfken G. B., Weber H. J. Mathematical methods for physicists. San Diego, CA, Academic Press, 1995. [Search CrossRef]
11. Housner G. W. Bending vibrations of a pipe line containing flowing fluid. Journal of Applied Mechanics-Transactions of the ASME, Vol. 19, Issue 2, 1952, p. 205-208. [Search CrossRef]