非线性输流曲管面内振动的样条小波有限元方法研究

曹建华， 黄颖青， 仝国军， 赵年顺

(1. 黄山学院 机电工程学院，安徽 黄山 245021; 2. 中国科学技术大学 工程科学学院 近代力学系，合肥 230027; 3. 河北水利电力学院 岩土工程安全与变形控制重点实验室，河北 沧州 061001)

由于小波具有多尺度和多分辨的优越特性，基于小波的数值方法，尤其是小波有限元方法，在工程和科学计算的应用研究吸引了众多国内外学者，并取得了巨大的进步。

在小波有限元研究方面，Williams等[1-2]做了很多基础性和开创性的研究工作，并研究第二代小波多分辨在有限元中的应用。国内众多学者在小波有限元方面做了很多工作。Chen等[3]应用二维多尺度区间样条小波单元于自适应有限元分析中，并用算例验证了其在奇异性问题上的有效性和可靠性。Wang等[4]又采用第二代小波构建小波有限元，并应用于偏微分的求解算例上。Xiang等[5]采用区间样条小波有限元，求解了转子动力学中转子轴承系统的频率问题。Xiang等[6]也用同样的方法，推导了圆锥厚壳的有限元矩阵，计算其模态数据，并用小波分解和支持向量机进行了探伤分析。

小波有限元与其他传统方法、新技术的结合，能够提高小波有限元方法的精度，缩短其计算时间。Shen等[7]结合小波有限单元和基于FFT(fast Fourier transform)的谱分析法，构建了高精度动刚度矩阵，研究了有裂纹或者脱落的杆和梁的波动特性。Zhang等[8]将小波有限单元和拉普拉斯变换相结合，提出了一种既能减少单元数又不用减小时间间隔的新方法，并将其应用于一维杆状结构中的超声波模拟。Hao等[9]采用区间样条小波构建复合板单元，基于图形处理单元(graphics processing unit,GPU)，实现一种小波板单元并行计算技术，应用于复合材料层合板的结构健康监测系统中，并指出基于GPU的并行计算比基于CPU的计算快了140多倍。

由以上可知，小波有限元计算可靠，应用广泛，但文献大多集中在线性计算方面，本文将以非线性输流曲管为研究对象，应用小波有限元求解其振动特性。输流管道在工业中应用广泛，受到学者们的关注。Misra等[10-11]采用传统有限元分别求解了基于轴线不可伸缩理论和轴线可伸缩理论的输流曲管振动问题。Jung等[12]采用非线性应变推导了轴线可伸长的输流曲管面内振动微分方程，并用传统有限元求解了系统频率。李宝辉等[13]采用波动法求解输流曲管的面内振动。考虑稳态组合力，Dai等[14]应用传递矩阵法分析了空间输流管道的流致振动。uczko等[15-16]采用伽辽金法分析了各类螺旋的三维输流曲管流致振动问题，也分析了流速和曲率对平面输流曲管的振动模态和频率的影响。Qu等[17]结合有限元和随机振动离散法，研究了受随机激励的液压曲管的振动特性。Zhou等[18]应用有限元研究了初始几何缺陷对微曲悬臂输流管道静平衡的影响，且提出其临界流速大小取决于静平衡构形。Ye等[19]采用伽辽金法、微分积分单元法和离散傅里叶法分析研究了微曲输流管道的超临界动力特性。基于尺度为4、阶数为6的区间样条小波，文献[20]构建了一维小波曲管单元，计算了轴线不可伸长的输流曲管频率。

本文详述基于曲管轴线可伸长的非线性输流曲管面内振动微分方程的推导过程，构建小波曲管单元，并应用小波有限元离散其微分方程，分析输流曲管的频率、模态和动力时间响应。

1 非线性输流曲管面内振动微分方程的推导

1.1 管道横截面上一点的应变公式推导

如图1所示，截取一段轴线弧长为ds、曲率半径为R的ds的曲管微元，在截面上放置一个局部坐标系，坐标轴为x和y轴，单位矢量分别为i和j，曲管微元离轴线为x处的弧长为L0。变形后，轴线弧长增大εEds，曲率半径为r，εE为轴线应变，离轴线为x处曲管微元的弧长为L1。则曲管微元离轴线为x处一点的应变εM[21]为

图1 曲管微元变形前和变形后示意图Fig.1 The diagram of the undeformed element and the deformed element of curved pipe

(1)

对于细长输流曲管管，x=R,x=r，因此，由式(1)可得管道横截面上一点的应变为

εM=εE+xΔκ

(2)

其中

(3)

图2 曲管轴线变形前和变形后示意图Fig.2 The diagram of the undeformed and deformed axial line of curved pipe

变形后的矢径r′表示为

r′=r+u(s)n+v(s)t

(4)

式中，u(s)，v(s)分别为轴线上一点沿切线单位矢量t和负法线单位矢量n方向的位移量，由图2可知，轴线的拉伸应变εE为

(5)

