基于遗传混合算法的二维耦合颤振方法

郑史雄， 朱进波， 唐 煜， 郭俊峰

(1.西南交通大学土木工程学院，四川成都610031;2.西南石油大学土木工程与建筑学院，四川成都610500)

颤振理论认为桥梁颤振是一种气动弹性不稳定现象，当气流经过不规则桥梁断面时，流体与结构产生相互激励作用，随着来流风速的增大，气流对结构的输送能量趋近于结构振动所耗散的能量，表现为运动耦合方程中系统振动阻尼由正值逼近0，而当风速超过一定值后，系统振动阻尼变为负值，桥梁振动幅度加大并趋于运动发散，线性颤振发生.针对颤振现象的研究，Scanlan等［1-2］首先基于风洞模型试验，引入颤振导数的概念，描述非定常气动自激力，应用于振动微分方方程中，建立了经典的自激力框架模型［3-4］.在此模型基础上，许多研究者通过求解该框架模型，研究分析结构的颤振机理，先后提出了不同的二维颤振分析方法.

总体上，二维颤振分析方法可以分成两类:一类求解方法是二自由度复模态特征值解法［5-6］，其基本思想是将二维颤振问题转变为求解复特征值问题，以某阶模态下的阻尼比为0时，作为系统发散的依据，认为颤振发生，对应的风速为颤振临界风速.该方法通过设置风速级数量和控制复频率的收敛性误差，可以得到方程组的精确解.为深入了解颤振机理，分析颤振导数在颤振过程中的贡献，Matsumot等［7-8］对一系列简单二维断面进行了系统研究，提出了另一种求解方法，即颤振分步分析法(step by step analysis);项海帆和杨咏昕等［9-10］在此基础上提出了描述颤振自由度参与程度的合理方法，能够定量描述颤振形态.该方法以某一自由度牵连运动阻尼比为0作为发散的根据，基于不同自由度间的激励-反馈原理来解耦系统方程，能够研究桥梁断面的颤振驱动机理.

上述两类方法都是频率迭代方法，即在某级风速下，选取初始的频率值代入频率迭代方程，求得该级风速下的振动频率与系统的阻尼，然后又逐级变化风速搜索直到出现阻尼比为0时的风速与频率值.丁泉顺等［11］已经指出上述两类方法结果的一致性.近年来，研究者分别对两类方法进行了优化发展，如葛耀君等［12］提出的二维颤振直接分析法就较好地提高了传统复模态特征值解法的计算效率.

本文将遗传混合算法应用于二维二自由度耦合颤振分步分析中，提出了基于遗传混合算法的分析方法.该方法将二维颤振分析转变为求解非线性方程组的问题，引入具有自组织、自适应和自学性能的遗传算法［13］，不需要人为选取初始频率，自动搜索每级风速下全局最优的各自由度振动频率解，避免了迭代算法陷入局部收敛甚至不收敛的情况.同时，为了避免在搜索最优解时的误差，再引入良好的最优算法——L-M 算法［14］——进行局部收敛修正.算例的计算结果验证了该方法的可靠度和适用性.

1 二维二自由度耦合颤振分析方法

基于Scanlan的经典自激力框架模型，首先给出自激力作用下二维桥梁节段仅做竖向和扭转两个方向的颤振运动方程:

式中:h、α分别为桥梁断面在竖向和扭转方向的位移;mh和I分别为结构竖向广义质量和扭转广义质量惯矩;ξh0和ξα0分别为结构竖向和扭转两个自由度运动的结构阻尼比;ωh0和ωα0分别为结构竖向和扭转两个自由度运动的结构固有频率，其中，广义质量、广义质量惯矩、结构固有频率都为常数;ρ为空气密度;v为来流平均风速，一般可以视来流为均匀流，风速大小根据实际情况确定;B为桥梁实际断面宽度;i=1，2，…，4)分别为无量纲气动升力导数和气动力矩导数，被称之为颤振导数;K为无量纲的折减频率，K=Bω/v，ω为系统振动圆频率.

求解上述方程时，二维二自由度耦合颤振分步分析方法将方程解耦为扭转牵连运动方程和竖向牵连运动方程，其振动频率分别对应为系统扭转频率和系统竖向运动频率.此处扭转牵连运动不仅包含扭转自由度运动，还包括竖向自由度运动，竖向牵连运动也是如此.

当二维桥梁节段运动系统做单自由度扭转运动时，运动方程形式变得格外简单，如式(2).

式中:Mse(α，α)为扭转运动自身所产生的气动力矩.

