非线性子系统的大振幅时域自激力模型

刘十一，葛耀君

（1.土木工程防灾国家重点实验室（同济大学），200092上海；2.上海飞机设计研究院，201310上海）

非线性子系统的大振幅时域自激力模型

刘十一1，2，葛耀君1

（1.土木工程防灾国家重点实验室（同济大学），200092上海；2.上海飞机设计研究院，201310上海）

针对强风作用下的大跨度桥梁或高速飞行中的机翼可能出现大振幅颤振响应，而现有时域自激力模型又无法模拟大振幅下的非线性气动力问题.提出一种新的非定常大振幅自激气动力模型及其参数拟合方法，新模型通过附加非线性微分方程组及附加气动力自由度来模拟气动力记忆效应及振幅非线性特性.对于同一主梁断面，使用一组模型参数即可模拟不同折算风速和不同振幅下的自激力，模型参数可通过风洞试验或CFD数值模拟结果拟合获得.结果表明，新模型能再现自激力阻尼特性随折算风速和振幅的变化，通过单频振动拟合得到的新模型能再现多频振动下的非线性自激力时程.

颤振；大振幅自激力；非线性子系统；记忆效应；参数拟合

模拟大跨度桥梁在强风作用下的颤振过程或颤抖振响应时，需对主梁自激力建模.平板自激力理论解表明，自激力特性与折算风速（或频率）有关，即包含气动力记忆效应.Scanlan等最早提出了实用的钝体桥梁自激力经验模型［1］，只适用于模拟小振幅情况下的自激力.由于大振幅下的非线性气动力不满足频率叠加原理，必须在时域中计算.现有时域模型中，基于阶跃响应函数的时域卷积［2］和基于有理函数近似的状态空间方法［3-4］均可模拟非定常线性自激力.但目前尚无时域模型能完整再现非定常非线性自激力，一种广泛采用的非线性气动力模拟方法是将气动力分为低频和高频分量［5-6］，低频分量可采用非线性准定常模型计算，高频分量则采用非定常线性模型计算，最终的总体气动力为高频和低频分量的叠加；另一种改进的高、低频叠加法［7］采用了流变模型来模拟气动力的高频分量.以上方法虽能模拟高频分量的非定常特性，但无法模拟大振幅产生的非线性.Diana等［8］提出了基于非线性多项式的气动力模型，吴腾等［9］采用人工神经网络来模拟气动力的非线性效应，这两种方法均能模拟大振幅下的非线性效应，但未模拟气动力记忆效应，属准定常非线性模型.本文提出一种新的非定常非线性时域自激力模型.新模型中，记忆效应通过附加自由度来模拟［10］，振幅非线性效应则通过非线性的附加微分方程组（即子系统）及其他非线性表达式来模拟.新模型所使用的微分方程组和表达式均根据自激力特点而构造，其中包含待拟合的模型参数.对于特定断面，模型参数可通过风洞试验或CFD数值模拟拟合获得.本文提出的自激力模型参数与频率、振幅无关，能反应自激力特性随折算风速（或频率）的变化，能模拟自激力特性随振幅的非线性变化.

1 自激力模型表达式

本文将均匀来流下的主梁自激力分为3部分：1）静风力部分；2）准定常部分，用于模拟当前运动对自激力的影响；3）非定常部分，用于模拟气动力记忆效应，由非线性子系统生成.在新模型中，这3部分均为非线性.将来流风速用向量u表示为

主梁中心处的瞬时相对风攻角θ∗和瞬时相对风速大小u∗分别为

式中：x、y分别为主梁中心的横坐标和纵坐标；α为主梁扭转角，以逆时针为正.

新模型将θ∗、u∗及其对时间导数作为输入变量.为区分不同的运动产生的气动力效应，将和分解为

为简化气动力表达式，下文使用瞬时相对风速u∗、主梁宽度B和空气密度ρ对所有物理量进行量纲约化，约化后的物理量用下标～表示，约化方程式为

式中fH、fV和M分别为主梁局部坐标系中的阻力、升力和升力矩.

将主梁局部坐标系中的无量纲气动力floc～表示为

新的自激力模型可表示为

式中：φ为附加气动力自由度；fst为约化静风力；fqs为约化准定常自激力部分；fus为约化非定常自激力部分，与φ有关；g=0为附加非线性微分方程组，即非线性子系统，它定义了φ的演化规律.由式（9）计算可得局部坐标系下的约化气动力，经坐标变换后可转化为全局坐标系下的气动力，即

