任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究

振动工程学报Journal of Vibration Engineering

Vol. 33 No. 1

Feb. 2020

任意拉格朗日-欧拉描述的薄平板大幅扭转

振动气动特性研究

祝志文林君福唐意

2!王钦华1

(1.汕头大学土木与环境工程系，广东汕头515063# 2.中国建筑科学研究院，北京100013)

摘要：采用任意拉格朗日-欧拉（ArbitraryLagrangian-Eulerian, ALE)描述法，推导了二维流变区域不可压流体流 动的控制方程。基于刚性断面边界特点，简化得到了薄平板强迫振动绕流场控制方程，采用有限差分法开展了薄 平板大幅强迫扭转振动的CFD模拟。研究发现，薄平板小幅扭转振动气动力呈线性特性，但大幅扭转振动气动力 非线性显著，且气动力非线性随来流折算风速的增大而变得显著。从一个周期内气流做功来看，虽然大幅扭转振 动时薄平板转动轴到后缘是耗散气流能量的气动稳定区，但当折算风速较高后，上游侧的半个薄平板表面变为气 动不稳定区，此时气流对整个薄平板所作总功由负变为正，薄平板将从气流中吸收能量。研究认为，当折算风速较 高时，大幅扭转振动的薄平板将进人气动不稳定状态。

关键词：润激振动#薄平板#大幅振动#气动不稳定#任意Lagrange-Euleria)描述 中图分类号：TU311. 3

文献标志码：A

文章编号：1004-4523(2020)01-0012-12

DOI：10. 16385/j. cnki. issn. 1004-4523. 2020. 01. 002

不仅可跟踪流体质点的运动轨迹，而且能准确描述

引言

处于大气边界层中的大跨度桥梁受自然风的作 用。在某些条件下，桥梁或构件可能出现大幅度的 振动，如拉索、吊杆和主梁的大幅涡激振动[13]，风雨 振[45]和尾流驰振6、邻近或进人颤振状态的主梁大 幅颤振[7]等。这些大幅振动都是构件表面作为流体 边界的流动域发生显著变形的流体-固体耦合作用 问题）];对钝体外形的桥梁构件，或构件因大幅振动 产生较大的相对攻角效应，其气动力系统可能会因 这种大幅振动呈现显著的非线性）]。

在连续介质力学中，流体微元是微观上远大于 分子运动平均自由程而宏观上远小于所研究的流动 域的一团流体，并在宏观上将流体微元看作一个质 点。流体力学描述流体微元集合的运动状态，并通 过三种方法描述流体的运动，其中一种是拉格朗日

法（Lagrangian),这种方法跟踪流场每一质点的运 动，并获得流场流动参数随时间的变化，其特点是计 算网格与流体质点固定并同步运动，不存在质点与 网格间的相对运动。与欧拉描述法相比，这种方法

流体的移动界面。但当流动出现大变形时，将直接

导致计算网格的畸变，计算失败。欧拉法（Euleri- an)通过将感兴趣的流动区域离散为位置固定而数 量有限的网格点，进而描述流动参数的空间分布随 时间的变化。因其计算网格不随物体形位变化而改 变，所以能容易地处理流动区域的大变形，但描述流 变界面的运动需要引人非常复杂的映射关系，可能 导致较大的计算误差。任意拉格朗日-欧拉法就是 将拉格朗日法和欧拉法组合起来，取长补短。其特 点是描述流体运动的空间点位置既不固定，也不随 流体运动，而是以某种方式随时间变化，这种质点的 变化方式能反映流动域的实际变化，因而可准确地 描述流体的运动界面）0]，并维持单元的合理形 状）]，因而可解决只用欧拉法或只用拉格朗日法解 决不了的复杂运动边界的大变形流动问题）]。

桥梁构件大幅振动中，构件周围空气所占据的 实际空间一般是随时间变化的，并具有运动的界面。 可见，离散点空间分布固定的欧拉方法显然不适合 描述此类问题，因流动区域扭曲而引起计算网格的 畸形也会导致计算失败，因而合理描述应该采用任

收稿日期：2018-04-09 ;修订日期：2019-07-12

基金项目：国家重点基础研究发展计划（2015CB057701)国家自然科学基金资助项目（51878269,51778103)汕头大学科

研启动经费项目（NTF18014)

第\"期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究13

意拉格朗日-欧拉法。文献[11-12]通过CFD模拟 获得的振动刚性模型上的气动力和强迫振动位移时 程，建立气动力离散时间气动模型，开展了颤振导数 识别，但模型竖弯和扭转振动的幅值很小。文献

[13]开展了薄平板竖弯和扭转大幅振动CFD模拟 和气动力小波分析，发现在气动力时程中出现了比 薄平板振动频率高的二阶高频成分，也即大幅振动 的薄平板气动力为非线性，但因CFD基于Fluent 开展，而F+uent的动网格即使在均勻网格上也仅有 -------►X

0-------►x