式中:Lse(h，h)为竖向运动自身所产生的气动力;Lse(h，α)为由竖向和扭转运动间的激励-反馈效应引起的气动力.

因为颤振为自激发散振动，对于每个牵连运动方程，可以转换为自由运动的形式，故将式(3)和式(4)牵连运动方程等式右端移到等式左端，结合为自由运动的形式，如式(5).

至此，成功地将一个耦合的振动系统解耦，系统是否发散取决于某一牵连运动方程的阻尼比.通过频率与阻尼比迭代求解方程组(5).

系统竖向运动频率为

系统扭转牵连阻尼比为

系统竖向牵连阻尼比为

式(6)～(9)清楚地用颤振导数表达出了系统阻尼与系统刚度，且能够反映不同状态下系统阻尼所包含的运动系统自身阻尼和各分项气动阻尼的影响贡献.

无量纲系数为

不同耦合运动间的相位差角为

式(11)可描述为将系统运动分解为两种牵连主运动，当系统以某一单自由度运动为结构主运动时，其运动方程称为主运动方程，此时自激力框架下的另一方程称为约束耦合方程，系统的某一单自由度运动会产生耦合气动力，并以外荷载的形式通过约束耦合方程强制激起桥梁节段的另一自由度的运动，形成的耦合气动力再反馈作用给主运动系统，构成了各自由度运动之间的激励-反馈机制［15-16］.

在式(1)～(11)的基础上，通过逐级搜索的方式找寻颤振临界点，具体的做法如下:风速从初始风速开始搜索，每级风速下求解式(6)～(7)中的系统扭转牵连运动频率和竖向牵连运动频率，计算并记录系统的各自由度牵连阻尼比.随着风速的增加，系统某一自由度阻尼比由正转负时，即认为系统发生颤振，当前风速的大小即该系统发生颤振的临界风速.值得一提的是，上述方法计算每级风速下的系统运动频率是通过给定圆频率初值代入式(6)和式(7)中，以不动点的迭代形式求得.根据上述原理，将二维二自由度耦合颤振分析法程序化，其计算流程见图1.

图1 传统的分步分析法流程Fig.1 Flowchart of the traditional method

2 传统数值解法的收敛性问题

2.1 构建非线性方程组

采用上述方法分析时，每级风速下的各自由度振动频率都需要经过式(6)和式(7)来迭代求得，以此得到启发，将颤振问题的分析转变为求解下列二元非线性方程组的问题:

方程组(12)中颤振导数、各自由度相位差角和无量纲系数都是关于系统振动圆频率的函数，故可以表示为

式中:X=(ωh，ωα).

求解上述形式的非线性方程组时，大量数值迭代法被提出，其主要思想是通过给定的初值代入迭代公式中逐次逼近准确值.文献［17］中给出了很多求解迭代的方法，有基于不动点迭代的简单迭代法，即原始耦合分析方法中使用的迭代法，本文中称为原始方法，如式(6)和式(7)所示;有为了增加方程组求解收敛速度而提出的基于泰勒展开级数思想的牛顿法;同时也涌现出大量的改善方法，如修正牛顿法、拟牛顿法等.牛顿法及拟牛顿法等具有比简单迭代法更高的收敛效率，完全适用于本文的求解过程，此类方法统一称为传统方法.

2.2 传统数值解法的收敛性问题

使用传统方法解非线性方程组时，首先要弄清楚方程组的性质.如线性方程组Ax=b，在理论上只要A－1存在，则解存在且唯一.而对于非线性方程组而言，情况变得尤为复杂，其可以有多解、唯一解或无解，最为重要的是迭代的收敛性与所取的闭区域有关，即传统方法仅具有局部收敛性.

通过一个较简单的非线性方程组来说明上述问题.

方程组(14)的真解为 x*=(0.7，0.7)，若取初值 x0=(0.1，0.1)作为初始点代入牛顿法［18］中收敛到点(1.550 7，－0.737 8)，迭代结果如图 2 所示;若初始点取 x0=(0.9，0.9)或 x0=x*+0.5，则能够收敛到精确解，迭代结果如图3所示.

图2 初值不在局部收敛域中的牛顿迭代法结果Fig.2 Results of the Newton iteration method with initial value not in the local convergence domain

上述类似方法都可以解决二维耦合颤振分析方法中每级风速下频率的计算问题，不同迭代法的优劣性也会因Jacobian矩阵的计算和函数求值的复杂性而不同，但都仅具有局部收敛性，初值选取不利会导致最终的迭代结果远离真解的情况.而用选取不动点迭代方法时，将式(14)迭代形式变为

图3 初值在局部收敛域中的牛顿迭代法结果Fig.3 Results of the Newton iteration method with initial value in the local convergence domain