式（9）中的fst、fqs、fus均为θ∗的非线性表达式，g同时为θ∗和φ的非线性表达式，静风力部分fst为

式中a0～a3为与频率无关的参数向量.

准定常自激力部分fqs为

式中：rαθ∗( )、rmθ∗( )、rumθ∗( )为瞬时相对攻角θ∗的多项式向量，bi、ci、ei（i=0，..，3）为与频率无关的参数向量.

非定常自激力部分fus为

式中：Rφθ∗( )为多项式矩阵，Ei（i=0，…，3）为与频率无关的参数矩阵.

附加非线性微分方程组g=0.

式中：φ3=；Kφ(θ∗)、Kφ3(θ∗)为多项式对角矩阵；kαθ∗( )、kmθ∗( )、kumθ∗( )为多项式向量；Fi、Gi（i=0，2，4）为与频率无关的参数对角矩阵；hi、li、pi（i=0…3）为与频率无关的参数向量.

为保证微分方程组（16）中的气动力自由度φ稳定，式（17）中的Kφθ∗( )和Kφ3θ∗( )只包含θ∗的偶数次项，且这两个多项式对角矩阵中的元素都大于0.同样，微分方程组（16）中只包含了φ的奇数次项.这样就能保证非线性子系统的自由度φ始终是稳定的.

将式（11）、（12）、（14）、（16）代入方程组（9）即可获得非线性子系统自激力模型的完整表达式.式（11）、（13）、（15）、（17）中的所有参数均可通过风洞试验或CFD数值模拟结果拟合获得.模型参数数量取决于所使用的附加微分方程数量.

2 模型参数拟合方法

根据已知的主梁位移时程和气动力时程，可拟合获得所有模型参数.拟合的基本流程为：1）假设一组模型参数的初始值.2）根据已知的主梁位移时程和当前的模型参数值，采用数值方法求解附加非线性微分方程组，然后计算出气动力时程.3）根据计算的气动力时程和已知的气动力时程之差，采用一种数值优化算法计算模型参数的修正量.4）修正模型参数，如果参数收敛，则结束，否则转到第2步.

此方法中最关键的是第3步，即如何根据计算的气动力时程和已知气动力时程之差确定模型参数的修正量.由于模型参数数量大（大于100个），且每个迭代步计算量较大（需要求解非线性微分方程组），因此必须使用收敛较快的非线性数值优化算法.本文使用LM算法计算每个迭代步的参数修正量.LM算法是收敛最快的数值优化算法之一，广泛用于神经网络训练等复杂的非线性优化问题.

包含附加微分方程组的气动力模型可表示为

式中：f t()为模型输出的气动力；x t()为输入的主梁位移；向量pf、pg分别表示f、g包含的模型参数.

定义模型拟合的目标函数为

指定pf、pg的初始值后，本文使用LM算法优化pf、pg以使目标函数Π取得局部最小值.LM算法需要各离散时刻的气动力计算值f对模型参数pf、pg的偏导数矩阵.而f与气动力自由度φ有关，故需要各离散时刻的φ对模型参数pg的偏导数矩阵.φ为附加非线性微分方程组的解，由于φ没有解析表达式，本文采用四阶龙格库塔法计算φ的离散时程，而φ对参数的偏导数矩阵则采用链式法则在数值积分过程中同时计算.计算偏导数的基本思路：因为φ的数值积分过程由基本代数操作组成，而各代数操作的偏导数可直接计算，所以积分获得的φ的离散时程点的偏导数可通过链式法则间接计算获得.四阶龙格库塔法的1个积分步由4个子步组成，子步的偏导数计算流程如图1所示.

图1 一个数值积分子步中的偏导数计算流程

需要指出，本文给出的模型拟合算法不仅适用于新的非线性自激力模型，还适用于传统的线性状态空间模型，是一种通用的模型拟合算法.由于新的拟合算法寻求的是时域最小二乘解，实际上对原始信号中的高频随机噪声起到了时域滤波的作用，因此对原始信号的噪声有较强的容忍能力.此算法只能获得目标函数的局部最小值，因此拟合结果与参数初始值有关.当参数初始值在合理范围内时，一般都能获得理想的拟合结果，对式（11）、（13）、（15）、（17）中的模型参数，本文推荐初始值为

F0∈ 0.2，1[ ]，其他 =0.

3 参数拟合与模型验证