根据Hutchinson的课件(https://scholar.harv ard.edu/files/jiaweiyang/files/ curved_beams.pdf)，可以推出曲率的变化率Δκ为

(6)

根据拉格朗日应变公式和式(5)，拉伸应变ηE为

(7)

对于小应变、适度转动的曲管变形，εE≪1，可得

(8)

则轴线应变εE为

(9)

因此，由式(3)可知管道横截面上一点的应变公式εM为

即为[22]

(10)

式(10)与Jung等所引用的公式是一致的。

1.2 输流曲管面内流固耦合振动微分方程

本节基于式(10)，简述Jung等推导输流曲管面内流固耦合振动微分方程的过程。如图3所示两端固定的细长输流曲管，等截面且轴线为半圆，其轴线半径为R。假设流经管道的流体流速为U的柱塞流体，OXY坐标系为固定在机架上的惯性坐标系，θ为曲管截面与X轴正向的夹角。在截面上放置一个局部坐标系，坐标轴为x和y轴，单位矢量分别为i和j，则输流曲管面内径向位移ur(x,θ,t) 和周向位移uθ(x,θ,t)表达式分别为

图3 轴线为半圆的细长等截面输流曲管Fig.3 The geometry of semi-circular curved pipe conveying fluid

ur(x,θ,t)=u(θ,t)uθ(x,θ,t)=v(θ,t)+xφ

(11)

式中：t为时间；u与v分别为管道轴线上的点沿径向和周向位移；φ为管道振动时截面的面内转角，其表达式为式(8)。

管道中一点位移可以表示为

rp=(R+u+x)i+(v+xφ)j

(12)

管道一点的速度为

(13)

其管道中流体的速度为

(14)

其中t如图3所示，含义与图2一致。输流曲管的动能为

(15)

式中，mp和mf分别为单位长度管道和流体的质量。输流曲管其势能为

(16)

式中，εM采用非线性von Karman应变，其表达式如式(10)所示。采用线性应力表达式为

(17)

对于输流曲管，哈密顿原理表述为如下

(18)

其中，非保守力做功为

(19)

经过一系列变分，可得

(20)

(21)

式中，Q为轴力，其表达式为

(22)

两端固定的细长输流管道的边界条件为

θ=0和π,u=u′=v=0

(23)

2 样条小波有限元离散输流曲管面内振动微分方程

2.1 区间样条小波尺度函数

(24)

式中：j为尺度；m为0和1的多重节点数。基于上述序列点，定义在ti点处j尺度，m阶数的B-样条函数表达式为

(25)

(26)

因此, 在区间[0,1]上样条小波的尺度函数可写成向量形式

(27)

式中：ξ∈[0,1]； 且2j0≥2m-1，j≥j0。

在小波数值计算中，进行积分生成各类矩阵，主要有两类方法：一是直接用数学软件中函数生成样条小波尺度函数，应用积分函数生成单元矩阵，存储起来，然后在程序中直接调用即可； 二是本文采用的方法，编程生成样条函数，再构建样条小波尺度函数，利用高斯积分生成单元矩阵。在线性分析中矩阵不需要更新，两种方法均能胜任，但在非线性分析中，迭代过程中需要更新单元矩阵，第一种方法需要不停地调用数学软件的积分函数，花费太多时间，而第二种方法是采用高斯方法积分，故能够适用非线性计算。

2.2 样条小波有限元离散输流曲管面内振动微分方程

本文采用样条小波有限元离散输流曲管的非线性微分方程。在样条小波有限元中，本文采用尺度j1为3、阶数m1为4区间样条小条尺度函数(BSWI43)作为横向位移插值函数，采用尺度j2为2、阶数m2为3区间样条小条尺度函数(BSWI32)作为轴向位移的插值函数，两种区间样条小条尺度函数分别如图4所示。

图4 区间样条小波尺度函数Fig.4 The scaling function of BSWI43 and BSWI32 in the interval [0, 1]

输流曲管的轴向位移、横向位移表达式分别为

(28)

定义单元长度为le，单元径向位移和单元周向位移分别为

(29)

其中，ξi=(i-1)/2j,i=1,2，…,n+1=2j+1。

将式(28)代入式(29)，可以得到

(30)

其中

将式(30)代入式(28)

(31)

式中，Nu(ξ)=Φu(Re,u)-1和Nv(ξ)=Φv(Re,v)-1分别为径向位移和周向位移的样条小波有限元的形函数。令

(32)

式中， “′”和“″”分别为对ξ的一阶和二阶导数。根据传统有限元的离散过程，单元离散方程表示为如式(33)所示

(33)

其中

(34)

(35)

(36)

(37)

式中，“·”和“··”分别为对时间的一阶和二阶导数。将所有单元集合起来，可得全局离散常微分方程如下

(38)

式中，q为所有单元ue，ve位移的集合；Mg，Cg，Kg(q)分别为管道的全局质量矩阵、全局阻尼矩阵和全局刚度矩阵。由于Kg(q)为q的函数，因此是一个非线性问题。

3 非线性问题的求解方法

3.1 曲管的自然频率

为了求得静平衡位置处曲管的自然频率，需要将式(38)线性化，即为

(39)

式中，KTg为在静平衡位置处的切线刚度矩阵。

(40)

式中，q0为静平衡时的曲管变形，由式(41)求出

Kg(q0)q0=F

(41)

计算过程都是在施加边界条件之后进行的。求解式(41)非线性问题，主要采用直接迭代法和牛顿-拉斐森迭代法[24]，读者可参考相关文献，在此不再赘述。考虑q=QeΩt，式(39)改写为

(42)

计算式(42)的特征值，即可求出自然频率，定义无量纲频率和流速分别为

(43)

3.2 切线刚度矩阵的求解

计算切线刚度矩阵KTg的数值方法有前向差分法、中心差分法、复数步进法，差分法对步长要求有一定的要求，太小并不一定能提高精度，而复数步进法却没有这方面的限制[25-26]。复数步进法的原理是：将复变函数f(x+ih)在x处进行泰勒展开，f(x+ih)的表达式如下

(44)

式中：x+ih中的i为虚数单位；h为沿虚轴的小扰动步长。忽略高阶项，取其虚部系数，可得f(x)的导数为

(45)

根据上述原理，求解切线刚度矩阵的复数步长法表达式如下

(46)

式中：q+iheJ中的i为虚数单位； Im[·]为复数的虚部系数；Fint,I(U+iheJ)为以复数向量为变量的函数向量Fint的第I个分量。 由于没有差分法中的加减，精度有所提高。

4 数值算例

根据第3章的分析，作者编写了一套用于输流管道流固耦合振动分析的小波有限元程序，其中包含：生成区间样条小波尺度函数的基础程序、用于计算样条小波有限元单元矩阵的分段高斯积分程序和相关有限元的程序，应用于本文中输流曲管流固耦合面内振动分析。本章将小波有限元计算的曲管频率、模态和时间积分与文献、传统有限元进行对比，验证精度。

表1 当时输流曲管面内无量纲频率对比结果Tab.1 The comparison of the dimensionless natural frequencies of fluid-conveying curved pipe with

图5 输流曲管前3阶径向和周向模态Fig.5 The first three modes of curved pipe conveying fluid

图6 不同单元数计算下的输流曲管前3阶无量纲频率实部随流速变化曲线Fig.6 The first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity solved by different numbers of wavelet-based elements

图7 输流曲管前3阶无量纲频率实部随流速变化曲线比较Fig.7 The comparison of the first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity

以上数值结果是采用复数步进法生成切线刚度矩阵的，分别应用复数步进法、中心差分法、前向差分法和向后差分法计算输流曲管频率随流速变化曲线，如图8所示，曲线相差不大，中心差分法、前向差分法和向后差分法的计算结果保持一致。由局部放大图可知，复数步进法的计算结果稍微小一些。

图8 不同计算切线刚度矩阵方法的输流曲管前3阶频率 实部随流速变化曲线比较Fig.8 The comparison of the first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity solved by different method of computing the tangent stiffness

以上计算均是关于输流曲管的静平衡位置的计算结果，接着计算不同流速对于曲管静平衡位置的影响，如图9所示。由图9可知，随着流速增大，输流曲管静平衡时径向位移u和周向位移v的绝对值均随之增大，这与式(41)相符：流度越大，力F越大，位移也就越大。

图9 不同流速下输流曲管的静平衡位置Fig.9 The effects of different flow velocities on the static equilibrium of a fluid-conveying semi-circular pipe

图10 非线性输流曲管中点径向位移的时域动力 响应Fig.10 The time history of the midpoint’s radial displacement of curved pipe conveying fluid

图11 不同流速下非线性输流曲管中点径向位移的 时域动力响应Fig.11 The time history of the midpoint’s radial displacement of curved pipe conveying fluid with different flow velocities

图12 非线性输流曲管中点径向位移的时域 动力响应Fig.12 The time history of the midpoint’s radial displacement of curved pipe conveying fluid

表2 不同流速下时域动力响应计算时间比较Tab.2 The comparison of the computational time of the time response solved by wavelet-based finite element method and traditional finite element method

5 结 论

本文在前人文献的基础上，详细推导曲管截面上一点的应变公式，并简述了输流曲管面内振动非线性微分方程的推导过程，采用尺度为3、阶数为4的区间样条小条尺度函数作为径向位移的插值函数，尺度为2、阶数为3的区间样条小条尺度函数作为周向位移的插值函数，推导了区间样条小波曲管单元相关矩阵，并将其用于离散输流曲管面内振动微分方程，求解了输流曲管的频率、模态和时域响应，并与文献、传统有限元对比，进行验证。根据数值结果讨论，主要结论如下：

(1) 中心差分法、前向差分法和向后差分法计算切线刚度矩阵的频率结果保持一致，与复数步进法相比，数值相差不大，说明了复数步进法的精确性和可靠性。

(3) 与传统有限元对比，采用小波有限元计算的时域动力响应相差不大，且数值结果表明流速越大，位移时域响应就越大。

(4) 无论是频率计算还是时域动力响应计算，小波有限元计算时间都比传统有限元的少，收敛速度快。

综上所述，与传统有限元相比，小波有限元在非线性输流曲管面内振动计算结果准确，计算时间少，有着一定的优势。