方程组(15)为Jacoobi迭代格式，选取多个初值都做发散变化不曾收敛，其收敛过程如图4所示.

图4 Jacoobi迭代格式迭代结果Fig.4 Iteration results of Jacoobi iterative format

方程组(14)也可以转变为

方程组(16)为Gauss-Seidel格式，选取真解附近的初始值(0.71，0.70)，迭代过程出现周期震荡性的不收敛，如图5所示.可见简单迭代法的收敛性受迭代形式与格式的制约.

图5 Gauss-Seidel格式迭代结果Fig.5 Iteration results of Gauss-Seidel iterative format

对于不复杂的非线性方程组，传统的数值解法计算局限性都已经很明显了，而本文中要求解的方程组(12)很不规则，并且参数随着风速变化，想要从理论上研究其在某个取值区域的收敛性质是十分困难的.严格来说，在方程性质尚不明确的条件下，采用局部收敛的迭代格式(传统方法)进行求解是存在一定风险的，可能并不能收敛到真解.为规避这一风险，研究全局收敛的数值求解方法显得很有必要.

3 遗传混合算法

结合遗传算法和传统方法的优点提出了遗传混合算法，规避了传统方法避免局部收敛甚至发散的情况.遗传算法具有群体搜索和全局收敛性，不需要人为选取初始点，具有自行适应性，克服了传统算法对初始点的敏感性问题.同时引入L-M算法对遗传算法计算结果下的优良个体进行局部搜索，克服遗传算法局部搜索缓慢的问题.

3.1 遗传算法

结合求解的非线性方程组，介绍遗传算法［19］的基本步骤.首先结合方程组(13)将求解方程组等价于求解如下的最小值优化问题:

式中:D为方程组解的区间.当f(ω)=0时，对应的ω即为方程组的解.通常利用目标函数来定义适应度函数，而适应度为非负值，适应度值越大对应的解越合理.所以对于求极小值问题，将目标函数变换为适应度函数，如式(18).

式中:C为正常数.

对上述优化问题实现遗传算法求解时，需要先进行解的编码与解码.由于实数编码容易过早收敛，本文采用具有稳定性高、种群多样性等优点的二进制编码，即将问题的解即系统扭转频率值分别空间映射到位串空间Bl上，其中，下标l为位数，Bl取0或1，每串编码上的每一位称为基因，所在位置称为基因位.在此基础上进行遗传算法操作，得到的解经过解码还原成实数解并进行适应度评价.

遗传算法的基本步骤及解释:

(1)随机产生N个个体构成二进制编码的初始群体，表示为 P(0)={ω1，ω2，…，ωN}，其中ωi={ωhi，ωαi}，i=1，2，…，N，ωi被形象地称为染色体，其上基因的组合构成染色体不同的表现型.

(2)根据目标函数f(ω)计算当前群体中各个体的适应度，即ffit(ωi).

(3)判断算法终止条件是否满足，若满足则转(8).

(4)根据各个体的适应度执行选择操作，本文采用基于各个体适应度概率决定它们继续繁衍还是消亡的轮盘式选择方法.

个体的选择概率为pi=ffit(ωi) ∑ffit(ωi)，，生成一个［0，1］内的随机数 r，若 p1+p2+…+pi－1＜r＜p1+p2+…+pi，则选择个体 i.

(5)按交叉概率pc执行交叉操作.将选择操作后的某两个父代个体的部分片段加以替换重组形成新的子代，在二进制编码中选择交叉点，体现为部分字符串的互换.

(6)按变异概率pm执行变异操作.以一定的概率选择父代个体的某一基因座，通过改变该基因座的基因变化为新个体，在二进制编码中体现为某段字符串的取反，即若原来为0，则变为1，若为1则取0.

(7)若已得到由N个新个体构成的新一代群体，则转(2)，否则转(4).

(8)输出搜索结果，终止.

综上，遗传算法［19］可以形式化描述为

式中:l为二进制编码长度，取决于解的取值范围和精度要求;

s为选择、交叉和变异的策略方法;

g为遗传算子，通常包括交叉算子和变异算子;

p为遗传算子的执行概率，通常包括交叉概率和变异概率;

f为适应度函数;

t为迭代次数即最大进化代数或者终止准则.

依照上述流程，根据相关经验选取种群规模、染色体长度、最大进化代数、交叉概率和变异概率等遗传参数，能够有效地找出非线性方程组(12)的全局最优解.

3.2 L-M 算法