一阶精度，因此研究结果的准确性需要验证。文 献[14]基于Volterra理论开展了薄平板非线性气 动力系统识别和CFD模拟，通过竖弯振动幅值为 1L和5%薄平板宽度的薄平板强迫振动CFD分 析，发现气动力的非线性效应并不显著，原因可能是 薄平板的振动幅度偏小。下面就桥梁气弹特性分析 常采用的薄平板为研究对象，针对二维不可压流动 特征，首先推导能适应任意流变区域的流体运动控 制方程，再根据刚性断面振动的特点，给出合适的空 间点运动方式和相应的流体运动控制方程。并以薄 平板为例研究其大幅振动的气动特征。

1

二维流体运动的任意拉格朗日-欧 拉描述方程

对二维流动问题，设Dphys(X，Y，）为在直角坐 标系OXY平面上随时间变化的真实物理域，同时 在另一个直角坐标系平面上设定一个虚拟的 计算域D_p(x，（）。这样，可将任意时刻真实物理 域的流动DPhys (X，Y，〇映射到一个不随时间变化的 虚拟计算区域D_p(x，（），如图1所示。由于真实 物理域Dpliys (X，Y，，）的非定常特性，其对计算域 D_p(x，（）的映射关系可能将随时间而变化。由 此，可建立任意时刻，，物理域内空间点（X，Y)和计 算域内点（工，（）的映射关系

(X = X(x，y，，）

(y = Y(x ,(,，）

通过对计算域D_p('，（）上感兴趣的物理量的描 述，达到描述时变物理域DPhys(X，Y，）上相应物理 量的目的。

设，时刻与计算区域点A 1('，（)对应的物理域 的点为私（'，（，〇,且物理域流体微元沿OX轴和 OY轴的速度分别为和V

，（，），如图2

所示。在下一时刻，+ !，与D_p ('，（）内A:('， ()点对应的物理域点变为('，（，+ /，），而DPhys (X，Y，)内找点处的流体微元此时运动到点.# ('+

什A埘刻物理域Dphys(U什

计

算

域

乃

_办，少）

图1

物理域和计算域的映射关系

Fig. 1

Mapping between physics domain and computa­tional domain

△'，（ P △(，，+ !〇,对应的速度为 U (' + △'，（ P △(，+/，），V (' + △'，（ P △(，，P △〇。而此时与 Dphys(X，Y，)内.2 点对应的 Aomp ('，（)点为.1('P

△x，（P △()。显然，从，到，P △，Dphys (X，Y，，）和

内有如下变化关系：

Acomp ('，（）内：Ai ('，（）---\"C1 (xP △'，y P△()

Aphys(X，Y，，）丄内

：

(x，（，，）--△ \"C2 (xP.

△'，

经过△〖时间后，APhys(X，Y，）上流体微元的速

度增量为

(AU = U(xP △' , (P △( ?，P △〇 一U(x，y，)

(△V = y (xP^x , (P^y ,，P △，~V(x，y，t)在(x，y，)处对式(2)等号右端第一项分别作泰 勒展开，并忽略二阶以上导数项，式(2)变为

△U (x，y，t) T $$x

—△xP $$— y

△( P

△t

3)

△V (x，y，，) T -$x —△xP $$—△ y

y P △t

这样，流体微元沿OX和OY轴的加速度分量分 别为：

/x=lim△-△-U t = —dt ---U lim$U$-------- △x x△t limUy$y，一〇△t

⑷

/y

tm△△t—$△--V = —$V ---$V △x t

$t

1 1 m$--------x△

t

1 1 $Vm —

一△y

$y，$ △

(5)

式中^x和△y是在△》时间内，流体微元从Apllys X，Y，t)上A运动到C，在A_p(x，y)上相对应点 的位移。

如将

(x，y)上点（x，y)所对应 Aphys(X，Y，

t)上的点称为参考点，因Apliys(X，Y，t)和A_p(x，

14振动工程学报第33卷

3〇的映射关系随时间发生变化，因此Dphys(X，Y“） 上的参考点位置也随时间发生改变。设时刻f流体 微元和参考点重合于A ( x，^，〇，！ Z时间后流体微 元运动到位置C(x + !x，y + ，+ !〖），参考点为A(x，3/，〇 (或 A(X，Y))位移到位置 B(x，3^+!〇(或

Y+AY))在 Dpliys(XY 上 AB 和 C

，

，

，

dU $XC?U $YdU

aX ) ^ 0 (U — ^) ^ 0 (V — ^) ‘

\\

,TT

、

| ,T/ 、

dt $t $X $t $Y

V \\ ( dX) dV I ( dY) dV

d t dt d X dtdY

(9)

对如图3所示的流体微元，设其密度为…如忽 略其彻体力（空气可以忽略），根据牛顿第二定律可

分别建立微元在OX和OY方向的平衡方程

三点位置如图2所示。

图2流体微元的相对运动

Fig. 2 Relative motion of fluid infinitesimal element

从图2可见，在以时间内，Dphys(X，Y，，）上流体微元相对于参考点的位移为，！T);流体微元在，时间内的位移为Ac，显然有sb=Ac—A\"B;设流体微元沿OX和OY轴的速度分量分别为卩和V，则Ac=(^!，V!〖）；Ab为以时间内的参考点位移，且A\"B=(AX，AY)。基于流体运动速度 和参考点速度，

可分别表示为

ArX = (U-$^)At)

ArY = (V-$Y )， (6dt

根据Dphys (X，Y，t )和LF_p (x，：y)的映射关系， W时间内D_p(x，^)上对应点的位移可用物理域 上流体微元的位移来描述，即

dx!x = $X!!rXX 00 $$YXArY

(7)

Ay = ^$X$X

ArX + (

AY

将式（6)代入式(7)可得

△x ) $Xx

(U —

0 石7(V — 77)/$Y

，

△

-X~0Y(v

$Y(8)y)Bu)!t再将上式代入式（3)和（4)并整理，有

/X )

_ d〇x I d#xY~dX dY(10)

f/Y

_ d#XY 0 d\"YdX 0 dY

对于不可压缩牛顿流体，流体内部的本构关系为

— 2dU ^°x=2$$x~p

\"Y = 2$ dV — 3

(11)

txy

：TYX

=$($U

式中

p

为流场内部压力;$为流体分子黏性系数。

如黏性系数为常数，将式（11)代入式（10)，有％

dU \\ ,T T 3 X、dU \\ ,T r dY

，十(U—眾)寇PV—沄)U

1 sp 丄,d2U丄d2U、_!dX+\"(dX2+dY2}

!12)

Vd

(U-1

,^VdXP(V—f

#

!$yPu(^^p$Y2 )

!13)

式中％=$为流体运动黏性系数。

图3

流体微元受到的应力

Fig. 3 Stress acting on fluid infinitesimal element

由于流场不可压，连续方程有如下形式

UdX

0dY

dV)〇

!14)

式（12)-(14)为采用任意拉格朗日描述法得到 的DPhys(X，Y，t)内流体运动的控制方程。与经典的 二维N-S方程相比，式（12)-(14)中的U(x，y，t)和

第1期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究15

V(:r，；y，f)不是物理域上不动点的流体流动速度，而 是与D_P (z，：y)内的点(z，：y)相对应的时变DPhys (X， y，d内的点（X，y)处流体微元的运动速度。需要指 出，式(12)-(14)方程中对流项系数不同于经典的二 维N-S方程，此为相对于运动参照系的流动速度。

市mT u'亦描并八r — $X$y 3X 3Y

未用Jacobl，又换，并令J=za^ — $7z，以及

dx1Y

其中

1 Y

Ur-( U-^-)

dtJ dy

(V-Ydt

)

1 X

J dy

)

(2 1)

1 YY

yr --(U-dX )+ ( y-dt$tJ dx

1 X

J dx

(2 2 )

基于式（18)-(22)可在二维计算域上实现对任 意流变区域内流体的运动描述，也即为描述任意流 变区域不可压流动的基本方程，其中映射关系式（1) dX =—了dydxdY =丄J X

dydy

丄Y

(15)

dXJ dx

dy1dX

dY

Jdx

可将式（12)-(14)变换到D_p(z，；y)，得到描述 二维任意流变区域的控制方程，以式（12 ) - (14 )为 例，如定义坐标变换量

一 (X)

-(()_

______XXY Y

$z dy 0 dx dy

(16)

7-

.(X)

$x

由链式法则推导一阶偏导数并定义，有

Z\\$2

x

X

2Xy+(y)+pX

p:

dY d 2 X 心 32X , d2X、

dy 3x2—2'xy+(y)—

3^ X $ 2 Y… d 2 Y)dy~32Y , dx 2\"—2'dxy+yy)

Q-

dY d X心 d2X ,

d2 X),Jhc('a dx—2'xy+(y)+3X( d2Y _ 3ZY ,d2Y、

x(ax—2'dxd二 ^

y +(d y 2(17)

将式（15)-(11)代人式（12)-(14)，经整理可分别得到动量方程和连续方程为：dU ,JJr 3U , yr dU ,1/Yd3 — dY3) —

dt

dx dy

5p dy dx

dx dy

uAxyU T 0

(18)

—PJ1 — pyr dV ,1 ( — dXdp, dXdp)_ dt dx dy Jp dy dx dx dy

uAxyV ：0

(19)dUdY _dUd¥—dVdXp dV_dX_dx dy dy dx dx dy dy dx

(0)

需根据所研究的流动区域的具体变形形式来确定。 当物理域变形特征简单时，对应的映射关系可简化， 此时，式（18)-(22)可能会得到明显简化，比如桥梁 断面的气动弹性问题，因桥梁断面的绝对刚性假设， 使得映射关系明显简化，如下文所述。

2

薄平板强迫振动绕流场控制方程

式(1)和(18)-(22)能描述任意变形物理域的流

体流动，比如带移动边界的流体动力学问题。移动边 界可为内流场或外流场流体外部边界的变形或移动，

也可为外流场绕流物体在流体力作用下其形位的改 变，这些均能导致与绕流物体接触的流体边界的移动 和变形。在进行桥梁断面一类绕流场模拟时，由于作 用在桥梁断面上的气动压力和摩阻力很小，因此桥梁 断面的变形很小，可近似将桥梁断面看作刚体。这 样，桥梁断面绕流场内边界的运动可用刚性桥梁断面 的整体运动来描述，也即采用刚体运动的描述方法。 在桥梁气动弹性分析中，通常用桥梁断面的竖弯和绕

某一参考轴的转动来描述桥梁断面上任意点的运动。

对振动的刚性断面，如果其计算域的远场边界 不变，将能恒定地给出远场边界条件，但由于计算域 内、外边界之间的物理区域在不断发生变化，如果进 行数值模拟，则需在每一时间步形成一次网格，如果 计算的时间步数量很大，由网格生成所带来的附加 计算量将很大，且难以保证每一步上重新生成的网

格的高质量。相比较，如果在离刚性断面足够远处 划定外边界，并对计算域划分高质量网格，数值模拟 时将计算域网格和刚性断面一起同步振动，就可只 在计算开始时形成一■次网格，而在每一■时间步上给 出相应的边界条件，如此处理，不仅使得数值模拟的 网格质量明显提高，也能省去每一时间步上网格的 重新生成。这样，需求解的物理区域只随时间发生 位置改变，形状并不改变。由于计算域事先给定，计 算域和物理域的对应关系也能明确。

取绝对坐标系OXY平面内任意点为物理域之

16振动工程学报第33卷

参考点，刚性断面的振动方式为相对于参考轴的平 动和绕参考点的转动。如取计算域的坐标原点在物 理域的映射点为参考点，并在此参考点上建立计算 域的相对坐标系奴^，当相对坐标系随刚性断面作 同步运动时，计算域D_p(x，^)和物理域Dphys(X，

存在下述映射关系

(X = X〇 0 'cos* + ysin*

Y _dY0

dt d，

~xcos*X *—ysin*X

*二

(28)

V0—+xsin*—+ycos*

将上述变换关系式代入式（1D)-(22)，有：

U + Ur d^PVr dy + — (d3cos*+|^sin*)'

$，$x$yP ^

d2U , d2U,

0(29)

dx2 dy2

u =Y$ —xsin* + ycos*

(23)

式中 *=*(，）为坐标系在'轴正向与坐标系在

dVPUr -PVr —十含(一xsin*pycos*)

，

xy

OX轴正向的夹角，如图4所示。

图4

物理域随时间的形位变化

Fig. 4

Deformation and displacement of physics domain with time

因计算域的形位并不发生变化，上述二个坐标 系存在如下对应关系％

:cos*

$Y =

—ysin*

$2X _

(24)

$2Y =

:0

d2Xdxdy$X_

$y

sin* Y

cos*

$2X_

$y20(25)

$2Y_

y

0$2X

dxdy

)0

此时，式（17)中坐标变换量简化为％ = 1，/?=0 y=1 J = 1，P = 0，Qt〇，且有

A _ d2 , d2

^xy_dX# 0dy

(26)

X，

X

d

d，

~xsin*X *+ycos*X * =

U0 —+xsin*P+ycos*

(27)

30)

dUt;-

x cos*+I ^;;dU— y sin. *—. | d V n /Q1)

x sm*+T；;—ycos*_ 0 (31)Ur _ (U—U0)c〇s* —(V—V0)sin*—+y (32)

Vr _ (U—U0)sin* —(V—V0)cos*++x (33)

式中U。

$X0

$，

，V0$Y，

0

+ _ *。

如将流体绝对速度之分量U和V分别投影到相对坐标系上，则得到的流体绝对速度沿x和y轴的分量〃和w分别为

9 _Ucos* — Vsin*

\\v _ Usin* 0 Vcos*

34)

也即物理域内流体微元的绝对速度在运动坐标系上的投影。

将式(34)作逆变换

(U = 9cos* 0 t;sin*[V = — usin* 0 vcos*

(35)

并分别计算U，V对相应变量的偏导数，代入式 (29)-(33)，并作如下运算：

式（29) X cos* 0 式（29) X (— sin*)

式（29) X sin*0 (29) Xcos*此时，将得到xy坐标系下流体运动的控制方程：

9 i du | du | 1 dp

，

+Ur dx+V^ dy+!$

x_

$2u $2u

dx2 dy2

36)

v i dv 丨 dv , 1 dp

%($d2xv2 +$d-Ty2

)++u37)$u $v

$x$(

38)

式中 ur _Ur

，vr _Vr

。此时式（32)，（33)变为：

ur _ u — U0 cos* 0 V0 sin* — +y (39)vr=v — U0 sin* — V0 cos* 0 cc

(40)

式中〜和Vr

分别是流体微元的绝对速度U和V 相对于运动坐标系的相对速度，在形成网格系统后，

第\"期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究17

即为流体微元相对于网格的运动速度。

对传统的桥梁气动弹性问题，刚性截面一般只 给定竖向平动和扭转二个自由度，如图5所示，式 (39)，（40)可简化为：

ur =u — (— — 0 sin* 0 +y ) (41)vr =v一(—0 cos*一cox )

(42)

通过选择合适的时间和空间步离散格式，对式

(36)-(38)和式（41)-(42)进行离散，在合理计算域 网格基础上，基于给定的初始和边界条件下，就可开 展刚性断面强迫振动的数值模拟。

图5

来流中竖弯和扭转振动的刚性断面

Fig. 5 Rigid section in incoming flow undergoing heave

and pitch vibration

3

薄平板扭转振动绕流场控制方程的 数值求解

3.1

CFD数值实现

对图6所示的长度为26的无厚度扭转振动刚

性薄平板（为便于查看，图示有厚度），其绕流控制方 程式本文采用非定常不可压流计算的二阶Projec­tion 算法 ）5* ， 即将上述控制方程为依次求解动 量方程得到中间速度场，由中间速度场求解压力场 的Poisson型方程，再由获得的压力场更新中间速 度场，得到修正的速度场方程，如此反复求解，反复

修正，并依据修正速度场在连续方程式（37)中产生 的残余值的大小来判断当前时间步流场变量的收敛 情况。当修正速度对连续方程的残余值满足收敛条 件时，求解进入下一时间步。

由二阶Projection算法得到的偏微分方程，空 间离散采用基于交错网格的有限差分法。在关于中 间速度场的动量方程求解上，左端动量项采用二阶 外插的Adams-Bashforth格式，右端黏性项采用二 阶中心差分格式;压力场方程左端项采用二阶中心 差分格式；时间导数采用二阶中心格式，得到的差分

图6扭转振动的薄平板

Fig. 6

Thin plate moving under pitch motion

方程均采用三对角矩阵算法求解，并通过Fortran 程序编程实现。

因刚性薄平板厚度为零，计算域网格采用直角 坐标系下的结构化网格，物面上下第一个网格高度 为0. 0016并沿垂直物面方向网格尺度逐渐增大， 但增长率不大于1. 05,以捕捉边界层的复杂流动，

并适度控制网格规模和非定常计算量。薄平板沿流

向均匀划分32个网格，但从薄平板前缘往计算域入 口，以及从薄平板后缘往计算域出口，采用不大于 145的网格增长率逐渐增大网格。计算域总网格 数为168X168,图7是计算域整体和薄平板周围的 网格划分，其中红色方框内有模拟的薄平板。

设图6所示的薄平板绕垂直OXY平面的薄平 板转动轴转动，振动方式为正弦振动，即

ait) =a0sin(2\"fat)

(43)

式中fa为薄平板扭转振动的频率;a(t)和a。分别 为扭转振动位移角和位移角幅值(°)，以薄平板前缘 抬头为正向，如图6所示。

从控制方程式（34)-(36)和式（39)，（40)的推导 可知，网格系统和薄平板刚性固定，同步振动。因此 在每一时间步，基于式（41)可得到网格和薄平板的 同步振动速度；薄平板表面采用无滑移条件，因此薄 平板物面流体的运动速度等于所在薄平板表面扭转 速度；对计算域的入口和上、下边界，因来流平行于

