主梁断面气动自激力的全过程数值模拟方法

摘要

主梁断面气动自激力的全过程数值模拟方法是为了解决主梁断面气动自激力的时域化是桥梁抖振时域分析的难点，实现气动自激力时域化方法的可行性和有效性，在已有动力方程推导的基础上，创建了一种专用于模拟气动刚度和气动阻尼的单元——Aero－dyn18单元，并给出了单元的相关参数、特性以及构造示意图。通过在桥梁结构有限元初始模型的基础上添加该Aero－dyn18单元，并采用CFD技术确定出Aero－dyn18单元的相关参数，就可以实现对主梁断面气动自激力的全过程数值模拟。显然与现有气动自激力时域化方法相比，该方法简单实用，大大提高了分析的效率，同时节约了成本，便于广大工程人员中推广应用。

说明书

技术领域

本发明涉及一种主梁断面气动自激力的全过程数值模拟方法，尤其适用于在进行桥梁抖振响应时域分析时模拟主梁断面所受的气动自激力。

背景技术

随着桥梁结构跨径的不断增大，桥梁抖振分析显得越来越重要。抖振响应计算主要有频域和时域两类分析方法。早期国内外桥梁抖振分析主要在频域内进行，但由于频域方法在分析过程中只能计入一定数量的模态，得出的是结构响应值的统计特征，且只能进行线性分析，并不能够很好地应用于大跨度斜拉桥、悬索桥等非线性结构体系。

建立在数值积分基础上的大跨度桥梁抖振响应时域分析方法将激励转化为时间系列，通过动力有限元方法得到结构抖振响应的时程，可以非常方便地计入各类非线性因素的影响。由于考虑因素比较全面，时域分析方法是抖振计算的发展方向。但抖振响应时域分析计算量非常大，而且在三维随机脉动风场的模拟以及自激力的时域化离散处理等问题上仍有待进一步研究。

目前许多学者开始借助于现有大型通用有限元软件如ANSYS、ABAQUS、ADINA和MSC/NASTRAN等来进行大跨度桥梁的抖振时域分析，根据基于准定常气动力模型导出的主梁单元的气动刚度和气动阻尼矩阵，将自激力的影响以单元气动阻尼矩阵和单元气动刚度矩阵的形式加以考虑。这种方法充分利用了现有通用有限元分析软件强大的非线性时域分析功能，但其中存在两个主要问题：

(1)由于现有通用有限元软件不一定有专门的气动自激力单元，抖振分析时通常采用其它具有类似特性的矩阵单元来代替以进行自激力的模拟，导致计算工作量增加、出错概率加大、结果精度下降以及使用不方便等；

(2)主梁单元的气动刚度和气动阻尼矩阵中通常采用风洞试验中所得气动导数作为参数，不仅试验过程复杂，且花费较大、造成模型材料的浪费，现有计算流体力学(Computational Fluid Dynamics，CFD)技术在识别气动导数领域的最新研究成果没有得到充分应用。

发明内容

技术问题：本发明的目的是创建了一种用于主梁断面主梁断面气动自激力的全过程数值模拟方法，包括发明了一种专用于主梁断面气动自激力模拟的单元，并将其命名为Aero-dyn18单元，同时在确定单元参数时利用CFD技术在识别气动导数的最新研究成果，实现了气动自激力数值模拟全过程的自动化和时域化，方便快捷地完成考虑气动自激力影响的大跨度桥梁抖振响应时域分析。

技术方案：针对上述问题，本发明开发了一种专用于大跨度桥梁主梁断面气动自激力模拟的矩阵单元，准确便捷地实现了气动自激力的时域化。单元形式简单、各参数物理意义简单明了，便于在广大桥梁风工程技术人员中推广应用。为了实现气动自激力的时域化，根据基于准定常气动力模型和非线性动力学理论导出的与主梁所受气动自激力等价的单元气动刚度矩阵和气动阻尼矩阵，开发出了一种专用于气动自激力模拟的Aero-dyn18单元，并给出了单元的示意图、单元矩阵构造参数及其意义。解决上述问题所采用的技术方案流程如下：

第一步：根据结构设计图纸，采用可编程参数化设计语言建立大跨度桥梁的有限元计算初始模型；