遗传算法虽有很好的全局收敛性，但局部搜索能力较差，在接近准确解时，搜索进程变得缓慢甚至停止.本文引入基于约束最小二乘思想的L-M算法对遗传算法求出的近似解进行局部修正.首先与方程组(17)相同，确定求解目标函数:式中:F(ω)=(f1(ω)，f2(ω))T，ω=(ωh，ωα)T.

设已知第 k次的迭代点ωk，则由泰勒公式可知:Δ

fu(ω)≈fu(ωk)+(fu(ωk))T(ω－ωk)，u=1，2.所以 F(ω)≈F(ωk)+A(ωk)(ω－ωk).其中:

记 Ak=A(ωk)，则有

式中:dk=ω－ωk.

记 φ(ω)=［Akdk+F(ωk)］T［Akdk+F(ωk)］，那么求解φ(ω)的极小点则近似原问题的极小点.而为克服雅克比矩阵的奇异致线性搜索得不到进一步下降，故将目标函数(20)转为如式(23)的信赖域模型.

式中:hk为信赖域半径.

这个方程的解可由解如式(24)得到.

适当调整 λk，保持 AT(ωk)A(ωk)+λkI正定，从而

若 (AT(ωk)A(ωk))－1AT(ωk)F(ωk) ≤ hk，则 λk=0;否则 λk＞0.由于 AT(ωk)A(ωk)+λkI正定，从而式(25)产生的方向一直是下降方向，最终能够搜索到精确解.

该算法结合传统高斯-牛顿法和梯度下降法的优点，同时增加对角元改善雅克比矩阵奇异性，避免了传统的Gauss-Newton算法不迭代情况，其快速的局部搜索能力得以充分发挥，具有很强的适用性.

3.3 遗传混合算法的应用

遗传算法拥有强大的全局搜索能力，其理论比较成熟且得到了数学上的证明，已然逐步应用到工程实践领域，而L-M算法克服了雅克比矩阵奇异的问题且拥有较快的收敛速度，两种方法组合的遗传混合算法成为一种可行高效的方法.基于遗传混合算法的二维二自由度耦合分步分析方法的思路如下:

(1)以系统振动频率为未知数，根据式(10)和式(11)分别表示各个相位角与无量纲量.

(2)以折算频率表示各颤振导数.若风速已知，则颤振导数仅与系统振动频率有关.

(3)输入风速v.

(4)求解关于系统频率的二元非线性方程组(12).首先应用遗传算法全局搜索最优近似解，其次将近似解作为L-M算法迭代的初始点做局部精细搜索.

(5)将求解的系统频率精确解代入式(8)和式(9)中计算系统阻尼比.若阻尼比都大于0，则风速增加一个定值Δv，返回(3)，重新计算;若某一阻尼比等于0，则程序终止，对应风速为颤振临界风速.

4 数值算例

为验证上述颤振直接分析方法的可靠性和适用性，用该方法对一平板的耦合颤振问题进行分析，并将分析结果与传统不动点迭代分析方法的结果进行比较.

二自由度平板断面基本参数如下:竖向振动的固有频率为0.120 1 Hz，扭转振动的的固有频率为0.280 1 Hz，断面宽度 B=38.8 m，每米长度广义质量m=22 600 kg/m，广义质量惯性矩 Im=3 667 100 kg·m2/m，空气密度 ρ=1.225 kg/m3，结构竖向和扭转振动阻尼比均取为0.005，平板的颤振导数采用Theodorsen解［20］.分别用上述两种方法进行运算，颤振分析结果列于表1，随风速变化的系统振动频率与阻尼比如图6～7所示.

从图6～7和表1中数据可知，在取定的多个检测风速节点处，本文方法与传统方法的计算结果误差都低于0.1‰，结果几乎一致，对应的颤振临界风速完全一致.而本文建立的遗传混合算法，由于具有全局收敛的先天优势，故更适宜应用于大跨度桥梁的颤振分析中.

表1 桥梁断面颤振分析结果对比Tab.1 Comparison results of the bridge flutter analysis

图6 随风速变化的系统振动圆频率比较Fig.6 System vibration frequency with the changing wind speed

图7 随风速变化的系统牵连阻尼比比较Fig.7 System damping ratio with the changing wind speed

5 结 论

通过对传统的二维二自由度耦合分步分析方法中数学问题的研讨，发现应用其中的数值迭代方法存在局部收敛性的问题.

将二维二自由度耦合颤振分析转变为求解关于系统振动频率的非线性方程组问题，使得求解的方式多元化.引入了不需要自行选取初始值的遗传混合算法进行全局搜索和局部精细收敛运算，建立基于遗传混合算法的二维耦合颤振分步分析方法.