本文对一个均匀来流下的箱型主梁断面进行了CFD强迫振动绕流模拟，并通过表面压强积分计算出三分力时程.主梁断面形状及尺寸如图2所示. CFD模拟方法为LES，采用Smagorinsky静态湍流模型，雷诺数等于10 000（相对主梁宽度）.CFD模拟了不同攻角下的静风力时程，及不同折算风速、不同振幅下的强迫振动自激力时程.静风力包含11个工况，攻角为-10°～10°，步长为2°.单频扭转强迫振动包含24个工况，折算风速为3、5、7、9、11、13，振幅为±1°、±3°、±6°、±9°.单频竖弯强迫振动包含4个工况，折算风速为3、5、7、11，振幅为±0.2H（H为主梁高度）.此外还模拟了两个多频扭转强迫振动工况，用于模型验证.本文折算风速定义为Ur=U／（Bf）=2πU／（Bω）.本文用CFD模拟获得的11组静风力时程和28组单频强迫振动时程拟合出新模型的一组模型参数.拟合的模型中包含4个附加非线性微分方程组，即有4个附加气动力自由度，模型共包含140个参数.模型拟合经过约200次迭代达到收敛.图3～5给出了3种单频强迫振动下的阻力、升力和升力矩时程，图中包含CFD计算及新自激力模型拟合结果.需要指出，CFD计算获得的气动力时程包含主梁特征紊流产生的高频分量，无法在自激力模型中模拟.由图3～5可见，大振幅正弦强迫振动产生的自激力呈现明显的非正弦形状，验证了气动力非线性的存在.为进一步分析自激力的非线性特性，本文定义等效扭转气动阻尼c-α为

由此可计算出任意一组扭转强迫振动位移和气动力时程的等效扭转气动阻尼.若自激力为线性，等效扭转阻尼应与振幅无关.图6给出了CFD模拟得到的等效扭转气动阻尼随折算风速和主梁振幅的变化情况，由图可见折算风速及振幅对气动阻尼均有很强影响.对于当前断面，低折算风速下的气动阻尼随振幅增加而增加，但高折算风速下的气动阻尼随振幅增加而减小.此外，扭转振幅为1°和3°时气动阻尼几乎相等，由此可见，当振幅小于3°时，当前断面自激力呈线性.

图2 箱型主梁断面形状及尺寸

图3 扭转振幅±9°，折算风速13模型拟合结果

图4 扭转振幅±6°，折算风速5时模型拟合结果

图5 竖弯振幅±0.2H，折算风速3时，模型拟合结果

图7对比了多组等效扭转气动阻尼，分别由CFD模拟、拟合的新模型和拟合的线性状态空间模型计算获得.线性状态空间模型使用Minimum-state格式［11］，使用CFD模拟的单频强迫振动结果拟合.图中包含不同振幅（3°、6°、9°）下的扭转阻尼随折算风速的变化曲线.由于线性状态空间模型不能反映振幅对气动阻尼的影响，因此只有一条曲线.拟合的新模型能较好地反映扭转气动阻尼随折算风速和振幅的变化，模拟了气动阻尼的非定常和非线性特性.

图6 折算风速、扭转振幅对无量纲等效扭转气动阻尼的影响

本文用拟合的新模型计算了两个多频扭转强迫振动下的自激力时程，并将其与CFD模拟结果对比，如图8、9所示.由图可见拟合的新模型能准确再现多频大振幅强迫振动产生的自激力时程.因此使用单频强迫振动试验拟合非线性自激力模型是可行的.

图7 不同折算风速、不同振幅下等效扭转气动阻尼的拟合结果比较

图8 包含折算风速分布为7、11，振幅为±3°的扭转频率成分的多频预测结果

图9 包含折算风速分别为5、7、11，振幅为±3°的扭转频率成分的多频预测结果

4 结 论

1）提出一种新的非定常非线性时域自激力模型及其拟合算法，新模型中，气动力记忆效应由附加气动力自由度模拟，大振幅非线性效应由与瞬时相对风攻角和气动力自由度有关的非线性表达式模拟.新模型包含一组与频率、振幅无关的参数，这些参数可通过已知的位移和气动力时程拟合获得，拟合过程基于LM数值优化算法，并使用链式法则计算微分方程数值解的偏导数矩阵.

2）使用CFD强迫振动模拟获得了一个桥梁断面在不同折算风速和不同振幅下的气动力时程，并拟合了新模型的一组参数.针对CFD模拟结果，新模型能较好地拟合不同折算风速和振幅下的非线性气动力时程，并能准确反映扭转气动阻尼随折算风速和振幅的变化.使用单频强迫振动拟合的新模型能准确再现多频强迫振动产生的气动力时程.