第二步：基于计算流体力学技术识别出大跨度桥梁主梁断面的气动导数，并将其以可调用数组方式进行存储；

第三步：根据用于抖振分析的风速数据以及识别出的气动导数来确定所创建的气动自激力单元——Aero-dyn18单元中的参数；

第四步：将已确定参数的Aero-dyn18单元添加到第一步所得有限元计算初始模型中，获得已计入气动自激力的有限元计算模型；

第五步：在已计入气动自激力的有限元计算模型上添加静力风荷载及抖振力时程，进行抖振响应非线性时程分析。

第三步所述的气动自激力单元——Aero-dyn18单元可以通过输入9个实常数的方式来确定刚度矩阵K和阻尼矩阵C，以模拟主梁断面的气动刚度和气动阻尼矩阵；该单元具有两个节点，每个节点有6个自由度；Aero-dyn18单元没有固定的几何形状，在两节点i和j之间可任意变换外形；不对称情况下，Aero-dyn18单元矩阵为12×12矩阵，共有9个不同的系数：

$\left(\begin{array}{cccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& C_1& C_2& C_3& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& C_4& C_5& C_6& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& C_7& C_8& C_9& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& C_1& C_2& C_3& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& C_4& C_5& C_6& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& C_7& C_8& C_9& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$

第四步所述的将已确定参数的Aero-dyn18单元添加到第一步所得有限元计算初始模型的过程中，每个桥面主梁节点处需添加两个Aero-dyn18单元E，包括一个刚度单元和一个阻尼单元；单元E的一个节点为主梁节点i或j，另一个节点可固定在任意位置。

有益效果：目前，由于没有专门的气动自激力单元，抖振分析时通常采用其它具有类似特性的矩阵单元来代替以进行自激力的模拟，导致计算工作量增加、使用不便等一系列问题。该专利在对现有类似矩阵单元进行改进的基础上，创建了一种专用于自激力模拟的矩阵单元，较为简单便捷地实现了气动自激力的时域化模拟，解决了大跨度桥梁的风致抖振响应时域分析中的一个关键问题。同时在进行气动参数识别时，采用了先进的CFD技术而非耗费大量人力物力财力的风洞试验，实现了主梁断面气动自激力的全过程数值模拟。由于时域分析和CFD技术是桥梁风工程发展的方向，本文方法在保证计算精度的基础上，大大降低了大跨桥梁抖振响应时域分析的难度，提高了分析的效率，节约了宝贵的社会资源，因此在未来的桥梁风工程领域具有广泛的应用前景，必将产生显著的社会和经济效益。

附图说明

图1主梁断面气动自激力全过程数值模拟流程图

图2风荷载作用下桥梁主梁断面的气动自激力示意图，

图3Aero-dyn18单元示意图，

图4Aero-dyn18单元的系数矩阵，

图5采用Aero-dyn18单元进行气动自激力有限元模拟示意图。

具体实施方式

根据上述技术方案，在计算大跨度桥梁的非线性抖振响应时，主梁断面气动自激力的时域化模拟可通过在主梁单元两端添加Aero-dyn18单元来实现，其中Aero-dyn18单元所需气动参数采用CFD技术进行识别，该方法具体应用于大跨度桥梁抖振分析中的实现过程包括如下5个步骤：

1)根据设计图纸建立大跨度桥梁的初始有限元计算模型；

2)采用CFD技术识别主梁断面的气动导数；

3)根据风速数据和主梁断面的气动导数，按照式(3)、(4)和(5)来共同确定用于模拟刚度矩阵和阻尼矩阵的Aero-dyn18单元的参数；

4)在初始有限元模型的基础上添加Aero-dyn18单元，得到考虑了气动自激力的用于抖振分析的有限元模型；

5)在用于抖振分析的有限元模型基础上，添加静力风荷载和抖振力时程进行抖振响应分析求解。

具体如下：

根据Scanlan教授提出的气动自激力表达式，主梁单位长度上所受的气动升力L_se、气动阻力D_se和气动扭矩M_se可分别表示为竖向位移h、水平位移p和扭转位移α的函数，采用无量纲气动导数H_i^*、P_i^*、A_i^*(i＝1，2，…，6)来表达：

