基于虚功原理的流固耦合面力和位移传递方法

仲继泽, 谢志强, 沈渡, 徐自力

(1.中国舰船研究设计中心船舶振动噪声重点实验室, 430064, 武汉;2.西安交通大学机械结构强度与振动国家重点实验室, 710049, 西安)

在工程实际中,流固耦合问题有时会引起结构破坏,造成巨大经济损失。1940年,美国Tacoma Narrows bridge在19 m/s的风中受到涡街流动导致的周期性气动力的作用而产生涡激振动,最终坍塌[1]。由此可见,流固耦合问题具有十分重要的研究价值。

对流固耦合问题的研究[2-4]通常需要通过力和位移传递方法在流场与结构网格之间传递力和位移信息[5]。目前,文献报道的流固耦合面的信息传递方法主要有曲面拟合、表面跟踪和局部插值3种[6]。曲面拟合方法对复杂曲面会产生病态的插值矩阵使得数值计算无法进行[7],所以曲面拟合方法的应用并不广泛。表面跟踪方法通过映射将流场网格节点坐标投影到临近的结构单元上,这种方法对高曲率曲面的计算精度很低[8]。局部插值方法对高曲率曲面网格插值时,会引起曲面出现锯齿以及网格负体积的问题[9],从而产生较大的计算误差。上述这些方法本质上都是插值算法,会引入额外的插值误差。

本文提出了一种用于流固耦合面的力和位移传递方法,该方法在结构网格发生变形的同时会带动流场网格的运动,从而避开了流场与结构网格之间位移反复插值的问题。对于结构与流场网格这一整体系统而言,可以直接采用流固耦合面处流场网格节点的气动力向量作为激励,而不需要将流体压力插值到结构网格上,所以本文方法不会引入额外的插值误差。通过理论分析发现,本文方法能够自动满足力和位移传递的能量守恒要求。采用本文方法,通过流固耦合计算了弹性梁流固耦合振动响应,预测了Agard Wing 445.6颤振边界,并与表面跟踪方法以及局部插值方法的结果进行了比较。

1　能量守恒的力和位移传递方法

1.1　位移传递

对于流固耦合面处的任意一个流场网格节点A,都存在一个与A距离最近的结构网格单元,结构变形前后流场网格节点位置的变化如图1所示。由图1可知,在结构未发生变形时,结构网格单元存在附近一点a,其位置与A完全重合,点a从属于该结构网格单元,并随着结构的变形产生位移。当结构与流场网格相互独立时,结构发生变形之后,由于点a随结构变形产生位移而流场网格节点A保持不动,所以点a、A的位置不再重合。如果在结构变形过程中,人为设定流场网格节点A、a的位移相同,那么结构变形之后,点A、a的位置仍然保持重合,这就相当于将流场与结构网格粘贴在一起,使之成为一个整体系统。

图1　结构变形前后流场网格节点位置的变化

对于结构与流场网格整体系统而言,结构为弹性体结构,参照弹性体动网格方法[10],将流场网格充满虚拟弹性物质,该物质允许流体在其内部自由流动,最后通过有限元方法构建结构与流场网格整体的振动控制方程。这一整体系统的振动方程规模巨大,直接求解会占用大量计算资源,为解决这一问题,可以采用本课题组的快速动网格方法[11-12]计算结构与流场网格的节点位移。因此,流场网格节点位移可以绕过插值过程直接得到,从而达到了将结构网格节点位移向流场网格节点传递的目的。

1.2　气流力传递

本文将结构与流场网格作为一个整体系统,流场网格被粘贴在结构网格上,所以无须将流固耦合面处流场网格面的流体压力插值到结构网格节点上去,相反可以直接采用流固耦合面处流场网格面的流体压力作为结构与流场网格整体系统的激励。在实际的计算过程中,可以通过力等效的方法计算出流固耦合面处流场网格节点的气动力向量,计算过程为

(1)

式中:sin为流固耦合面;Fnode为流固耦合面流场等效气动节点力向量;p1为分布在流固耦合面的微元面ds上的流场压力;ds为流固耦合面处微元面的面积向量。

流固耦合面处流场网格动网格节点M如图2所示,节点M是流场网格单元面1、2、3、4的公共节点,图2中由虚线所围成的四边形是以节点M为形心的微元面,其面积向量用dsM表示。根据式(1),节点M的气动节点力向量FM可以采用下式进行计算

(2)