绝对参考系OXY的OX轴，需在每一时间步上分 别投影到运动参考系xy的两个坐标轴方向，成为 计算域对应边界的速度边界条件;计算域出口设置 为纽曼边界条件。对压力边界条件，在薄平板物面 和计算域入口、上下边界也采用纽曼边界条件。

定义薄平板流动的雷诺数为

Re =pU0B/ $

(44)

式中B = 2^。考虑作用在振动薄平板上的气动力

18振动工程学报第33卷

分别开展不同折算风速的CFD模拟，得到作用 在扭转振动薄平板上的气动力，利用式（46)和（47) 基于最小二乘法识别薄平板的4个颤振导数，并将 结果和Theodosen解析解作了比较，如图8所示。 可见二者具有一致的趋势线，高折算风速二者出现 较小偏差，可能是Theodosen解析解是基于无黏、 无流动分离的势流理论，但本文流动是真实的有黏 空气，且薄平板振动会导致流动分离，而黏性效应和

流动分离可能随折算风速的增大而增大。研究表 (a) Whole computational(a)计算域

domain

薄平板周围

(b) Around thin plate

图I

计算域网格布置

Fig. 7 Mesh arrangement

对雷诺数不敏感[16]，从降低不可压流非定常迭代求 解的计算量考虑，本文所有模拟的雷诺数均为 1200。

