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

祝志文 林君福 唐意 王钦华

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

关键词：涡激振动;薄平板;大幅振动;气动不稳定;任意Lagrange-Eulerian描述

中图分类号：TU311.3 文献标志码：A 文章编号：1004-4523（2020）01-0012-12

DOI：10.16385/j.cnki.issn.1004-4523.2020.01.002

引言

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

在连续介质力学中，流体微元是微观上远大于分子运动平均自由程而宏观上远小于所研究的流动域的一团流体，并在宏观上将流体微元看作一个质点。流体力学描述流体微元集合的运动状态，并通过三种方法描述流体的运动，其中一种是拉格朗日法（Lagrangian），这种方法跟踪流场每一质点的运动，并获得流场流动参数随时问的变化，其特点是计算网格与流体质点固定并同步运动，不存在质点与网格问的相对运动。与欧拉描述法相比，这种方法不仅可跟踪流体质点的运动轨迹，而且能准确描述流体的移动界面。但当流动出现大变形时，将直接导致计算网格的畸变，计算失败。欧拉法（Euleri-an）通过将感兴趣的流动区域离散为位置固定而数量有限的网格点，进而描述流动参数的空间分布随时问的变化。因其计算网格不随物体形位变化而改变，所以能容易地处理流动区域的大变形，但描述流变界面的运动需要引入非常复杂的映射关系，可能导致较大的计算误差。任意拉格朗日一欧拉法就是将拉格朗日法和欧拉法组合起来，取长补短。其特点是描述流体运动的空问点位置既不固定，也不随流体运动，而是以某种方式随时问变化，这种质点的变化方式能反映流动域的实际变化，因而可准确地描述流体的运动界面，并维持单元的合理形状，因而可解决只用欧拉法或只用拉格朗日法解决不了的复杂运动边界的大变形流动问题。

桥梁构件大幅振动中，构件周围空气所占据的实际空问一般是随时问变化的，并具有运动的界面。可见，离散点空问分布固定的欧拉方法显然不适合描述此类问题，因流动区域扭曲而引起计算网格的畸形也会导致计算失败，因而合理描述应该采用任意拉格朗日一欧拉法。文献[11-12]通过CFD模拟获得的振动刚性模型上的气动力和强迫振動位移时程，建立气动力离散时问气动模型，开展了颤振导数识别，但模型竖弯和扭转振动的幅值很小。文献开展了薄平板竖弯和扭转大幅振动CFD模拟和气动力小波分析，发现在气动力时程中出现了比薄平板振动频率高的二阶高频成分，也即大幅振动的薄平板气动力为非线性，但因CFD基于Fluent开展，而Fluent的动网格即使在均匀网格上也仅有一阶精度，因此研究结果的准确性需要验证。文献[14]基于Volterra理论开展了薄平板非线性气动力系统识别和CFD模拟，通过竖弯振动幅值为1%和5%薄平板宽度的薄平板强迫振动CFD分析，发现气动力的非线性效应并不显著，原因可能是薄平板的振动幅度偏小。下面就桥梁气弹特性分析常采用的薄平板为研究对象，针对二维不可压流动特征，首先推导能适应任意流变区域的流体运动控制方程，再根据刚性断面振动的特点，给出合适的空问点运动方式和相应的流体运动控制方程。并以薄平板为例研究其大幅振动的气动特征。

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

基于式（18）-（22）可在二维计算域上实现对任意流变区域内流体的运动描述，也即为描述任意流变区域不可压流动的基本方程，其中映射关系式（1）需根据所研究的流动区域的具体变形形式来确定。当物理域变形特征简单时，对应的映射关系可简化，此时，式（18）-（22）可能会得到明显简化，比如桥梁断面的气动弹性问题，因桥梁断面的绝对刚性假设，使得映射关系明显简化，如下文所述。

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

式（1）和（18）-（22）能描述任意变形物理域的流体流动，比如带移动边界的流体动力学问题。移动边界可为内流场或外流场流体外部边界的变形或移动，也可为外流场绕流物体在流体力作用下其形位的改变，这些均能导致与绕流物体接触的流体边界的移动和变形。在进行桥梁断面一类绕流场模拟时，由于作用在桥梁断面上的气动压力和摩阻力很小，因此桥梁断面的变形很小，可近似将桥梁断面看作刚体。这样，桥梁断面绕流场内边界的运动可用刚性桥梁断面的整体运动来描述，也即采用刚体运动的描述方法。在桥梁气动弹性分析中，通常用桥梁断面的竖弯和绕某一参考轴的转动来描述桥梁断面上任意点的运动。

对振动的刚性断面，如果其计算域的远场边界不变，将能恒定地给出远场边界条件，但由于计算域内、外边界之问的物理区域在不断发生变化，如果进行数值模拟，则需在每一时问步形成一次网格，如果计算的时问步数量很大，由网格生成所带来的附加计算量将很大，且难以保证每一步上重新生成的网格的高质量。相比较，如果在离刚性断面足够远处划定外边界，并对计算域划分高质量网格，数值模拟时将计算域网格和刚性断面一起同步振动，就可只在计算开始时形成一次网格，而在每一时问步上给出相应的边界条件，如此处理，不仅使得数值模拟的网格质量明显提高，也能省去每一时问步上网格的重新生成。这样，需求解的物理区域只随时问发生位置改变，形状并不改变。由于计算域事先给定，计算域和物理域的对应关系也能明确。