理论推导与算例分析均表明，本文所提出的基于遗传混合算法的二维耦合颤振分步分析方法具有精度高、无条件收敛的优点.

［1］ SCANLAN R H，TOMKO J J.Airfoil and bridge deck flutter derivatives［J］. Journal of Engineering Mechanics，1971，97(6):1717-1737.

［2］ HUSTON D R.Flutter derivatives from 14 generic deck sections［C］∥ Proc. ofASCE Structure Congress.Orlando:ASCE，1987:281-292.

［3］ 陈政清，于向东.大跨桥梁颤振自激力的强迫振动法研究［J］.土木工程学报，2002，35(5):34-41.CHEN Zhenqing，YU Xiangdong.A new method for measuring flutter self-excited forces of long-span bridges［J］.China Civil Engineering Journal， 2002，35(5):34-41.

［4］ 李志国，王骑，伍波，等.不同攻角下扁平箱梁颤振机理研究［J/OL］.西南交通大学学报，优先出版.(2017-09-27) ［2017-10-28］.http://kns.cnki.net/kcms/detail/51.1277.U.20170927.2111.034.html.［5］ WILDE K，FUJINO Y，MASUKAWA J.Time domain modeling of brige deck flutter［J］. Structural Engineering Earthquake Engineering，1996，13(2):93-104.

［6］ XIANG H，ZHANG R.On mechanism of flutter and unified flutter theory of bridges［M］∥Wind Engineering Into the 21st Century.Rotterdam: ［s.n.］，1999:1069-1074.

［7］ MATSUMOTO M，HAMASAKI H，YOSHIZUMI F.On flutter stability of decks for super long-span bridge［J］.Structural Engineering EarthquakeEngneering， 1997，14:185-200.

［8］ MATSUMOTO M. Flutter instability and recent development in stabilization of structures［J］.Journal of Wind Engineering and Industrial Aecodynamics，2007，95(9/10/11):888-907.

［9］ 杨詠昕.大跨度桥梁二维颤振机理及其应用研究［D］.上海:同济大学，2002.

［10］ 杨詠昕，葛耀君，项海帆.大跨度桥梁典型断面颤振机理［J］.同济大学学报，2006，34(4):455-460.YANG Yongxin，GE Yaojun，XIANG Haifan.Flutter mechanism of representative sections for long-span bridges［J］.JournalofTongjiUniversity， 2006，34(4):455-460.

［11］ 丁泉顺，朱乐东.桥梁主梁断面气动耦合颤振分析与颤振机理研究［J］.土木工程学报，2007，40(3):69-74.DING Quanshun， ZHU Ledong. Aerodynamically coupling flutter analysis and flutter mechanism for bridge deck sections［J］.ChinaCivilEngineering Journal，2007，40(3):69-74.

［12］ 邵亚会，葛耀君，柯世堂.超大跨度悬索桥二维颤振频域直接分析方法［J］.哈尔滨工业大学学报，2011，43(8):119-123.SHAO Yahui，GE Yaojun，KE Shitang.Straight forward method for two-dimensional flutter analysis of superlong-span suspension bridge in frequency domain［J］.Journal of Harbin Institute of Technology University，2011，43(8):119-123.

［13］ 罗亚中，袁端才，唐国金.求解非线性方程组的混合遗传算法［J］.计算力学学报，2005，22(1):109-114.LUO Yazhong，YUAN Duancai，TANG Guojin.Hybrid genetic algorithm for solving systems of nonlinear equations［J］. Chinese Journal of Computational Mechanics，2005，22(1):109-114.

［14］ 张鸿燕，耿征.Levenberg-Marquard算法的一种新解释［J］.计算机工程与应用，2009，45(19):5-9.ZHANG Hongyan，GENG Zheng.Novel interpretation for Levenberg-Marquardtalgorithm［J］.Computer Engineering and Applications，2009，45(19):5-9.

［15］ 项海帆，葛耀君，朱乐东，等.现代桥梁抗风理论与实践［M］.北京:人民交通出版社，2005:206-216.

［16］ 葛耀君.大跨度悬索桥抗风［M］.北京:人民交通出版社，2011:259-269.

［17］ 朱小飞.迭代法解非线性方程组［D］.合肥:合肥工业大学，2014.

［18］ 晁玉翠.求解非线性方程组的修正牛顿法［D］.哈尔滨:哈尔滨工业大学，2007.

［19］ 屈颖爽.基于拟牛顿法和遗传算法求解非线性方程组的混合算法［D］.长春:吉林大学，2005.

［20］ 陈政清.桥梁风工程［M］.北京:人民交通出版社，2005:67-73.