定义薄平板扭转振动的折算风速为

Vr=U〇/(fi)

(45)

对本文计算的全部折算风速、远方来流风速、薄平板 长度、网格划分和时间步长等均保持不变，仅通过改 变模型振动频率，通过式(45)获得不同的折算风速。

3.2

CFD方法验证

薄平板颤振导数有Theodosen解析解，通常用

于检验气弹分析CFD方法和程序的正确性。对单 自由度扭转强迫振动的薄平板，其非定常气动力与 相应颤振导数的表达式为[9]%

Fl= ^PU〇2b) (kH( <*〇 + k2H(*) (46)

M =1 pU〇 (2 b2 ) ( kA2 <8 0 k2 A3\" *)

(47)

式中 b = B/2;k = 2=b/U〇为折算频率；H(和 A(

2，3)为颤振导数。

明，本文CFD方法能较准确得到振动薄平板上的气 动力，数值方法是正确的。

图8识别的薄平板扭转颤振导数与Theodorsen解析 解的对比

Fig. 8 Identified flutter derivatives in comparison with

Theodorsen analytical solutions

4

薄平板大幅扭转强迫振动的能量 特征

4.1振动能量的定义

在来流中绕扭转轴作扭转强迫振动的薄平板，

第\"期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究19

定义薄平板表面的压力系数为

=pl(0• 7pU〇2)

振动时，其气动力系统具有良好的线性特性;另外，从

(48)

不同折算风速下的曲线来看，升力和扭矩间的相位差 随折算风速不同而有变化。

式中.3为压力系数，,0为来流风速。

定义在一个振动周期内，气流压力在薄平板上 微元dxz的做功为^3，,，对应的薄平板表面全部微 元一个周期内的总功为W3,也即

N N rTW3 = )i)W1 3，, = 2| 3l,Jdxlll^dt =

N0 L

i = 1 = 1

dXi3i,j^li

(49)

j式中K为作用在微元d'上的压力对薄平板扭转 轴的力臂&为薄平板扭转振动的角速度；N和L分 别为薄平板上下表面的单元总数和一个振动周期内 的计算时间步数;3j为第J个时间步作用在薄平板

微元d'上的压力。显然，气流对薄平板或表面微 元作正功表示薄平板从气流中吸收能量；如果气流 对薄平板作负功，表示薄平板消耗气流能量。

为方便比较，一个周期内扭转振动的总功采用 下述公式无量纲化

式中G为一个周期内竖弯振动的无量纲总功；T 为薄平板振动的周期。

如G为正，表示气流对薄平板作正功，薄平板 从气流中吸收能量，薄平板振动将被激发，也即薄平 板为气动不稳定状态;反之，如果G为负，气流对薄 平板作负功，薄平板消耗气流对其的做功，薄平板将 是气动稳定的。

4.2小幅扭转振动

定义作用在薄平板上的气动力系数为

(CL=FL/(0.5pU〇2B)

7

[Cm=M/(0. 5^,0+#)

式中Cl和CM分别是与升力Fl和扭矩M对应 的升力系数和扭矩系数。

定义无量纲时间为

$ =U

B

(52)

为对比薄平板大幅振动的气动力特征，下面给出 了薄平板小幅扭转振动（振幅0. 2[的气动力时程。

图9(a)和(b)分别是折算风速\\^ = 4和24时，一^个周 期内的气动力系数时程曲线。可见无论是升力还是 扭矩系数时程，其曲线具有明显的谐波特性，其波形 频率等于扭转强迫振动频率，可知当薄平板小幅扭转

(a) V=4

s

u-4->3py;ofuJo oyuIBUJ<3XPO

(b) V=24

图9

小幅扭转振动气动力系数时程

Fig. 9

Time histories of aerodynamic coefficients of thin plate undergoing small-amplitude pitch motion

图1〇为小幅扭转振动薄平板在一个完整周期

内4个关键时刻的压力等高线，对应从平衡位置出 发后达到最大正攻角位置（图10(a))，然后回到平 衡位置（图10(b))，再达到最大负攻角位置（图10 (c))，最后回到零攻角的平衡位置（图10(d))。可 见从一^个状态到另外一^个状态，压力等高线呈现有 规律地变化，薄平板前缘为高风压区，后缘为低风压 区。等高线分布规律的改变反映了薄平板的扭转振 动方向。

图11(a)和（b)分别为在一个扭转振动周期内， 气流在振动薄平板上表面各点、以及整个薄平板所 作的无量纲总功。可见薄平板前缘（L.E.)局部为 稳定区外，其后到转动轴（A. R.)为不稳定区，也即 气流对薄平板作正功，薄平板从气流中吸收能量，因 此为激励区;扭转轴到薄平板后缘（T. E.)气流做负 功，因而是耗散气流能量的稳定区，因此薄平板上激 励区与稳定区共存。图11 (b)给出了