取绝对坐标系OXY平面内任意点为物理域之参考点，刚性断面的振动方式为相对于参考轴的平动和绕参考点的转动。如取计算域的坐标原点在物理域的映射点为参考点，并在此参考点上建立计算域的相对坐标系oxy，当相对坐标系随刚性断面作同步运动时，计算域Dcomp（x，y）和物理域Dphys（x，y，t）存在下述映射关系

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

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

3.1CFD数值实现

对图6所示的长度为26的无厚度扭转振动刚性薄平板（为便于查看，图示有厚度），其绕流控制方程式本文采用非定常不可压流计算的二阶Projec-tion算法，即将上述控制方程分裂为依次求解动量方程得到中问速度场，由中问速度场求解压力场的Poisson型方程，再由获得的压力场更新中问速度场，得到修正的速度场方程，如此反复求解，反复修正，并依据修正速度场在连续方程式（37）中产生的残余值的大小来判断当前时问步流场变量的收敛情况。当修正速度对连续方程的残余值满足收敛条件时，求解进入下一时问步。

由二阶Projection算法得到的偏微分方程，空问离散采用基于交错网格的有限差分法。在关于中问速度场的动量方程求解上，左端动量项采用二阶外插的Adams-Bashforth格式，右端黏性项采用二阶中心差分格式;压力场方程左端项采用二阶中心差分格式;时问导数采用二阶中心格式，得到的差分方程均采用三对角矩阵算法求解，并通过Fortran程序编程实现。

因刚性薄平板厚度为零，计算域网格采用直角坐标系下的结构化网格，物面上下第一个网格高度为0.001b，并沿垂直物面方向网格尺度逐渐增大，但增长率不大于1.05，以捕捉边界层的复杂流动，并适度控制网格规模和非定常计算量。薄平板沿流向均勻划分32个网格，但从薄平板前缘往计算域入口，以及从薄平板后缘往计算域出口，采用不大于1.05的网格增长率逐渐增大网格。计算域总网格数为168×168，图7是计算域整体和薄平板周围的网格划分，其中红色方框内有模拟的薄平板。

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

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

从控制方程式（34）-（36）和式（39），（40）的推导可知，网格系统和薄平板刚性固定，同步振动。因此在每一时问步，基于式（41）可得到网格和薄平板的同步振动速度;薄平板表面采用无滑移条件，因此薄平板物面流体的运动速度等于所在薄平板表面扭转速度;对计算域的人口和上、下边界，因来流平行于绝对参考系OXY的OX轴，需在每一时问步上分别投影到运动参考系oxy的两个坐标轴方向，成为计算域对应边界的速度边界条件;计算域出口设置为纽曼边界条件。对压力边界条件，在薄平板物面和计算域人口、上下边界也采用纽曼边界条件。式中B=26。考虑作用在振动薄平板上的气动力对雷诺数不敏感，从降低不可压流非定常迭代求解的计算量考虑，本文所有模拟的雷诺数均为1200。

3.2CFD方法验证

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

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

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

4.1振动能量的定义

在来流中绕扭转轴作扭转强迫振动的薄平板，定义薄平板表面的压力系数为

为对比薄平板大幅振动的气动力特征，下面给出了薄平板小幅扭转振动（振幅0.2°）的气动力时程。图9（a）和（b）分别是折算风速Vr=4和24时，一个周期内的气动力系数时程曲线。可见无论是升力还是扭矩系数时程，其曲线具有明显的谐波特性，其波形频率等于扭转强迫振动频率，可知当薄平板小幅扭转振动时，其气动力系统具有良好的线性特性;另外，从不同折算风速下的曲线来看，升力和扭矩问的相位差随折算风速不同而有变化。

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

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

4.3大幅扭转振动

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

图13为不同折算风速下大幅扭转振动的薄平板，在一个振动周期内气流在振动薄平板上表面各点所作的无量纲总功。与图11对比可见，在计算的所有折算风速范围上，转动轴到薄平板后缘气流做负功，因而是薄平板上耗散气流能量的稳定区，且在低折算风速（Vr=4-16）下，薄平板表面的能量消耗特征与小幅度扭转振动类似。但当折算风速高于20后，迎风侧的半个薄平板上每个点在一个周期内气流做的功均为正，即为不稳定区，气流对薄平板作正功，薄平板从气流中吸收能量，因此薄平板上的激励区与稳定区共存，而一个周期内气流对整个薄平板表面所有点的总功之和决定了一个周期内气流对薄平板所作的总功，该总功的正负反映了扭转振动薄平板的气动稳定性。

图14给出了一个扭转振动周期内气流对薄平板所作总功随折算风速的变化，可见当折算风速小于16时，一个周期内气流对薄平板所作总功为负，也即薄平板消耗气流的能量;但当折算风速大于20后，气流对薄平板所作的总功由负变为正，此时一个周期内扭转振动薄平板从气流中吸收能量，这表明大幅度扭转振动的薄平板已进入气动不稳定状态。

5结论

本文推导了二维流变区域的任意拉格朗日一欧拉描述方程，基于刚性薄平板简化得到了薄平板强迫振动绕流场控制方程，采用数值方法开展了薄平板大幅强迫扭转振动的气动特征研究，可得到下述结论：

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

2）采用有限差分法开展了薄平板扭转振动的CFD模拟。从小幅扭转与大幅扭转振动气动力时程的对比发现，大幅扭转振动气动力呈现非线性特征，折算风速越高，气动力的非线性特征越显著。

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