式中:p2为微元面dsM上的流体压力;m为以节点M为顶点的流场网格单元面的个数;ni为流场网格单元面i的顶点的个数;ps为流场网格单元面i上的流体压力;si为流场网格单元面i的面积向量。流固耦合面处所有流场网格节点的节点力向量的集合构成了整体系统的气动力向量Fin,可以用FM表示为

(3)

图2　流固耦合面处流场网格动网格节点M

1.3　能量守恒性验证

为保证流固耦合计算的准确性,通过流固耦合面的力和位移传递过程必须满足能量守恒,即等效计算得到的气动力向量对结构的虚功必须等于结构表面流体压力对结构的虚功,即

(4)

式中:δxnode为流固耦合面上节点的虚位移向量;δxds为流固耦合面上微元面ds的虚位移向量。

对于流固耦合面上任意一个流场网格节点M,气动节点力向量在该点处的虚功可表示为

(5)

对于微元面dsM,流体压力在该微元面上的虚功为

(6)

微元面dsM的形心位于节点M处,有δxds=δxM,根据式(2)、(5)、(6),可得

(7)

式(7)表明,在任意一个流场网格节点M处,力和位移的传递方法都能满足能量守恒的要求,所以力和位移的传递方法可以保证结构与流体之间能量传递的守恒性。

结构在流固耦合面处受到的惯性力作用主要是虚拟弹性体运动产生的惯性力和流体运动产生的惯性力,其中虚拟弹性体的密度为0,其惯性力也为0,因此虚拟弹性体的惯性力可以忽略。流体运动产生的惯性力以流体压力形式作用在结构上,本文方法对流体压力信息进行了传递,考虑了惯性力的作用,适用于动力学问题的研究。

2　流固耦合计算方法

本文方法将结构与流场网格作为一个整体系统,通过求解该系统的振动方程,计算结构和流场网格节点的位移,采用本课题组的快速动网格方法[11]以减少动网格计算时间。为进一步加快流固耦合计算速度,在研究中采用了本课题组的时空同步流固耦合算法[12]。同时,为了加快流场迭代的收敛速度,采用稳态流场的计算结果作为流固耦合中瞬态流场的初值。本文采用Fluent内置的Coupled Method计算流场,每一时间步的时空同步流固耦合计算流程如图3所示。

图3　时空同步流固耦合计算流程

3　算例验证

3.1　弹性梁流固耦合振动计算

采用一个标准的弹性梁流固耦合算例来验证本文方法。该标准算例是由Turek等[13]提出的,Turek等给出了弹性梁以及流体域的几何尺寸,如图4所示,同时提供了弹性梁的密度为1.0×104kg/m3、弹性模量为1.4×106Pa、泊松比为0.3、流体密度为1.0×103kg/m3、动力黏度为1.0 Pa·s等参数。

图4　弹性梁及流体域的几何参数

弹性梁流场边界设置是:流体域左端面为流速进口边界,流体沿着垂直于左端面的方向流入,流速沿y方向变化,呈抛物线分布;流体域右端为压力出口边界;流体域上、下两个端面为无滑移壁面边界,弹性梁外表面和圆形壁面也是无滑移壁面边界;弹性梁外表面与流体域连接,设置为流固耦合面。弹性梁发生颤振时,流场出口压力为0,流场进口流速分布满足

(8)

本文采用的弹性梁流场网格如图5所示,共有4节点矩形单元19 627,节点总数为20 034。

图5　弹性梁流场网格

首先进行稳态流场计算,为下一步的瞬态流场计算提供初始值。稳态流场计算得到流场流速分布如图6所示。在流体域左端面的进口处,流速呈对称的抛物线分布,最大流速点位于流体域的中轴线处,弹性梁的中轴线距流体域上端面0.21 m,距流体域下端面0.2 m。进口处的最大流速点位于弹性梁中轴线的上方,所以弹性梁上部流场流速的分布与下部流场流速的空间分布是相似的,上部流场流速较高。由图6可知,在圆形壁面的右端,弹性梁的上下两侧,存在一个低流速区域,对应的流线分布图中显示该低流速区域存在两个旋涡。

图6　弹性梁稳态流场流速分布

图7　弹性梁颤振时流场流速分布

流固耦合计算得到弹性梁颤振时附近流场流速的分布,结果如图7所示。弹性梁附近的流场中共有4个旋涡,在弹性梁下部接近中间的位置存在1个最大的旋涡,在弹性梁的上部,圆形壁面右侧存在1个较小的旋涡,另外的2个旋涡已经从弹性梁最右端脱落下来。由此可见,弹性梁附近流场旋涡的周期性形成、发展及在弹性梁右端脱落,使得作用在弹性梁流固耦合面的流场压力发生周期性波动,在周期性波动的压力作用下使得弹性梁发生了颤振。