二4-40，一-

个扭转周期内气流对薄平板所作的总功。可见在所

20振动工程学报第33卷

4.3大幅扭转振动

强迫薄平板发生大幅度的扭转振动，如扭转振 幅角为20[图12(a)-(d)分别为折算风速K = A，

Fig. 10 Pressure contour around thin plate in one cycle undergoing small-amplitude pitch motion at —r

图11小幅扭转振动薄平板一个周期内气动力总功随

折算风速的变化

Fig. 11

Aerodynamic work acting on thin plate in onecycle versus reduced wind speed when undergo­ing small-amplitude pitch motion

计算的折算风速范围内，一个周期内气流对扭转振 动的薄平板均做负功，也即振动薄平板消耗气流的 能量，因此小幅扭转振动的薄平板是稳定的。

s

yupIDy:j :30ou(?3oyP>>lOuJ3

<:240 280 320 360 400 440 480

(b) V=24

^U3p^j;3l*H1-0uouc3op>,oj3<

500 600 700 800

(c) t

V=44

syuspy;J 30ooyua3lup>>oJ3

< :700 800 900 1000 1100 1200 1300

t(d) V=6S

图12大幅扭转振动气动力系数时程

Fig. 12

Time histories of aerodynamic coefficients of thinplate undergoing large-amplitude pitch motion

第\"期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究21

24,44和68时，一个周期内的气动力系数时程曲 线。与图9相比，可见随着折算风速的提高，升力 和扭矩系数时程逐渐失去谐波特性，特别是高折 算风速。由于扭转振动为正弦强迫位移，可见当 薄平板大幅扭转振动时，其气动力系统不再为线 性，且随着折算风速的提高，其非线性特性越来越 显著。

图13为不同折算风速下大幅扭转振动的薄平 板，在一个振动周期内气流在振动薄平板上表面各 薄平板所作的总功，该总功的正负反映了扭转振动 薄平板的气动稳定性。

图14给出了一个扭转振动周期内气流对薄平 板所作总功随折算风速的变化，可见当折算风速 小于16时，一个周期内气流对薄平板所作总功为 负，也即薄平板消耗气流的能量#旦当折算风速大 于20后，气流对薄平板所作的总功由负变为正， 此时一个周期内扭转振动薄平板从气流中吸收能 量，这表明大幅度扭转振动的薄平板已进入气动 点所作的无量纲总功。与图11对比可见，在计算的 所有折算风速范围上，转动轴到薄平板后缘气流做 负功，因而是薄平板上耗散气流能量的稳定区，且在 低折算风速d = 4-16)下，薄平板表面的能量消耗 特征与小幅度扭转振动类似。但当折算风速高于 20后，迎风侧的半个薄平板上每个点在一个周期内 气流做的功均为正，即为不稳定区，气流对薄平板作 正功，薄平板从气流中吸收能量，因此薄平板上的激 励区与稳定区共存，而一个周期内气流对整个薄平 板表面所有点的总功之和决定了一个周期内气流对

Fig. 13 Aerodynamic work acting on thin plate in one

cycle when undergoing large-amplitude pitchmotion

不稳定状态。

Fig. 14 Aerodynamic work on thin plate versus reduced wind speed when undergoing large amplitudepitch motion

*结论

本文推导了二维流变区域的任意拉格朗日-欧

拉描述方程，基于刚性薄平板简化得到了薄平板强 迫振动绕流场控制方程，采用数值方法开展了薄平 板大幅强迫扭转振动的气动特征研究，可得到下述 结论：

1)

基于任意拉格朗日-欧拉描述的二维流体

动方程，通过建立计算域和物理域的随时间变化的 映射关系，可在不随时间变化的二维计算域上，采用 欧拉描述法实现对二维物理域上任意时变区域的流 体流动的描述。对薄平板扭转强迫振动，相应的映 射关系和控制方程明显简化。

2)

