本发明涉及双圆柱建筑模拟
技术领域：
，具体涉及一种串列双圆柱尾流驰振非定常气动力数学模型建立方法。
背景技术：
：实际工程中建筑物通常是柱群结构，如海上结构、海底管道、大群烟囱、热交换管、桥墩、拉索、桅杆、化学反应塔、临近的摩天大厦和输电导线。流体力是设计这些工程结构的重要因素，因此理解临近柱群结构对流致振动的影响有重要的工程实际意义。双圆柱由于形式简单，是研究圆柱阵列流致振动的理想模型。当两个圆柱相互靠近时，上游圆柱产生的旋涡脱落会诱发第二个圆柱振动。双圆柱的布置形式及其间距是影响尾流驰振的重要因素。在现有技术中存在以下设计和实验：通过风洞实验研究表明，尾流驰振发展阶段的振幅特点为随折算风速的增加而接近线性增加。对等直径双圆柱的尾流干扰研究中，两圆柱之间的距离对尾流驰振的诱发有着较大的决定作用，流经上游圆柱的尾流区会对下游导线形成三个区域，分别为近距离失稳区、稳定区和远距离失稳区，尾流驰振只会在特定的间距区间才会发生。以等直径双圆柱间的流体形态与相对位置的不同为依据，将尾流迹对下游圆柱的影响分为尾流干扰和近距离干扰两种类型，并将l/d＝3.5和t/d＝2作为两个关键距离，其中，l为两圆柱在流向的垂直距离，t为两圆柱在横流向的垂直距离，d为圆柱直径。通过风洞实验结果对并列等直径导线的尾流区域进行了详细划分，近距离干扰区对应的导线间距比为1～1.1，尾流驰振现象发生的导线间距比为1.1～3.8，导线间距比超过3.8的为远距离干扰区。对处于固定圆柱后的弹性支撑的下游圆柱(仅允许横向振动)进行实验，在串列距离l/d＝3-5.6观察到振幅随折算风速不断增加的尾流驰振现象。有研究人员利用风洞对串列布置间距为l/d＝4.3-8.7的等直径双圆柱进行全尺度实验，在l/d＝4.3时观察到下游圆柱的尾流驰振现象，该现象于l/d＝6.5时消失。同时还有研究人员对不同间距比的双圆柱进行风洞试验，根据上下游圆柱是否发生驰振与涡激共振，划分了7个流致振动响应模式。尾流振动作为工程中常见的复杂振动形式，除了间距的影响，许多研究者对它的气动特性和振动响应有影响的其他因素进行了研究分析。例如通过试验研究了尾流下弹性支撑的圆柱吊索的干扰驰振，发现尾流驰振的临界风速随着sc数的增大而增大。例如系统地研究了输电线的结构阻尼和气动阻尼对尾流振动响应的影响。例如通过风洞试验对下游拉索的尾流驰振进行了研究，结果表明尾流作用下的下游拉索振动响应会受雷诺数和风速的影响。这些研究成果对揭露尾流驰振的特性具有一定促进作用，但是未深入揭露尾流驰振形成机理与诱发的原因。传统尾流试验多关注结构响应、表面压力分布和流场形态等特性，对非定常驰振力非线性特性等研究较少。现有技术中，最早的设计是基于准定常假设利用线性分析方法判断了输电线发生尾流驰振失稳的可能性。其次是基于波传递理论采用传递矩阵方法对分裂导线尾流驰振的振动特性开展了研究，对尾流驰振振动响应进行预测。基于上述两项研究成果，有研究人员完善了尾流驰振的准定常力学模型，建立了尾流驰振的非线性振动方程。采用风洞试验方法，测得了尾流输电线平均气动力系数的空间分布规律，并采用runge-kutta法求解，得到了尾流驰振的幅值。上述尾流驰振模型虽然都具有一定的合理性，但是建立在准定常理论基础之上，气动力模型本身仍是定常的，具有一定的局限性。因此，需要基于非定常气动力时程数据建立非定常气动力非线性数学模型。技术实现要素：有鉴于此，本发明提供了一种串列双圆柱尾流驰振非定常气动力数学模型建立方法，建立串列双圆柱尾流驰振模型，利于现有工程建设，提高工程建造稳定性和可靠性。其技术方案如下：一种串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于，具体步骤为：s1：采用计算流体动力学软件，模拟圆柱流向间距l/d＝4的情况下，上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，得到串列双圆柱的尾流驰振几何模型并对该几何模型进行网格划分；其中，l为两圆柱在流向的垂直距离，d为圆柱直径；s2：根据串列双圆柱的尾流驰振几何模型，设立至少一种工况，分析串列双圆柱的尾流驰振气动力的频谱规律；s3：根据步骤s1中的上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，建立含有未知参数的串列双圆柱的尾流驰振气动力数学模型；其中气动力表示为自激力和涡激力两种形式；s4：采用基于能量等效的最小二乘参数识别方法，拟合出自激力中气动阻尼力参数和涡激力参数；，得到最终的串列双圆柱的尾流驰振气动力数学模型；s5：对步骤s4得到的串列双圆柱的尾流驰振气动力数学模型进行对比。通过构建上述步骤，在l/d＝4最易发生尾流驰振的情况下，采用计算流体动力学(cfd)模拟了双圆柱，且双圆柱风的流向上，模拟了前柱固定，后柱横向振动串列双圆柱的尾流驰振。给出了气动力模型的数学参数表达式，求取未知参数后，得到最终串列双圆柱的尾流驰振气动力数学模型。该模型经验证和比对后，合理性好，精确度高，有利于以后建筑过程中更为稳定的修建工程。再进一步的技术方案，在步骤s1中，所述串列双圆柱的尾流驰振中，下游圆柱横向振动采用尾流圆柱横向振动方程求解，该下游圆柱横向振动方程为：其中，m为下游圆柱质量；y、和分别为垂直于流向的位移、速度和加速度；c和k分为为结构阻尼和刚度常数；fy为圆柱受到的垂直方向的力。本发明中，模拟工况设置时，二维圆柱绕流的非定常不可压缩rans方程如下：其中ui和p分别表示平均速度和平均压力；ρ为流体密度；μ为分子黏度；sij为平均应力张量；为雷诺应力张量。模拟工况设置时使用了剪切应力传递模型，即sst-kω模型。模型涉及两个运输方程，即湍动能和耗散率方程，其中采用了非稳态分离算法(usa)。关于模型具体详见文章：menterfr.two-equationeddy-viscosityturbulencemodelsforengineeringapplications[j].aiaajournal.1994,32.耦合速度表达式时用simple算法处理，并且针对不稳定情况采用了二阶隐式方案。二阶隐式方案用于k-ω传递方程和动量方程中的对流项。再进一步的，步骤s1中，所述串列双圆柱的尾流驰振几何模型中，采用的矩形计算域为60d×40d，d为圆柱直径；矩形计算域中，上游圆柱中心距离入口边界为25d，下游圆柱中心距离出口边界为31d；矩形计算域中，入口边界设为速度入口边界条件；矩形计算域中，流域出口边界设为压力出口边界条件；矩形计算域中，上下边界定义为对称边界条件；矩形计算域中，上游圆柱和下游圆柱均采用圆柱且圆柱表面采用无滑移壁面。再进一步的，步骤s1中，所述对几何模型进行网格划分时采用非均匀四边形结构化网格，上、下游圆柱表面的表面网格采用o型网格；由于壁面附近的流场变化剧烈，所述对几何模型进行网格划分时加密圆柱近壁面的网格以及控制壁面第一层网格高度，使上、下游圆柱壁面第一层网格高度y ＜1；来降低网格对数值计算的影响。为了模拟圆柱非稳态振动，矩形计算域中，入口下游27d和34d处的交界面采用了滑移网格技术，在横向边界处使用动态层铺法更新网格。再进一步的，在步骤s3中，涡振和驰振力都以自激力为主，涡激力仅在振动发展阶段起作用。而尾流驰振力中两个频率的涡激力在振动各个阶段一直存在，且成分占比很大，必须考虑。根据步骤s1中提及的上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振中，下游圆柱横向振动振动方程可得到上述带未知自激力参数和未知涡激力参数的串列双圆柱的尾流驰振气动力数学模型。含有未知自激力参数和未知涡激力参数的串列双圆柱的尾流驰振气动力数学模型为：其中，fy为结构的总气动力；为自激力的表达式，其中，自激力是自激振动导致的，自激振动是由于系统的力与运动状态的相互反馈的结果，因此通常将这种带有自激性的振动力可以描述为结构速度与位移的函数；其中，fv(fv)为涡激力的表达式；涡激力是上游圆柱的涡脱或者下游圆柱自身涡脱导致的，而涡脱频率就是漩涡脱落的表征，因此通常将这种由涡脱引起的表面压力可以描述为涡脱频率的函数。分别对y和进行高阶泰勒序列展开，可得到尾流驰振力中的自激力多项式表达式如下：式中gij为气动参数。g00这一项与描述响应状态的y和都没有关系，可以将该项认定为静力项，它对尾流驰振动力响应没有影响，故忽略该项。由之前的气动力频谱可以看出，对整个频率成分贡献的较大的基本在5阶倍频之内，故可以将泰勒展开式中超过y5和的项截断。基于能量等效原理，由三部分组成：气动阻尼力项(改变了结构的阻尼)、气动刚度力项(改变了结构的刚度)和纯力项(既不改变结构的阻尼也不改变圆柱的刚度)。则，自激力可以写成下式忽略纯力项，无量纲自激力表达式如下对于横向振动，气动刚度项仅轻微改变结构频率，可以不考虑，则上述公式可精简为：ρ为流体密度，u为来流的速度，d为圆柱直径，y为垂直于流向的位移，为垂直于流向的速度；α01为一阶无量纲气动阻尼参数，α03为三阶无量纲气动阻尼参数，α05为五阶无量纲气动阻尼参数，均为未知参数；首先，最大的有两个涡脱频率fv1和fv2，且由于涡激力始终是与涡脱频率息息相关的，故可先将其表示为涡脱频率的函数，则涡激力的表达式可变为下式：由于上式中的相位角与都是需要拟合求解的参数，无量纲后可以改写为下式：式中：δ11,δ12,δ21,δ22为无量纲后的涡激力参数；st1为上游圆柱斯托哈尔数；st2为下游圆柱斯托哈尔数；ρ为流体密度；ur为折算风速，fs为固有频率；v11v12v21v22分别表示无量纲前的涡激力参数；故，只考虑稳定振幅的尾流驰振力简化模型可表示如下：再进一步的，步骤s5的具体内容为：s51：根据得到的串列双圆柱的尾流驰振气动力数学模型，采用runge-kutta法反算位移和速度；s52：定义下游圆柱能量的参数名称：ωs为下游圆柱的自然圆频率；ws(t)为下游圆柱自身结构阻尼在当前时间的累积功，消耗能量为负值；wa,v(t)为利用拟合的参数反算的下游圆柱的气动阻尼力和涡激力在当前时间的累积功，提供能量为正值，ee(t)为下游圆柱在当前时间的机械能，其中包括动能与势能；s53：将步骤s51得到的位移、速度带入公式(1)-(3)中，重构得到气动阻尼力和涡激力的累积功；s54：将气动阻尼力和涡激力的累积功与所述串列双圆柱的尾流驰振几何模型的模拟值进行对比。与现有技术相比，本发明的有益效果：建立了非定常气动力的串列双圆柱尾流驰振数学模型，适用双圆柱建筑尾柱驰振，适用范围广，为以后含有双圆柱的建筑力学分析提供了有利的数据支撑，并且该模型经分析对比后，可反推演练，和其他被认可的现有技术相互支撑。确定了其正确性和可靠性。附图说明图1为计算域和边界条件示意图；图2为计算域网格和圆柱壁面网格划分结果示意图；图3为网格依赖性研究示意图；图4为三种模拟工况的无量纲振幅随折算风速变化结果示意图；图5为尾流驰振发展阶段ur＝9.179对应全时程的升力频谱变化示意图；图6为尾流驰振发展阶段ur＝14.49对应全时程的升力频谱变化示意图示意图；图7为工况一不同折算风速的升力频谱示意图；图8为上下游圆柱的st数随ur的变化规律示意图；图9为在ur＝13.889气动参数重构的升力做功以及机械能时程示意图；图10为在ur＝13.889气动参数重构的位移时程效果示意图；图11为在ur＝12.08涡激力和气动阻尼力参数重构分别做功曲线以及重构力曲线示意图；图12为在ur＝12.08各项气动力对位移的作用效果示意图图13为在ur＝12.08参数重构的自激力和涡激力时程示意图；图14为在ur＝12.08参数重构的升力与模拟值对比示意图；图15为基于升力时程拟合得到工况一的无量纲涡激力幅值参数示意图；图16为气动力参数识别结果示意图；图17为气动力模型预测振幅与模拟值的对比示意图；图18为串列双圆柱尾流驰振非定常气动力数学模型建立步骤流程图；图19为串列双圆柱的尾流驰振气动力数学模型进行对比步骤流程图。具体实施方式以下结合实施例和附图对本发明作进一步说明。参见图1-12提供串列双圆柱尾流驰振非定常气动力数学模型建立方法，具体步骤为：s1：采用计算流体动力学软件，模拟圆柱流向间距l/d＝4的情况下，上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，得到串列双圆柱的尾流驰振几何模型并对该几何模型进行网格划分；其中，l为两圆柱在流向的垂直距离，d为圆柱直径；l/d＝4时，易发生尾流驰振；在本实施例中，二维圆柱绕流的非定常不可压缩rans方程如下：其中ui和p分别表示平均速度和平均压力，ρ为流体密度，μ为分子黏度，sij为平均应力张量，为雷诺应力张量。本发明使用了剪切应力传递模型，即sst-kω模型。该模型涉及两个运输方程，即湍动能和耗散率方程。其中采用了非稳态分离算法(usa)。关于模型具体详见文章：menterfr.two-equationeddy-viscosityturbulencemodelsforengineeringapplications[j].aiaajournal.1994,32.耦合速度表达式时用simple算法处理，并且针对不稳定情况采用了二阶隐式方案。二阶隐式方案用于k-ω传递方程和动量方程中的对流项。在步骤s1中，所述串列双圆柱的尾流驰振中，下游圆柱横向振动采用尾流圆柱横向振动方程求解，该下游圆柱横向振动方程为：其中，m为下游圆柱质量；y、和分别为垂直于流向的位移、速度和加速度；c和k分为为结构阻尼和刚度常数；fy为圆柱受到的垂直方向的力。在本实施例中，步骤s1中，所述串列双圆柱的尾流驰振几何模型中，参见图1可以看出，采用的矩形计算域为60d×40d，d为圆柱直径；最大阻塞率为5％，矩形计算域中，上游圆柱中心距离入口边界为25d，下游圆柱中心距离出口边界为31d；矩形计算域中，入口边界设为速度入口边界条件；矩形计算域中，流域出口边界设为压力出口边界条件；矩形计算域中，上下边界定义为对称边界条件；矩形计算域中，上游圆柱和下游圆柱均采用圆柱且圆柱表面采用无滑移壁面。参见图2可以看出，本实施例中，所述对几何模型进行网格划分时采用非均匀四边形结构化网格，上、下游圆柱表面的表面网格采用o型网格；所述对几何模型进行网格划分时加密圆柱近壁面的网格以及控制壁面第一层网格高度，使上、下游圆柱壁面第一层网格高度y ＜1；实现降低网格对数值计算的影响。矩形计算域中，为了模拟圆柱非稳态振动，入口下游27d和34d处的交界面采用了滑移网格技术，在横向边界处使用动态层铺法更新网格。对于圆柱绕流这种具有周期性现象的模拟，时间步取单位周期内至少有100个时间步。为了验证本文的数值模拟方法的可行性与网格划分策略的精确性，先模拟了两个固定圆柱的气动特性，模拟的圆柱间距选为l/d＝6,t/d＝0-4。雷诺数re＝3.48×104，re＝ρud/μ，ρ为流体密度，μ为流体粘性系数,u为来流的速度，处于亚临界范围，湍流度为1％。通过保证y ＜1，壁面第一层网格高度为0.013mm。通过加密圆柱周围网格，选用了4种不同数量的网格进行模拟，将得到平均阻力系数与现有实验结果进行对比，详见文献wuw,huangs,barltropn.currentinducedinstabilityoftwocircularcylinders[j].appliedoceanresearch.2002,24(5):287-297.网格数量从1.3×105增加到3.7×105，结合图3中(a)可看出，平均阻力系数不断接近现有实验结果。图3中(b)可以看出，阻力系数时程曲线表现为收敛。当网格数量由3.7×105增加至4.5×105时，两种网格对应的平均阻力系数非常接近，阻力系数时程曲线在稳定段几乎重合，此时网格对结果无太大影响。故后续动态圆柱的模拟选用与3.7×105网格数相同的网格划分策略，即圆柱壁面第一层网格高度y ＜1，壁面附近的o-block增长率为1.05，其他区域网格增长率为1.08。网格长宽比控制在5以内，并且相邻网格尺寸接近。s2：根据串列双圆柱的尾流驰振几何模型，设立至少一种工况，分析串列双圆柱的尾流驰振气动力的频谱规律；如表1，包括三种工况，基于这三种工况，进行了三种下游圆柱的sc数较小工况的尾流驰振模拟。表1三种模拟工况尾流圆柱的自身特性参数和sc数的数据表质量比m*阻尼比ζ固有频率fssc工况一84.350.24％6.9hz3.39工况二113.240.28％6.78hz4工况三142.120.26％7.12hz4.64表1中，质量比m*＝m/ρd2，m为下游圆柱单位长度质量。sc＝2m*(2πζ)，ξ为下游圆柱自身结构阻尼比。折算风速ur＝u/fsd。三种模拟工况的无量纲振幅a/d随折算风速变化的散点图详见图4。工况一的ur＝9.179和ur＝14.49对应全时程的位移时程曲线以及位移时程频谱详见图5和图6所示。其中，从图5可以看出，ur＝9.179时，对应的位移时程频谱表现的非常粗糙，运动表现为类似准周期运动。从图6可以看出，而当ur＝14.49时范围时，对应的位移时程频谱表现的较为光滑，运动表现为类似周期运动。两种位移现象的分界点为折算风速ur＝13.5，而该风速在模拟结果中恰巧是振幅增加速率与准定常理论一致的风速，说明在ur＝9.179时是受非定常效应影响较为严重的风速范围，而在ur＝14.49非定常效应在逐渐减弱。各风速对应的位移时程频谱主频都为自振频率，但还分别包含了很微弱的其他频率。尾流驰振发展阶段ur＝9.179和ur＝14.49对应全时程的升力频谱详见图7可知：基于所有的风速的气动升力频谱成分主要包括7个主要频率；其中，1-5阶自振频率，随着折算风速的增大，自振频率保持不变，将1-5阶自振频率对应的气动力定义为自激力；6-7阶频率随着折算风速呈线性增加，将6-7阶频率对应的气动力定义为涡激力。s3：根据步骤s1中的上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，建立含有未知参数的串列双圆柱的尾流驰振气动力数学模型；其中气动力表示为自激力和涡激力两种形式；在步骤s3中，含有未知参数的串列双圆柱的尾流驰振气动力数学模型为：未知参数包括：未知自激力参数和未知涡激力参数；fy为结构的总气动力；为自激力的表达式，fv(fv)为涡激力的表达式；分别对y和进行高阶泰勒序列展开，可得到尾流驰振力中的自激力多项式表达式如下：式中gij为气动参数。g00这一项与描述响应状态的y和y都没有关系，可以将该项认定为静力项，它对尾流驰振动力响应没有影响，故忽略该项。由之前的气动力频谱可以看出，对整个频率成分贡献的较大的基本在5阶倍频之内，故可以将泰勒展开式中超过y5和的项截断。基于能量等效原理，由三部分组成：气动阻尼力项(改变了结构的阻尼)、气动刚度力项(改变了结构的刚度)和纯力项(既不改变结构的阻尼也不改变圆柱的刚度)。则，自激力可以写成下式忽略纯力项，无量纲自激力表达式如下对于横向振动，气动刚度项仅轻微改变结构频率，可以不考虑，则上述公式可精简为：ρ为流体密度，u为来流的速度，d为圆柱直径，y为垂直于流向的位移，为垂直于流向的速度；α01为一阶无量纲气动阻尼参数，α03为三阶无量纲气动阻尼参数，α05为五阶无量纲气动阻尼参数，均为未知参数；由于最大的有两个涡脱频率fv1和fv2，由于涡激力始终是与涡脱频率息息相关的，故可先将其表示为涡脱频率的函数，则涡激力的表达式可写为：两个圆柱的斯托哈儿数st随折算风速变化详见图8，可以发现，上下游圆柱的斯托哈尔数并不是一个定值。上游圆柱是一个固定的静态圆柱，它的斯托哈尔数值由0.209以接近线性缓慢减小到0.184，随着风速不断增加也有趋于稳定的趋势；而下游圆柱为一个动态圆柱，它的斯托哈尔数值是由0.302缓慢增加最后平稳在0.315，与振幅随风速的变化趋势一致。无论是上游静态圆柱还是下游动态圆柱，它们的斯托哈尔数虽然随折算风速的变化缓慢，但是都不为一个定值，为了使公式的涡激力计算具有更高的预测精度，必须精确考虑到斯托哈尔数的变化。又由于涡激力的表达式中，相位角与都是需要拟合求解的参数，无量纲后可以改写为下式：式中：δ11,δ12,δ21,δ22为无量纲后的涡激力参数；st1为上游圆柱斯托哈尔数；st2为下游圆柱斯托哈尔数；ρ为流体密度；ur为折算风速，fs为固有频率；v11v12v21v22分别表示无量纲前的涡激力参数；则可以得到：下游圆柱在振动过程中，对其提供能量的输入或消耗的有气动阻尼力、结构自身阻尼和涡激力。s4：采用基于能量等效的最小二乘参数识别方法，拟合出自激力中气动阻尼力参数和涡激力参数；，得到最终的串列双圆柱的尾流驰振气动力数学模型；其中，基于能量等效的最小二乘参数识别方法详情参见文献gaogz,zhuld.nonlinearmathematicalmodelofunsteadygallopingforceonarectangular2:1cylinder[j].journaloffluidsandstructures.2017,70:47-71.s5：对步骤s4得到的串列双圆柱的尾流驰振气动力数学模型进行对比。步骤s5的具体内容为：s51：根据得到的串列双圆柱的尾流驰振气动力数学模型，采用runge-kutta法反算位移和速度；s52：定义下游圆柱能量的参数名称：ωs为下游圆柱的自然圆频率；ws(t)为下游圆柱自身结构阻尼在当前时间的累积功，消耗能量为负值；wa,v(t)为利用拟合的参数反算的下游圆柱的气动阻尼力和涡激力在当前时间的累积功，提供能量为正值；ee(t)为下游圆柱在当前时间的机械能，其中包括动能与势能；s53：将步骤s51得到的位移、速度带入公式(1)-(3)中，重构得到气动阻尼力和涡激力的累积功；s54：将气动阻尼力和涡激力的累积功与所述串列双圆柱的尾流驰振几何模型的模拟值进行对比。以ur＝13.889为例，参见图9，为各项做功与机械能的对比关系，根据拟合得到的参数重构的位移时程与fluent计算的位移时程的对比细节详见图10。实际计算的升力累积做功时程曲线与拟合参数重构的吻合较好，两者最大相对误差仅为3％，ee(t)的时程曲线也与wa,v(t)-ws(t)非常接近，两者最大相对误差仅为6％，这些功曲线“抖动”的误差来源于驰振力模型忽略的高频振荡成分，这些成分的贡献几乎可以忽略。对于图10的位移时程曲线，虽然两条曲线之间存在着明显的相位差，这是有忽略了气动刚度项，但是参数重构的驰振位移曲线幅值和发展过程方面都与fluent计算的吻合度较好，两者的最大相对误差仅为1％。但是，参见图11，在ur＝12.08时，基于能量等效的最小二乘法拟合的参数重构的涡激力只为气动阻尼力的1.1倍，这与前述升力fft变换出来的频谱现象(涡脱频率占比非常高)相悖。由于涡激力在整个驰振振动过程中的累积做功几乎为0，下游圆柱运动能量的贡献主要是来源于气动阻尼力项。基于能量拟合时，涡激力对能量不敏感，使得拟合的涡激力项的参数误差较大。综上，气动阻尼力项对下游圆柱的能量输入起绝对控制作用，涡激力项参数对能量变化不敏感，对模型中涡激力项参数的拟合需要从力本身出发，即对升力中扣除气动阻尼力项后的力直接利用最小二乘法进行拟合，再进行涡激力参数的规律分析。通过能量拟合得到气动阻尼参数，通过从力本身拟合得到其他的气动参数。为了更清楚的说明涡激力项和气动刚度项对位移的作用，将利用式(4)参数重构的位移与仅不包含涡激力项参数以及仅不包含气动刚度力项参数重构的位移进行对比分析。如图12所示，在尾流驰振发展时间段14s-14.6s，式(4)与仅不包含气动刚度力项模型参数重构的位移幅值最大相对误差为1％，但是两者之间存在较大的相位差，而式(4)参数重构的位移幅值要高于仅不包含涡激力项参数重构的15％，但是两者之间存在的相位差较小。在尾流驰振稳定时间段41s-41.6s，三者的相位差的关系与尾流驰振发展时间段一致，不同的是三者位移幅值都相接近，最大相对误差都为1％。这说明在尾流驰振的全阶段，气动刚度力项不会改变位移幅值，仅起到了改变相位差的作用，这与前述等效线性系统计算结果一致；涡激力项的涡脱强迫振动仅在驰振初期提供较大扰动作用，使得驰振得以快速发展，但是对最终的稳态位移振幅不产生影响，这也证实了涡激力在振动全过程做功接近于0的假设。参数重构力的时程曲线如图13所示，涡激力的峰值为0.151n，而自激力(气动阻尼力 气动刚度力 纯力)峰值只有0.032n，其值只有涡激力的20％，与前述涡脱频率占比非常高现象一致，直接对力的本身进行拟合得到的涡激力参数可以更好的反映出下游圆柱实际的振动状态。并且由图14可知，仅气动阻尼力参数和涡激力参数重构的力与fluent模拟结果的升力无论是在相位上还是幅值上都吻合非常好，误差最大为4％左右，说明自激力中气动刚度力成分与纯力成分的占比非常小几乎可以忽略。综上所述，无论是稳态位移还是气动力的预测，涡激力项和气动阻尼力项组成式(5)的简化模型已经具有足够精度。参见图15上游圆柱的涡脱引起下游圆柱表面的无量纲涡激力参数振幅av1随着折算风速的增加接近线性递减，而下游圆柱自身涡脱引起表面的无量纲涡激力参数振幅av2始终保持着值为0.195-0.205的小幅波动状态，说明在尾流驰振区间范围内，上游圆柱的涡脱效应随折算风速的增加对下游圆柱产生的表面压力影响开始减弱，下游圆柱自身涡脱的激励效应随折算风速变化则显示出较好的定常性，这些规律的凸显对构建涡激力参数模型有着极大的促进作用。对模型的验证内容如下：工况一的气动参数存在较明显的规律性，通过参数重构的力和位移都与模拟值接近。尾流驰振气动力模型如果要被应用于工程预测，不应该只对工况一具有适用性，气动参数的值在固定间距下应该只与结构外形和折算风速相关，其应该对工况二和工况三都具有相似的气动规律，故对工况二和工况三的圆柱尾流驰振同样进行无量纲气动参数的识别，观察它们的分布规律，以便获得具有普适应的数学预测公式。涡激力无量纲参数与气动阻尼力无量纲参数识别结果如图16。由图16可知，三种工况的无量纲气动阻尼参数和无量纲涡激力参数具有一定的离散性，这是由于在建立气动力模型时对尾流驰振力采取了一定的近似，忽略了更高阶的非线性频率和一些随机涡脱频率成分对气动力的贡献。但是图中参数大体上落在一个区域，呈现出较为明确且相似的客观规律性。三种工况的气动参数呈现的规律较为相似，说明基于工况一建立的气动力模型不仅仅只对工况一有适用性，证实了本文的尾流驰振气动力模型具有一定的普遍性和可应用性。α01、α03、α05、av1和av2根据散点规律都可以选择对应的多项式进行拟合，拟合过程可以对三种工况的气动参数先求平均值，再利用这些平均值进行最小二乘拟合。如图16所示，拟合的曲线较为光滑，可以较大程度的体现无量纲参数随折算风速的气动规律。一阶无量纲气动阻尼参数α01，三阶无量纲气动阻尼参数α03和五阶无量纲气动阻尼参数α05在ur＝13.5时都有一个跳跃性，该风速可以理解为非定常效应影响的分割风速。因此，接下来将对这两个阶段的无量纲气动阻尼参数分别进行分析。ur＜13.5：α01和α05随着折算风速的增加而递减，而α03是随着折算风速的增加而递增。α01在风速范围内恒为正值，说明气动阻尼力中始终为负阻尼，而三阶无量纲气动阻尼参数α03和五阶无量纲气动阻尼参数α05在风速范围内有正有负，是从负阻尼变为正阻尼，则是从正阻尼变为负阻尼，但是α03和α05几乎在任意折算风速始终保持着一正一负的相反状态。在一个特定折算风速下，下游圆柱在起振阶段时，非线性总阻尼为负，和或者和和对下游圆柱源源不断的在输送驰振振动振幅发展需要的能量，使下游圆柱的尾流驰振振幅不断增加。随着振幅不断增加，或者以及结构自身提供的正阻尼逐渐增加直到与负阻尼相等，下游圆柱非线性总阻尼为0，此时下游圆柱从空气吸收的能量与耗散的能量相同，下游圆柱处于动态平衡状态并以稳定振幅作极限环振荡，α03和α05的相互制约是保证尾流驰振在每个风速拥有振幅自限性的关键因素。从整个驰振发展风速范围来看，这种涡振与驰振耦合的非定常效应不仅体现在位移时程曲线上，也明显的体现在无量纲气动参数随折算风速的强烈变化上。ur＞13.5：α01、α03和α05的值远小于ur＜13.5的无量纲气动阻尼参数值，从整体上看，α01、α03和α05随风速增加最终趋于一个定值，而α03稳定值为一个接近于0的正阻尼值，而α05稳定值为一个接近于0的负阻尼值，两者仍然保持着反号。虽然所有无量纲气动阻尼参数在小范围内仍然在随折算风速变化，但是此时的变化幅度已经非常小了，标志着非定常效应的减弱，随着折算风速的继续增加，所有参数趋近于一个常数，故振幅曲线在后续呈现线性增加是可以解释的。将工况一、工况二和工况三的无量纲气动参数散点整合成统一可使用的数学公式后，便可带入非定常力学方程直接对工况一、工况二和工况三进行非定常尾流驰振振幅的预测。如图17可以看出，工况一、工况二和工况三的预测振幅与fluent模拟计算的振幅最大相对误差主要集中小范围的小振幅区域，因为此时振幅对微小的差值变化十分敏感。工况一在整个预测区间的最大振幅差值为0.06d，在ur＞9的最大相对误差只有7％，工况二在整个预测区间的最大振幅差值为0.06d，ur＞9.45的最大相对误差只有10％。虽然工况三在整个预测区间的相对误差波动较大，但是最大振幅差值也只有0.07d，说明平均法拟合虽然不能完全消除圆柱的自身特性对气动力模型无量纲气动参数的影响，但是三种工况绝大部分折算风速范围的振幅预测结果都较为理想且满足工程精度要求。而且预测结果能较大程度还原振幅随折算风速变化的客观规律，如低折算风速下振幅增加的斜率变化、由非定常振动过渡到定常振动的平台以及高折算风速下振幅的近线性增加等。综上所述，尾流驰振气动力模型对工况二和工况三也具有较高的适用性，同时也进一步验证了气动力模型预测结果的可靠性与精确性。综上，经过上述方法，可以得到：(1)建立一种新的针对低sc数尾流驰振非定常气动力非线性模型。气动力表示为自激力和涡激力和的形式。基于能量等效原理将自激力分为气动阻尼力、气动刚度力和纯力项，其中气动阻尼项对尾流驰振稳态振幅起主要作用。将涡激力表达为结构涡脱频率的函数。结合固定上下游圆柱的气动力频谱，涡激力主要包含两个圆柱的涡脱频率。(2)识别了气动力模型的气动参数，验证了识别方法的可靠性以及气动参数的准确性。涡激力只在振动初期提供扰动，在整个振动过程贡献的能量几乎为0。对于气动阻尼参数的拟合选用基于能量等效的最小二乘法，对于涡激力参数的拟合选用从力的时程本身进行最小二乘法拟合。利用拟合的参数无论是各项力做功重构、位移重构还是力重构，其结果都表明，气动阻尼项的参数识别和涡激项的参数识别都具有较高的识别精度，验证了对各项气动参数所采用拟合方法的可靠性与准确性。(3)总结了工况一、工况二和工况三的无量纲气动参数共有规律。工况一、工况二和工况三识别的气动参数大体上落在一个区域，呈现出较为明确且相似的客观规律性，这些规律可用具体数学表达式表示成为统一的实用性公式。通过公式预测三种工况的振幅结果都较为理想且满足工程精度要求，而且能较大程度还原振幅随折算风速变化的客观规律，如低折算风速下振幅增加的斜率变化、由非定常振动过渡到定常振动的平台以及高折算风速下振幅的近线性增加等。尾流驰振非定常气动力非线性模型对工况二和工况三也具有较高的适用性，进一步验证了模型预测结果的可靠性与精确性。最后需要说明的是，上述描述仅仅为本发明的优选实施例，本领域的普通技术人员在本发明的启示下，在不违背本发明宗旨及权利要求的前提下，可以做出多种类似的表示，这样的变换均落入本发明的保护范围之内。当前第1页1 2 3
技术特征：1.一种串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于，具体步骤为：
s1：采用计算流体动力学软件，模拟圆柱流向间距l/d＝4的情况下，上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，得到串列双圆柱的尾流驰振几何模型并对该几何模型进行网格划分；
其中，l为两圆柱在流向的垂直距离，d为圆柱直径；
s2：根据串列双圆柱的尾流驰振几何模型，设立至少一种工况，分析串列双圆柱的尾流驰振气动力的频谱规律；
s3：根据步骤s1中的上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，建立含有未知参数的串列双圆柱的尾流驰振气动力数学模型；其中气动力表示为自激力和涡激力两种形式；
s4：采用基于能量等效的最小二乘参数识别方法，拟合出自激力中气动阻尼力参数和涡激力参数，得到最终的串列双圆柱的尾流驰振气动力数学模型；
s5：对步骤s4得到的串列双圆柱的尾流驰振气动力数学模型进行对比。
2.根据权利要求1所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：在步骤s1中，所述串列双圆柱的尾流驰振中，下游圆柱横向振动采用尾流圆柱横向振动方程求解，该下游圆柱横向振动方程为：
其中，m为下游圆柱质量；y、和分别为垂直于流向的位移、速度和加速度；c和k分为为结构阻尼和刚度常数；fy为圆柱受到的垂直方向的力。
3.根据权利要求1所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：步骤s1中，所述串列双圆柱的尾流驰振几何模型中，采用的矩形计算域为60d×40d，d为圆柱直径；
矩形计算域中，上游圆柱中心距离入口边界为25d，下游圆柱中心距离出口边界为31d；
矩形计算域中，入口边界设为速度入口边界条件；
矩形计算域中，流域出口边界设为压力出口边界条件；
矩形计算域中，上下边界定义为对称边界条件；
矩形计算域中，上游圆柱和下游圆柱均采用圆柱且圆柱表面采用无滑移壁面。
4.根据权利要求1所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：步骤s1中，所述对几何模型进行网格划分时采用非均匀四边形结构化网格，上、下游圆柱表面的表面网格采用o型网格；
所述对几何模型进行网格划分时加密圆柱近壁面的网格以及控制壁面第一层网格高度，使上、下游圆柱壁面第一层网格高度y ＜1；
矩形计算域中，入口下游27d和34d处的交界面采用了滑移网格技术，在横向边界处使用动态层铺法更新网格。
5.根据权利要求1所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：在步骤s2中，分析得到的串列双圆柱的尾流驰振气动力的频谱规律为：
基于所有的风速的气动升力频谱成分主要包括7个主要频率；
其中，1-5阶自振频率，随着折算风速的增大，自振频率保持不变，将1-5阶自振频率对应的气动力定义为自激力；
6-7阶频率随着折算风速呈线性增加，将6-7阶频率对应的气动力定义为涡激力。
6.根据权利要求5所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：在步骤s3中，
含有未知参数的串列双圆柱的尾流驰振气动力数学模型为：
未知参数包括：未知自激力参数和未知涡激力参数；
fy为结构的总气动力；为自激力的表达式，fv(fv)为涡激力的表达式；
其中，
ρ为流体密度，u为来流的速度，d为圆柱直径，y为垂直于流向的位移，为垂直于流向的速度；α01为一阶无量纲气动阻尼参数，α03为三阶无量纲气动阻尼参数，α05为五阶无量纲气动阻尼参数，均为未知参数；
式中：δ11,δ12,δ21,δ22为无量纲后的涡激力参数；
st1为上游圆柱斯托哈尔数；st2为下游圆柱斯托哈尔数；ρ为流体密度；ur为折算风速，fs为固有频率；v11v12v21v22分别表示无量纲前的涡激力参数；
即：
7.根据权利要求6所述的串列双圆柱尾流驰振非定常气动力数学模型建立方法，其特征在于：步骤s5的具体内容为：
s51：根据得到的串列双圆柱的尾流驰振气动力数学模型，采用runge-kutta法反算位移和速度；
s52：定义下游圆柱能量的参数名称：
ωs为下游圆柱的自然圆频率；
ws(t)为下游圆柱自身结构阻尼在当前时间的累积功，消耗能量为负值；
wa,v(t)为利用拟合的参数反算的下游圆柱的气动阻尼力和涡激力在当前时间的累积功，提供能量为正值，
ee(t)为下游圆柱在当前时间的机械能，其中包括动能与势能；
s53：将步骤s51得到的位移、速度带入公式(1)-(3)中，重构得到气动阻尼力和涡激力的累积功；
s54：将气动阻尼力和涡激力的累积功与所述串列双圆柱的尾流驰振几何模型的模拟值进行对比。
技术总结本发明提出一种串列双圆柱尾流驰振非定常气动力数学模型建立方法，具体步骤为：采用计算流体动力学软件，模拟上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，得到串列双圆柱的尾流驰振几何模型并对该几何模型进行网格划分；根据几何模型，设立工况，分析串列双圆柱的尾流驰振气动力的频谱规律；根据上游圆柱固定、下游圆柱横向振动的串列双圆柱的尾流驰振，建立含有未知参数的串列双圆柱的尾流驰振气动力数学模型；采用基于能量等效的最小二乘参数识别方法，拟合出自激力中气动阻尼力参数和涡激力参数，求解串列双圆柱的尾流驰振气动力数学模型。建立串列双圆柱尾流驰振模型，利于现有工程建设，提高工程建造稳定性和可靠性。
技术研发人员：晏致涛;杨小刚;钟永力;孙毅;刘欣鹏;李妍;赵爽;王灵芝;聂小春;巩黎明;曾玉洁;罗钧
受保护的技术使用者：重庆科技学院
技术研发日：2021.05.06
技术公布日：2021.08.03