本文计算得到的弹性梁颤振的振动幅值和周期如表1所示,本文方法的计算结果与文献[13]结果相对偏差约为0.5%。本文采用局部插值方法和表面跟踪方法,计算了弹性梁的流固耦合振动,得到了弹性梁颤振时的振动幅值和频率,与文献结果的偏差分别为1.2%和1.4%,比本文方法的计算偏差高出约0.8%。本文中弹性梁算例的计算采用同一套流场网格、同样的时间步长和相同的离散格式,因此本文方法和表面跟踪法、局部插值法的计算结果的差别是由数值方法不同带来的。弹性梁的流固耦合面比较平缓,局部插值方法和表面跟踪方法的高偏差是由插值误差引起的。由于本文方法没有额外的插值误差,所以计算准确性优于局部插值和表面跟踪方法。

表1　弹性梁振动幅值及周期

3.2　Agard Wing 445.6颤振边界预测

机翼前缘和尾缘一般都是高曲率曲面,因此采用国际标准颤振模型Agard Wing 445.6[14]来验证本文方法对高曲率流固耦合面的适用性。AGARD Wing 445.6的展弦比为1.65,根梢比为0.66,该机翼的1/4弦线后掠角为45°,具体几何尺寸如图8所示。Agard Wing 445.6的横截面为对称的NACA 65A004翼型,采用木质材料制造,展向弹性模量为3.15 GPa,厚度方向的弹性模量为0.416 2 GPa,各个方向的剪切模量为0.439 2 GPa,泊松比为0.31,密度为393.5 kg/m3,材料阻尼比为0.02。

图8　Agard Wing 445.6的几何尺寸

本文Agard Wing 445.6的计算网格数量为563 472,机翼表面以及横截面的流场网格分布情况如图9所示。首先进行稳态流场分析,然后采用稳态流场的结果作为瞬态流场的初值,对亚声速(Ma=0.678)、跨声速(Ma=1.072)和超声速(1.141)工况下的Agard Wing 445.6颤振现象进行流固耦合分析。

图9　Agard Wing 445.6的流场网格分布

(1)亚声速工况。当马赫数为0.678时,计算得到的Agard Wing 445.6翼尖截面的流速分布如图10所示。由图10可知:在翼尖前缘的顶点,流体流速为0,此处为流场滞止点;随着流体沿机翼两侧向尾缘流动,过了翼型中间位置之后,边界层的厚度逐渐增加,最终在翼尖尾缘形成一个流速较低的尾迹流动区域。从机翼颤振时的流场流速分布来看,流动在机翼前缘发生分离,分离后的流体沿机翼弦向向尾缘流动,并在机翼的上侧形成平均流速较低的分离区。由此可知,在Ma=0.678的亚声速工况下,机翼附近流场的流动分离是诱发机翼颤振的主要原因。

(a)稳态流场

(b)颤振时流场图10　Ma=0.678时Agard Wing 445.6截面的流速分布

(2)跨声速工况。当马赫数为1.072时,计算得到的Agard Wing 445.6翼尖截面的流速分布如图11所示。由图11可知:随着流体沿机翼两侧向尾缘流动,流体流速迅速增加到最大,然后又逐渐减小,并在距机翼前缘65%弦长的位置产生了激波;过了流场激波所在位置之后,边界层的厚度逐渐增加,最终在翼尖尾缘形成了一个流速较低的尾迹流动区域。从机翼颤振时流场的流速分布来看,机翼上侧流场的激波位于距翼型前缘的60%弦长处,流动从激波的位置产生分离,并向翼型尾缘流动,形成分离区;机翼下侧流场的激波位于距翼型前缘的65%弦长处,没有出现明显的流动分离现象。在机翼的振动过程中,机翼附近流场的激波在距翼型前缘60%～65%弦长的区间内周期性振荡,并诱发了流动分离现象。由此可知,在Ma=1.072的跨声速工况下,机翼附近流场的振荡激波及其诱发的流动分离是诱发机翼颤振的主要原因。

(a)稳态流场

(b)颤振时流场图11　Ma=1.072时Agard Wing 445.6截面的流速分布