采用有限差分法开展了薄平板扭转振动

CFD模拟。从小幅扭转与大幅扭转振动气动力时 程的对比发现，大幅扭转振动气动力呈现非线性

特征，折算风速越高，气动力的非线性特征越 显著。

3) 大幅扭转振动时，薄平板转动轴到后缘是能

22振动工程学报

第33卷

量耗散区。但当折算风速高于20后，迎风侧的半个 薄平板从气流中吸收能量，为气动不稳定区#丑转振 动一个周期内气流对整个薄平板作的总功由负变为 正，此时薄平板从气流中吸收能量，表明大幅扭转强 迫振动的薄平板已进人气动不稳的状态。

参考文献：

[1]

Chen W L，Li H，Ou % P，et al. Numerical simulation 的气动力非线性特性研究[J].中国公路学报，2018,

31(1)：74-81.

ZhuZhiwen，ShiYaguang，Li Jianpeng. Research on nonlinear aerodynamics of flat steel box girder based on Volterra theory[J]. China Journal of Highway and Transport， 2018， 31(1)：74-81.

[10] Lohner R，Yang C. Improved ALE mesh velocities for

moving bodies [J]. Communications in Numerical 1996,12 (10) ：599-608.Methods inEngineering，

of vortex-induced vibrations of inclined cables under different wind profiles [%]. Journal of Bridge Engi­neering， 2013， 18(1)：42-53.