$L_{\mathrm{ae}}=\frac{1}{2}{\mathrm{\rho U}}^2\left(2B\right)[{\mathrm{KH}}_1^{*}\frac{\stackrel{·}{h}}{U}+{\mathrm{KH}}_2^{*}\frac{B\stackrel{·}{\alpha }}{U}+K^2{H}_3^{*}\alpha +K^2{H}_4^{*}\frac{h}{B}+{\mathrm{KH}}_5^{*}\frac{\stackrel{·}{p}}{U}+K^2{H}_6^{*}\frac{p}{B}]---\left(1a\right)$

$D_{\mathrm{ae}}=\frac{1}{2}{\mathrm{\rho U}}^2\left(2B\right)[{\mathrm{KP}}_1^{*}\frac{\stackrel{·}{p}}{U}+{\mathrm{KP}}_2^{*}\frac{B\stackrel{·}{\alpha }}{U}+K^2{P}_3^{*}\alpha +K^2{P}_4^{*}\frac{p}{B}+{\mathrm{KP}}_5^{*}\frac{\stackrel{·}{h}}{U}+K^2{P}_6^{*}\frac{h}{B}]---\left(1b\right)$

$M_{\mathrm{ae}}=\frac{1}{2}{\mathrm{\rho U}}^2\left(2B^2\right)[{\mathrm{KA}}_1^{*}\frac{\stackrel{·}{h}}{U}+{\mathrm{KA}}_2^{*}\frac{B\stackrel{·}{\alpha }}{U}+K^2{A}_3^{*}\alpha +K^2{A}_4^{*}\frac{h}{B}+{\mathrm{KA}}_5^{*}\frac{\stackrel{·}{p}}{U}+K^2{A}_6^{*}\frac{p}{U}]---\left(1c\right)$

式(1)中，ρ为空气密度；U为平均风速；B为桥面宽度；K＝Bω/U为无量纲频率；ω为振动圆频率；气动导数H_i^*、P_i^*、A_i^*(i＝1，2，…，6)是无量纲风速$\tilde{U}=U/\left(\mathrm{fB}\right)$或无量纲频率的函数，它们的值与桥梁截面的几何形状有关。L_se、D_se和M_se的方向及B、α等参数的示意见附图2。

显然，主梁单位长度上的受到的气动力为分布荷载，将这些分布荷载等效为作用于单元两个节点的集中荷载：

$\left(\begin{array}{c}{\left(F_{\mathrm{ae}}^e\right)}_i\\ {\left(F_{\mathrm{ae}}^e\right)}_j\end{array}\right)=\left(\begin{array}{cc}{\left(K_{\mathrm{ae}}^e\right)}_{\mathrm{ii}}& {\left(K_{\mathrm{ae}}^e\right)}_{\mathrm{ij}}\\ {\left(K_{\mathrm{ae}}^e\right)}_{\mathrm{ji}}& {\left(K_{\mathrm{ae}}^e\right)}_{\mathrm{jj}}\end{array}\right)\left(\begin{array}{c}{X}_i^e\\ X_j^e\end{array}\right)+\left(\begin{array}{cc}{\left(C_{\mathrm{ae}}^e\right)}_{\mathrm{ii}}& {\left(C_{\mathrm{ae}}^e\right)}_{\mathrm{ij}}\\ {\left(C_{\mathrm{ae}}^e\right)}_{\mathrm{ji}}& {\left(C_{\mathrm{ae}}^e\right)}_{\mathrm{jj}}\end{array}\right)\left(\begin{array}{c}{\stackrel{·}{X}}_i^e\\ {\stackrel{·}{X}}_j^e\end{array}\right)---\left(2\right)$

式(2)中，X_i^e和X_j^e分别为单元e在i和j节点的位移向量；和分别为单元e在i和j节点的速度向量；K_ae^e和C_ae^e分别为单元e的气动刚度和气动阻尼矩阵。当采用集中气动力矩阵时，气动刚度矩阵和气动阻尼矩阵分别为：

$K_{\mathrm{ae}}^e=\left(\begin{array}{cc}{K}_{\mathrm{ae}1}^e& 0\\ 0& K_{\mathrm{ae}1}^e\end{array}\right)---\left(3a\right)$