(3)超声速工况。当马赫数为1.141时,计算得到的Agard Wing 445.6翼尖截面的流速分布如图12所示。由图12可知:在距机翼前缘80%弦长的位置产生了激波,过了流场激波所在位置之后,边界层的厚度逐渐增加,最终在翼尖尾缘形成了一个流速较低的尾迹流动区域。从机翼颤振时流场的流速分布来看,机翼上侧流场的激波位于距翼型前缘72%弦长处,流动从激波的位置产生分离,并向翼型尾缘流动,形成分离区;机翼下侧流场的激波位于距翼型前缘84%弦长处,没有出现明显的流动分离现象。在机翼振动的过程中,机翼附近流场的激波在距翼型前缘72%～84%弦长的区间内周期性振荡,并诱发了流动分离现象。由此可知,在Ma=1.141的超声速工况下,机翼附近流场的振荡激波及其诱发的流动分离是诱发机翼颤振的主要原因。

(a)稳态流场

(b)颤振时流场图12　Ma=1.141时Agard Wing 445.6截面的流速分布

本文计算得到的Agard Wing 445.6的无因次颤振流速如表2所示,本文方法的亚声速、跨声速颤振流速的计算结果与实验值[14]相对偏差约为1.3%,超声速颤振流速的计算结果与实验值的相对偏差为5.18%。本文采用局部插值方法和表面跟踪方法计算得到的Agard Wing 445.6的无因次颤振流可知,亚声速、跨声速颤振流速的计算结果与实验值的相对偏差分别为3.0%和3.4%,比本文方法的计算偏差高出约1.9%;超声速颤振流速的计算结果与实验值的相对偏差分别为8.46%和9.40%,比本文方法的计算偏差高出约3.75%。本文Agard Wing 445.6颤振算例的计算采用同一套流场网格、同样的时间步长和相同的离散格式。因此,这些数值方法不同带来本文方法和表面跟踪法、局部插值法的计算结果的差别。

表2　Agard Wing 445.6无因次颤振流速

机翼前缘和尾缘的曲面为高曲率曲面,由于本文方法没有额外的插值误差,而且具有对高曲率曲面的适用性,所以计算准确性优于局部插值方法和表面跟踪方法。

4　结　论

本文提出了一种自动满足能量守恒的流固耦合面力和位移传递方法,克服了目前主流的力和位移传递方法会引入额外插值误差的不足,而且本文方法对高曲率的流固耦合面也具有较好的适用性。本文对标准弹性梁流固耦合振动算例和Agard Wing 445.6颤振问题进行了流固耦合分析,并对比了本文方法和主流的表面跟踪法、局部插值法的计算准确性。研究发现,对于弹性梁附近流场旋涡的周期性形成、发展和脱落是弹性梁发生颤振的主要原因。本文方法计算得到的弹性梁颤振周期和幅值与文献结果相对偏差约为0.5%,而由于存在额外插值误差,导致局部插值方法和表面跟踪方法的计算结果与文献的相对偏差比本文方法的计算偏差高出约0.8%。

对Agard Wing 445.6的研究表明,亚声速时机翼附近的流动分离诱发了颤振,跨声速、超声速时机翼附近的振荡激波及其诱发的流动分离导致了颤振。本文方法计算得到的亚声速、跨声速颤振流速与实验值相对偏差约为1.3%,超声速颤振流速的计算结果与实验值相对偏差为5.18%。由于机翼前缘和尾缘曲面的高曲率和额外插值误差的综合影响下,局部插值方法和表面跟踪方法的亚声速、跨声速颤振流速计算结果的相对偏差比本文方法高出约1.9%,超声速颤振流速计算结果的相对偏差比本文方法高出约3.75%。

参考文献:

[1]BILLAH K Y, SCANLAN R H. Resonance, Tacoma narrows bridge failure, and undergraduate physics textbooks [J]. American Journal of Physics, 1991, 59(2): 118-124.

[2]王俊毅, 招启军, 肖宇. 基于CFD/CSD耦合方法的新型桨尖旋翼气动弹性载荷计算 [J]. 航空学报, 2014, 35(9): 2426-2437.

WANG Junyi, ZHAO Qijun, XIAO Yu. Calculations on aeroelastic loads of rotor with advanced blade-tip based on CFD/CSD coupling method [J]. Acta Aeronautica and Astronautica Sinica, 2014, 35(9): 2426-2437.

[3]仲继泽, 徐自力. 流固单向耦合的能量法及机翼颤振预测 [J]. 西安交通大学学报, 2017, 51(1): 109-114.

ZHONG Jize, XU Zili. Wing flutter prediction using an energy method based on one-way fluid structure coupling [J]. Journal of Xi’an Jiaotong University, 2017, 51(1): 109-114.