[2] Zuo D，Jones N P. Interpretation of field observations

of wind- and rain wind-induced stay cable vibrations [J]. Journal of Wind Engineering and Industrial Aero­dynamics， 2010， 98(2)：73-87.

[3] Zhu Q，Xu Y L，ZhuL D，et al. Vortex-induced vi­

bration analysis of long-span bridges with twin-box decks under non-uniformly distributed turbulent winds [J]. Journal of Wind Engineering and Industrial Aero­dynamics， 2018,172：31-41.

[4] Gu M， Du X Q. Experimental investigation of rain-

wind-induced vibration of cables in cable-stayed bridges and itsmitigation[J]. Journal of Wind Engineering and Industrial Aerodynamics，2005，93(1):79-95.

[5] Ge Yaojun，Chang Ying，Xu Linshan， et al. Experi­

mental investigation on spatial attitudes， dynamic characteristics and environmental conditions of rain- wind-induced vibration of stay cables with high-preti- sion raining simulator[J]. Journal of Fluids and Struc­tures ，2018,76:60-83.

[6] KitagawaT，OhtaH. Numerical investigation on flow

around circular cylinders in tandem arrangement at a subcritical Reynolds number[J]. Journal of Fluids and Structure，2008，24(5)：680-699.[7]

Josef Malik. Sudden lateral asymmetry and torsional oscillations in the riginal Tacoma suspension bridge [J]. Journal of Sound and Vibration , 2013,332(15)： 3772-37.