其中，$K_{\mathrm{ae}1}^e=\frac{{\mathrm{\rho U}}^2{K}^2{L}_e}{2}\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& P_6^{*}& P_4^{*}& {\mathrm{BP}}_3^{*}& 0& 0\\ 0& H_6^{*}& H_4^{*}& {\mathrm{BH}}_3^{*}& 0& 0\\ 0& {\mathrm{BA}}_6^{*}& {\mathrm{BA}}_4^{*}& B^2{A}_3^{*}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)---\left(3b\right)$

$C_{\mathrm{ae}}^e=\left(\begin{array}{cc}{C}_{\mathrm{ae}1}^e& 0\\ 0& C_{\mathrm{ae}1}^e\end{array}\right)---\left(4a\right)$

其中，$C_{\mathrm{ae}1}^e=\frac{\mathrm{\rho U}\mathrm{BK}{L}_e}{2}\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& P_5^{*}& P_1^{*}& {\mathrm{BP}}_2^{*}& 0& 0\\ 0& H_5^{*}& H_1^{*}& {\mathrm{BH}}_2^{*}& 0& 0\\ 0& {\mathrm{BA}}_5^{*}& {\mathrm{BA}}_1^{*}& B^2{A}_2^{*}& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)---\left(4b\right)$

式(3)和(4)中，L_e为主梁单元e的长度。需要指出的是，以上两式的形式与总体坐标系的方向有关。由以上两式可知，气动刚度和气动阻尼矩阵可由桥梁截面的几何形状完全确定，因此根据式(1)可知，当桥宽B以及全部气动导数均为已知时，单元e两个节点上所受到的自激力由各自节点的位移即可确定，但前提是必需创建可用于式(3)和(4)所示的气动刚度和气动阻尼矩阵模拟的单元。

本发明所创建的气动刚度和气动阻尼模拟专用单元——Aero-dyn18单元具有两个节点，每个节点有6个自由度，其示实意图见图3。图中，X、Y、Z轴和X_e、Y_e、Z_e轴分别表示单元e的整体坐标系和局部坐标系；i和j表示单元e的两个节点。该单元没有固定的几何形状，它可以通过实常数的方式输入刚度矩阵K和阻尼矩阵C，以模拟结构系统的气动刚度和气动阻尼矩阵。

Aero-dyn18单元矩阵为12×12矩阵，共有9个不同的系数，分别对应式(3)或(4)中的矩阵系数，其矩阵构造示意图见图4。将已定义好的气动刚度矩阵系数或气动阻尼矩阵系数代入Aero-dyn18单元矩阵，该单元就可以用来模拟桥面受到的气动自激力。

由于一个Aero-dyn18单元只能模拟气动刚度或者气动阻尼矩阵中的二者之一，而不能同时模拟这两者，因此在每个桥面主梁节点处添加两个Aero-dyn18单元E，其中包括一个刚度单元和一个阻尼单元，如图5所示。在图5中，单元E的一个节点为主梁节点i或j，另一个节点固定，由于单元长度可以任取，因此可任意设定固定节点的位置。例如，在节点i处，单元E1用于模拟气动刚度而单E3用于模拟主梁节点i处受到的等效气动阻尼，单元E1和E3共用两节点。单元E2和E4的功能与模拟方法与单元E1和E3相同。

由上述分析知，当模拟主梁的各BEAM4单元的长度相等时，用于模拟气动刚度和气动阻尼的Aero-dyn18矩阵可由式(5)得到：

$K^{E_1}=-{2K}_{\mathrm{ae}}^e,$$C^{E_3}=-{2C}_{\mathrm{ae}}^e$(5a)

$K^{E_2}=-{2K}_{\mathrm{ae}}^e,$$C^{E_4}=-{2C}_{\mathrm{ae}}^e$(5b)

式(5)中，和分别为单元E1和E2的刚度矩阵；和分别为单元E3和E4的阻尼矩阵；K_ae^e和C_ae^e分别为单元e的气动刚度矩阵和气动阻尼矩阵。因此，将式(5)的刚度矩阵和阻尼矩阵转化到总体坐标系下并组集成结构总气动力矩阵，并分别与原结构的刚度矩阵和阻尼矩阵进行叠加，就可得到用于抖振分析的有限元模型的系统刚度矩阵和阻尼矩阵。同时，为了避开现有分析方法中气动参数识别所必需进行的节段模型的风洞试验，气动刚度矩阵和气动阻尼矩阵中的气动参数采用CFD技术进行识别。