3）由于CFD模拟与真实气动力之间存在一定差别，新模型还有待通过风洞试验进一步校核.本研究仅限于均匀来流下的自激力非线性问题，真实桥梁风振还受来流紊流产生的抖振力影响，抖振力非线性及其与自激力的相互影响还有待进一步研究.

［1］GE Y J，TANAKA H.Aerodynamic flutter analysis of cable-supported bridges by multi-modeand full-mode approaches［J］.JournalofWind Engineering and Industrial Aerodynamics，2000，86：123-153.

［2］COSTA C，BORRI C，FLAMAND O，et al.Time-domain buffeting simulations for wind-bridge interaction ［J］. Journal of Wind Engineering and Industrial Aerodynamics，2007，95：991-1006.

［3］CHEN X，KAREEM A.Aeroelastic analysis of bridges：effects of turbulence and aerodynamic nonlinearities［J］. Journal of Engineering Mechanics，2003，129：885-895.

［4］CHEN X Z，MATSUMOTO M，KAREEM A.Time domain flutter and buffeting response analysis of bridges［J］. Journal of Engineering Mechanics-Asce，2000，126：7-16.

［5］CHEN X Z，KAREEM A.Nonlinear response analysis of long-span bridges under turbulent winds［J］.Journal of Wind Engineering and Industrial Aerodynamics，2001，89：1335-1350.

［6］DIANA G，FALCO M，BRUNI S，et al.Comparisons between wind tunnel tests on a full aeroelastic model of the proposed bridge over Stretto di Messina and numerical results［J］.Journal of Wind Engineering and Industrial Aerodynamics，1995，54-55：101-113.

［7］DIANA G，ROCCHI D，ARGENTINI T.An experimental validation of a band superposition model of the aerodynamic forces acting on multi-box deck sections［J］. Journal of Wind Engineering and Industrial Aerodynamics，2013，113：40-58.

［8］DIANA G，ROCCHI D，ARGENTINI T，et al.Aerodynamic instability of a bridge deck section model： linear and nonlinear approach to force modeling［J］.Journal of Wind Engineering and Industrial Aerodynamics，2010，98：363-374.

［9］WU T，KAREEM A.Modeling hysteretic nonlinear behavior of bridge aerodynamics via cellular automata nested neural network［J］.Journal of Wind Engineering and Industrial Aerodynamics，2011，99：378-388.

［10］LIU S，GE Y.A new unsteady aerodynamic model with unique parameters［C］／／Proceedings of BBAA7.Shanghai：China Communication Press，2012.

［11］WILDE K，OMENZETTER P，FUJINO Y.Suppression of bridge flutter by active deck-flaps control system ［J］. Journal of Engineering Mechanics，2001，127：80-89.

（编辑 魏希柱）

Nonlinear dynamic subsystem model for large-amplitude motion-induced aerodynamic forces of bridge decks

LIU Shiyi1，2，GE Yaojun1

（1.State Key Laboratory of Disaster Reduction in Civil Engineering（Tongji University），200092 Shanghai，China；2.Shanghai Aircraft Design and Research Institute，201310 Shanghai，China）

In simulation of long-span bridge and airfoil flutter，a large-amplitude motion-induced force model is required.However，existing time-domain motion-induced force models cannot simulate memory effects and amplitude dependency simultaneously.This paper proposes a new time-domain non-linear motion-induced force model，which employs a set of nonlinear differential equations and augmented aerodynamic degrees of freedoms to simulate memory effects and amplitude dependency.Motion-induced forces under different deck motion amplitudes，different reduced wind velocities and non-sinusoidal motions can be simulated using a single set of model parameters.A fitting algorithm for the new model is also proposed and the model parameters of a bridge deck are fitted using results obtained by CFD simulations.Numerical example shows that the new model can reproduce variation of aerodynamic damping characteristics with reduced wind velocity and amplitude.A model fitted using single-frequency forced vibration tests can reproduce forces generated by multi-frequency vibrations.

flutter；large-amplitude motion-induced force；nonlinear dynamic subsystem；memory effects；model fitting

U441.3

A

0367-6234（2015）09-0073-06

10.11918／j.issn.0367-6234.2015.09.014

2014-01-20.

国家自然科学基金（91215302）.

刘十一（1986—），男，工程师，博士研究生；葛耀君（1958—），男，教授，博士生导师.

刘十一，liushiyi_ha＠163.com.