[8] Hirt C W，Amsden A A , Cook J L. An arbitrary La-

grangian-Eulerian computing method for all flow speeds[J]. Journal of Computational Physics , 1974， 14:227-253.[]

祝志文，石亚光，李健朋.扁平箱梁基于Voltera理论

[1] 祝志文，顾明，陈政清.基于3211多阶跃激励CFD

模型的颤振导数识别研究[].振动工程学报，2008,

21(1)18-23.

Zhu Zhiwen， Gu Ming， Chen Zhengqing. Identifica­tion of flutter derivatives based on 3211 multi-step ex­citation of CFD model[J]. Journal of Vibration Engi­neering 2008 21(1)：18-23.

[2] 祝志文，顾明，陈政清.指数脉冲强迫激励CFD模型

运动的气动参数识别法[].振动工程学报，2007, 20

(2) :133-139.

Zhu Zhiwen， Gu Ming， Chen Zhengqing. System i­dentification of aerodynamic parameter based on CFD modeling and exponential pulse excitation[J]. Journal ofVibrationEngineering，2007，20(2) ： 133-139.[13] HuangLin，XuYou-Lin，Liao Hai-li. Nonlinear aero­

dynamic forces on thin flat plate： Numerical study[J]. Journal of Fluids and Structures，2014,44:182-194. [4]祝志文，袁涛，陈政清，等.基于CFD仿真的平板非

线性气动力系统特征研究[].湖南大学学报（自然科 学版），2017, 44(1)32-38.

ZHU Zhi-wen，Yuan Tao，CHEN Zheng-qing， et al. Investigation on characteristics of nonlinear aerody­namic system of thin plate based on CFD simulations [J]. Journal of Hunan University (Natural Sciences)，2017,44(1)32-38.

[15] Nenad Filipovic，Srboljub Mijailovic，Akira Tsuda，et

al. An implicit algorithm within the arbitrary Lagrang- ian-Eulerian formulation for solving incompressible fluid flow with large boundary motions [J]. Computer Methods in Applied Mechanics and Engineering，2006， 195(44-47)：6347-6361.

[16] Zhu Zhiwen，Chen Zhengqing， Gu Ming. CFD based

simulations of flutter characteristics of ideal thin plates with and without central slot [J]. Wind and Struc­tures ，2009,12(1) ： 1-19.

第1期祝志文，等：任意拉格朗日-欧拉描述的薄平板大幅扭转振动气动特性研究23

Investigation on aerodynamics of thin plate undergoing large-amplitude

torsional oscillation with Arbitrary-Lagrangian-Eulerian descriptions

ZHUZht-zven1，LIN Jun-fu1，TANG Yt2，WANG Qing-hua1

(1. Department of Civil and Environmental Engineering，Shantou University，Shantou 515063，China；

2. China Academy of Building Research，Beijing 100013, China)

Abstract: Based on the Arbitrary Lagrangian-Eulerian descriptions，the governing equations of incompressibl

dimensional deformed region are derived. Due to the boundary features of a rigid section，the governing equations of forced vi­bration flow around thin plate are simplified. The finite difference method is used to carry out the CFD simulations of the large- amplitude forced torsional vibration of a

thin

plate. The research

shows that

the aerodynamics of

plitude torsional vibration are linear. However，when the thin plate is forced to oscillate with large amplitude，aerodynamically significant non-linearity will downstream half plate is defined tion occurs. It is noted

present. Considering as stabilizing

the

zone

the

total

work acting on the moving thin by

energy

dissipation

plate

plate bwhen

characterized

that when reduced wind speed is high，the upstream half become

zone，with the total work being positive， which means the thin plate absorbs the energy from t he flow. It is concluded thatwhen the reduced wind speed is high，the aerodynamic instability occurs if the thin plate is forced to perform large-amplitudeoscillations.

Keywords： vortex-induced vibration； thin plate； large-amplitude oscillation； aerodynamic instability； ALE descriptions

作者简介：祝志文（1968-)，男，博士，教授，博士生导师&电话：（0754)86502985; E-mail:zhuzw@stu. edu.cn