[4]刘占生, 马瑞贤, 杨帆, 等. 基于流固耦合作用的柔性体流噪声降噪机理研究 [J]. 机械工程学报, 2016, 52(10): 176-184.

LIU Zhansheng, MA Ruixian, YANG Fan, et al. Study on noise reduction mechanism of flow induced noise for flexible body based on fluid-structure interaction [J]. Journal of Mechanical Engineering, 2016, 52(10): 176-184.

[5]周岱, 李磊, 邓麟勇, 等. 流固耦合问题的网格更新与信息传递新方法 [J]. 工程力学, 2010, 27(5): 83-90.

ZHOU Dai, LI Lei, DENG Linyong, et al. Novel methods for mesh update and data transfer technique of fluid-structure interaction [J]. Engineering Mechanics, 2010, 27(5): 83-90.

[6]DE BOER A, VAN ZUIJLENA H, BIJL H. Review of coupling methods for non-matching meshes [J]. Computer Methods in Applied Mechanics and Engineering, 2007, 196(8): 1515-1525.

[7]安伟刚, 梁生云, 陈殿宇. 一种局部动态数据交换方法在流固耦合分析中的应用 [J]. 航空学报, 2013, 34(3): 541-546.

AN Weigang, LIANG Shengyun, CHEN Dianyu. Local dynamic data exchange in fluid structure interaction analysis [J]. Acta Aeronautica and Astronautica Sinica, 2013, 34(3): 541-546.

[8]KIM Y H, KIM J E. New hybrid interpolation method for motion transfer in fluid-structure interactions [J]. Journal of Aircraft, 2006, 43(2): 567-569.

[9]杨敏, 王福军, 戚兰英, 等. 流固耦合界面模型及其在水力机械动力学分析中的应用 [J]. 水利学报, 2011, 42(7): 819-825.

YANG Min, WANG Fujun, QI Lanying, et al. Fluid-structure coupling interface model and its application in dynamic analysis of hydraulic machinery [J]. Journal of Hydraulic Engineering, 2011, 42(7): 819-825.

[10] TEZDUYAR T E. Stabilized finite element formulations for incompressible flow computations [J]. Advances in Applied Mechanics, 1991, 28: 1-44.

[11] 仲继泽, 徐自力, 陶磊. 基于虚拟弹性体的快速动网格方法 [J]. 西安交通大学学报, 2016, 50(10): 132-138.

ZHONG Jize, XU Zili, TAO Lei. An efficient dynamic mesh method based on pseudo elastic solid [J]. Journal of Xi’an Jiaotong University, 2016, 50(10): 132-138.

[12] 仲继泽, 徐自力. 采用快速动网格技术的时空同步流固耦合算法 [J]. 振动工程学报, 2017, 30(1): 41-48.

ZHONG Jize, XU Zili. Time-space synchronizing fluid structure coupling method using a fast dynamic mesh technique [J]. Journal of Vibration Engineering, 2017, 30(1): 41-48.

[13] TUREK S, HRON J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow [M]. Berlin, Germany: Springer, 2006: 371-385.

[14] YATES J E C. AGARD standard aeroelastic configurations for dynamic response: IWing 445.6 [R]. Washington, DC, USA: NASA, 1987: 1-74.

[本刊相关文献链接]

范增华,荣伟彬,王乐锋,等.压电微喷辅助液滴的多物理场耦合与实验.2016,50(11):56-61.[doi:10.7652/xjtuxb2016 11009]

仲继泽,徐自力,陶磊.基于虚拟弹性体的快速动网格方法.2016,50(10):132-138.[doi:10.7652/xjtuxb201610020]

姜涛,黄伟,王安麟.多路阀阀芯节流槽拓扑结构组合的神经网络模型.2016,50(6):36-41.[doi:10.7652/xjtuxb201606 006]

高炎,晏鑫,李军.燃气透平叶片尾缘开缝结构冷却性能的数值研究.2016,50(3):29-37.[doi:10.7652/xjtuxb201603005]

郭涛,管志成,孙光普,等.调频振子-液体联合水平减振的流固耦合机理研究.2016,50(1):28-33.[doi:10.7652/xjtuxb 201601005]

宋明毅,吴伟烽,李直.汽车空调压缩机气阀运动规律模拟.2015,49(12):144-150.[doi:10.7652/xjtuxb201512023]

季家东,葛培琪,毕文波.换热器内弹性管束流体组合诱导振动响应的数值分析.2015,49(9):24-29.[doi:10.7652/xjtuxb201509005]