耦合响应范文

耦合响应范文（精选5篇）

耦合响应 第1篇

事实上, 由于旋翼系统本身具有高柔性, 因此, 这些刚性条件下的气动力模型都不能完全符合真实情况;而且由于结构弹性和气动力耦合影响, 使得弹性模型动态失速情况更加复杂、难解。目前关于弹性模型动态失速的研究寥寥无几, 主要集中在另行建立弹性结构模型, 然后与半经验气动力模型、CFD模型进行耦合计算方面[5,6]。但在目前的研究中, 无论是半经验气动力模型/CSD耦合计算, 还是CFD/CSD耦合方法, 都是在时间域内不断在结构场和气动场之间迭代的松耦合方法。而实际上, 气动和弹性结构必然是同步互相影响的;而且由于动态失速现象的复杂性, 即使是细微的同步耦合影响也可能对结果造成很大差别。因此, 需要建立动态失速下弹性模型的气动和结构同步求解的紧耦合模型, 才能更准确地进行动态失速下气动弹性响应计算分析, 为弹性模型动态失速机理分析打下基础。

本文在经典气动弹性响应模型和动态失速模型的基础上, 新建二维弹性翼型动态失速下气弹响应的紧耦合模型;并进行了相应的紧耦合计算, 分析了刚性动态失速和弹性动态失速模型的区别, 为以后的弹性体动态失速问题研究提供新的紧耦合理论思路, 并为以后的三维弹性体动态失速紧耦合分析打下基础。

1 数学模型

1.1 经典气动弹性响应模型及数学表达

完整的二维气动弹性响应模型如图1所示。

图1中, V为来流速度, L为升力, M为力矩, K为扭转刚度, Kh为弯曲刚度, θ为刚性的攻角变化, 为气动弹性影响引起的攻角变化量。A.C.为气动中心, E.A.为结构刚心, C.G.为结构重心。

二维情况下, 气动弹性耦合影响主要为弹性影响下的攻角变化与升力之间的互相影响, 而只有K影响攻角增量, Kh则不产生影响。因此, 在以气动弹性耦合影响为研究目标时, 可简化模型, 不考虑Kh, 如图2所示。

根据图2得到经典的二维气动弹性响应数学表达式[7]:

式 (1) 中, IФ为翼段对刚心 (E.A.) 的惯量, IФ=I0+mσ2, σ为刚心到重心 (C.G.) 的距离, m为翼段质量, Ф为气动弹性耦合引起的攻角变量, e为气动中心 (A.C.) 到刚心的距离, L为升力。

1.2 修正的气动弹性响应模型及数学表达

根据二维旋翼动态失速情况建立响应计算的物理模型, 如图3所示。图2中, θ为刚性攻角变化。因此, 气动弹性耦合影响下的实际攻角为α=θ+Ф。公式 (1) 中, 第一项为惯性项, 第二项为弹性项, 第三项为合力矩。动态失速情况下, θ本身为变量, 通常实验中表达为:。θ变化也将引起作用在重心上的力的变化, 即惯性项改变。这点在根据图1 (图2) 的模型建立的公式 (1) 中无法显示, 因此, 需要对模型进行修正, 如图3所示。

图3中, θ作为刚性底座运动形式, 由此, 得到修正的气动弹性响应表达式:

公式 (2) 体现了动态失速情况下θ和Ф分别对气动弹性耦合的贡献。

由此, 可得完整的动态失速下的气弹表达式:

式 (3) 中,

一般情况下, 气动弹性响应应该为小量。这是因为也与结构扭转的情况有关;如果太大, 即结构扭转量很大, 意味着结构受到很大影响, 可能会缩短结构寿命甚至结构破坏, 这在结构设计中需要避免。因此, 对动态失速下的气动弹性响应问题, 作小扰动假设:Ф相对于θ为小量。

令Ф=Δθ, 则Cl (α) 作泰勒级数展开:

式 (4) 中, 下标b表示当前刚性底座的运动形态。

忽略公式 (4) 中高阶小量, 则公式 (3) 可化为

根据公式 (5) 可知, 若已知某状态下的刚性动态失速特性, 则可预测此时弹性体的附加变化量。整理公式 (5) , 得:

公式 (6) 为变系数常微分方程, 通解难以求得。因此, 可用Runge-Kutta方法求数值解。令:

则得:

式 (8) 中, a、b、c为定值, 分别为

2 计算模型

2.1 计算模型

计算模型选择翼型NACA0012, 基础模型运动θ=9°+5°sinωt, Ma$=0.297, c=0.61 m, k=0.051, ω=16.9 Hz, 见参考文献[8]中风洞实验7019。以上条件与NASA实验一致。假设模型长度为1 m, 则, 模型材料为铝材, 则模型质量m=98.8 kg, KФ=1.49105Nm/rad, IФ=1.110-5kgm2, 模型固有频率ωs=40.5 Hz。

2.2 气动力模型

应用LS法 (最小二乘法) 拟合参考文献中给出的实验Cl曲线的函数形式。将Cl拟合成为周期ωt的八阶傅里叶函数, 根据拟合函数求得的气动力与实验值基本一致, 如图4所示。

3 仿真分析

根据公式 (8) 计算二维弹性旋翼动态失速下的响应, 计算方法为二阶Runge-Kutta法。计算结果如图5~图7所示。

图5展示了由模型弹性引起的攻角增量随时间的变化。由于本文中弹性模型的固有频率与模型运动频率, 即弹性模型的激励频率, 相差较远。因此攻角增量的变化频率应与激励频率基本一致, 如图5所示。此时, 攻角增量的变化形态也会接近正弦波, 符合正弦激励下的气弹响应理论。此外, 攻角增量的变化范围在2度以内, 符合小扰动要求。

图6为弹性模型和刚性模型攻角随时间的变化对比。由图6可知, 考虑弹性后, 攻角变化范围比刚性结果要更大, 即, 实际弹性模型的动态失速现象可能会比纯刚性模型更严重。

图7为弹性模型和刚性模型升力随时间的变化对比。弹性影响导致攻角变化增大, 因此, 令动态失速下弹性模型的升力范围也产生了变化。

4 结论

本文根据已知的动态失速特性, 修正了经典的气动弹性响应模型, 初步建立了二维弹性翼型动态失速下气弹响应的紧耦合模型。应用此紧耦合模型对一经典动态失速风洞模型进行计算分析, 结果表明, 气动弹性耦合影响下的攻角增量频率、形态符合理论推导, 紧耦合响应计算有较高可靠性。并且同样的风速条件下, 弹性模型比刚性模型的动态失速现象会更明显, 实际弹性模型的攻角及升力变化范围可能增大, 弹性影响不可忽视。

摘要：根据经典气动弹性响应模型和刚性动态失速模型, 建立了二维弹性翼型动态失速下气弹响应的紧耦合模型;计算模型为典型动态失速风洞模型, 应用最小二乘法拟合实验气动力得到了气动力模型;应用此气动力模型进行气动弹性响应紧耦合计算, 结果表明, 同样的条件下, 弹性模型比刚性模型的动态失速现象会更明显, 实际飞行攻角及升力的变化范围可能增大, 弹性影响不可忽视。

关键词：动态失速,气动弹性,弹性响应,紧耦合

参考文献

[1] Friedmann P P, Hodges D H.Rotary wing aeroelasticity—a historical perspective.Journal of Aircraft, 2003;40 (6) :1019—1046

[2] McAlister K W, Lambert O, Petot D.Application of the ONERA Model of Dynamic Stall.NASA Technical Paper 2399, 1984

[3] Sheng W, Galbraith R A M, Coton F N.A modified dynamic stall model for low mach numbers.AIAA Paper 2007—0626, 2007

[4] Mulleners K, Kindler K, Raffel M.Dynamic stall on a fully equipped helicopter model, Aerospace Science and Technology, 2012;19 (1) :72 —76

[5] Tang D M, Dowell E H.Experimental investigation of three-dimensional dynamic stall model oscillating in pitch.Journal of Aircraft, 1995;32 (5) :1062—1071

[6] Ahmad J, Biedron R T.Code-to-code comparison of CFD/CSD simulation for a helicopter rotor in forward Flight.AIAA Paper, 2011:2011—3819

[7] 陈桂彬, 杨超, 邹丛青.气动弹性设计基础.北京:北京航空航天大学出版社, 2004:82—89

耦合响应 第2篇

摘要：基于提出的抖振力模型和建立的风车桥耦合振动模型，发展了一种可以考虑抖振力空间相关性的风车桥耦合振动分析方法，并编制了相应的计算程序.以江顺长江大桥为工程背景，测试了桥梁抖振力的空间相关性和考虑车桥耦合作用的车桥气动参数，分析研究了桥梁抖振力空间相关性对侧风作用下桥梁和车辆耦合动力响应的影响.研究结果表明：桥梁抖振力空间相关性对桥梁动力响应有显著影响，对车辆的动力响应也有一定的影响.

关键词：大跨度斜拉桥；风工程；风车桥系统；耦合动力分析；抖振响应；相关性

中图分类号： U441.3文献标识码：A

风荷载作用下汽车桥梁系统的振动、车辆行驶安全性和舒适性以及桥梁的疲劳损伤研究是桥梁工程的前沿课题之一.过去10年内，国内外学者在这方面进行了大量的研究并取得了大量的研究成果.夏禾， Xu和Guo， Cai和Chen、李永乐以及韩万水等均建立了风车桥系统空间耦合振动分析框架，编制了相应程序，但目前风车桥系统空间耦合振动分析中均假设抖振力的空间相关性与来流脉动风速的空间相关性相同.这与实际情况不符，很多风桥领域的研究已经证明这一点.

Larose认为，假设抖振力的空间相关性等于来流风的空间相关性是导致预测抖振响应不确定性的一个主要原因，特别是对于类似闭口箱型桥面这类结构更是如此.Kimura等人对扁平六角形和矩形断面进行了风洞试验，结果表明2种断面的抖振力中升力沿桥跨方向相关性要明显大于来流脉动风速的相关性.此外，许志豪，韩艳，赵传亮，张冠华和赵林等，李少鹏和李明水等[12]都对抖振力空间相关性进行了一定的研究，研究结果均表明：抖振力的空间相关性明显大于来流脉动风速的空间相关性.

本文提出了一套考虑桥梁抖振力空间相关性的风车桥耦合振动分析方法，并编制了相应的计算程序.以江顺长江大桥为工程背景，测试了桥梁抖振力的空间相关性，分析研究了桥梁抖振力空间相关性对侧风作用下桥梁和车辆耦合动力响应的影响.研究结果表明：桥梁抖振力空间相关性对桥梁动力响应有显著影响，对车辆的动力响应也有一定的影响.因此，在今后的风车桥系统耦合振动研究中有必要考虑桥梁抖振力的空间相关性.

2 算例分析

2.1计算参数

本文以江顺大桥为工程背景，江顺大桥是一座主跨700 m的钢混凝土混合斜拉桥，其跨径布置为：60 m （混凝土） +176 m + 700 m +176 m + 60 m （混凝土），主桥全长1 172 m，总体布置图如图1所示.表1给出了江顺大桥考虑车辆影响的桥梁断面三分力系数（取自文献[14]），图2给出了考虑车辆影响的江顺大桥气动导数（取自文献[15]）.桥梁抖振力的空间相关性可以用相关函数表示，相关函数通常表示成指数衰减函数

coh（f，Δ）=exp （-Cf·ΔU）. （9）

式中：C为衰减因子，即抖振力的空间相关系数，如式（6a～f）；f为抖振力的频率；Δ为两排测压孔之间的距离；U为平均风速.在主梁节段模型的不同断面上布置了测压孔，通过对各个截面进行压力积分得到各截面的抖振力，然后对抖振力进行相关分析并采用非线性最小二乘拟合得到抖振力的空间相关系数，见表1.文中对脉动风场的相关性没有进行相关测试，根据经验取顺风向和竖风向脉动风空间相关系数均为16.

采用ANSYS建立桥梁有限元模型，计算得到的桥梁结构前10阶动力特性见表2，结构阻尼比取0.005.首先分别取前10和20阶桥梁振型进行车桥耦合振动计算，计算得到的桥梁和车辆动力响应非常接近，故为了节省计算时间，桥梁振型只取前10阶进行相关计算.

本文采用侧面面积比较大的集装箱2轴车辆作为分析对象，车辆的模型及主要参数参考文献[13]，车辆的气动力系数参考文献[14]，如图3所示.

（a） 立面图（b） 标准钢箱梁断面图

图1江顺大桥的整体布置图（单位：cm）

Fig.1The sketch of Jiangshun Bridge （unit： cm）

顺风向的功率谱采用Kaimal谱，竖风向采用Lumley和Panofsky提出的功率谱[16]，采用Davenport相关函数，顺风向和竖风向脉动风速的空间相关系数均取16.频率采样点数为1 024，截止频率为2π，平均风速为20 m/s.采用谐波合成法[17]沿桥纵向共模拟了99个点的脉动风时程.假定桥面粗糙度是一零均值的平稳高斯随机过程，利用离散傅立叶逆变换计算得到路面不平度值，本文只考虑良好路面一种情况，粗糙度系数取20× 10-6 m3/circle.

2.2计算结果及分析

2.2.1抖振力空间相关性对桥梁动力响应的影响

计算中桥梁气动参数取表1和图2给出的气动参数，车辆气动参数取图3中拟合的车辆气动力系数，平均风速取5 m/s， 10 m/s， 15 m/s， 20 m/s， 25 m/s和30 m/s，车辆行驶速度为20 m/s，考虑一辆车辆，车辆距桥中心线的距离为15 m.受篇幅限制，本文只给出部分计算结果，图4给出了平均风速为20 m/s时桥梁位

图4抖振力空间相关性对桥梁位移响应的影响

Fig.4Effect of the coherence of buffeting forces on the displacement responses of the bridge

图5桥梁位移响应随风速变化情况

Fig.5The displacement responses of the bridge vs. the wind velocity

移响应沿桥纵向位移的变化情况.图5给出了桥梁跨中位移响应随风速变化的情况.

耦合响应 第3篇

1 高压涡轮航空机械的边界条件概述

对于某型航空发动机高压涡轮盘, 它由盘体和封严盘两部分组成。根据其结构特点及换热机理, 将轮盘简化为自由盘和旋转盘腔两种模型, 即将涡轮盘从前后安装边以下取为自由盘, 而以上部分取为旋转盘腔, 同时在轮盘外缘处考虑接触热阻, 盘两侧考虑冷却空气沿盘面流程的温升。

2 应力场的计算

应力场的计算, 主要采用的方法为航空机械耦合法, 这种方法的应用可以有效的计算出在航空机械载荷耦合作用下, 应力场的具体系数。而在计算的过程中, 可以根据热力条件来对温度载荷进行分析和处理, 在对温度载荷进行处理的过程中, 需要先对冷空气的温度以及燃气温度等具体数值进行确定, 在确定好相关的热力系数之后, 将其作为计算的条件, 通过该条件来实现对耦合结构荷载的计算和求解, 最后, 再利用相关的耦合结构单元对航空发动机涡轮的应力场进行具体的计算。在计算时, 所采用的计算方法可以是直接耦合场分析法, 利用这种方法可以准确的计算出应力场。在计算阶段中, 需要对温度和位移自由度的耦合场单元进行构造上的分析, 并建立相应的模型, 依据热力条件, 实现叠加式的温度计算和位移求解, 这样可以使得航空机械的载荷耦合响应分析结果更加的符合实际。

2.1 涡轮航空机械有限元模型的建立

在对航空发动机涡轮航空机械进行机械载荷耦合响应分析时, 可以采用有限元分析法, 建立有限元模型, 利用该模型, 将叶片榫头与涡轮轮盘的榫头相连接, 在两者结合的过程中, 也需要将严盘以及前后档板加入到连接中, 在模型建立的时候, 可以忽略涡轮上的通气孔, 这样可以做到建立一个简易的模型, 并且模型的样式为对称式。

2.2 边界条件和载荷的确定

在进行边界条件以及载荷的确定时, 需要从位移约束条件、热力边界等方面入手, 下面就对这些方面进行具体分析。

2.2.1 位移约束条件的确定。

位移约束条件的确定, 需要依据发动机的具体工作情况来制定约束条件, 所建立的涡轮航空机械模型为对称式的模型, 而且在模型中, 涡轮的前盘边界与前轴是相互连接的关系, 前轴设置为高压压气机, 将高压压气机和高压涡轮中的转子进行结构的连接, 使得前轴作为承力件进行使用。依据前轴设置的条件, 需要在涡轮的前盘进行Y方向的约束。

2.2.2 热力边界的确定。

针对温度载荷的作用影响, 要想使其能够得到有效的处理, 就需要以热力条件为基础, 在对冷却空气温度进行确定后, 将其作为耦合单元进行热力边界条件的计算依据, 为了保障计算的准确性, 还需要对轮缘处燃气温度以及换热系数来对人力边界进行处理。通常情况下, 在确定热力边界时, 需要将冷却气流引入到压气机中, 将压气机口的温度作为已知数值, 将冷却气流的温度设定在552.3k上, 而将轮缘处燃气温度值设定在1086K上, 而要想能够对热力边界进行确定, 就需要在对热力边界条件进行处理的过程中, 忽略冷气在引气管中所产生的温度变化, 同时要将轮缘温度控制在常温上, 另外, 还需要对冷气进行对流换热处理, 分段对轮盘的进行压力的施加, 只有做好以上处理, 才能够有效的保障热力边界确定的有效性。

2.2.3 荷载的确定。

在对荷载进行确定的时候, 需要依据航空发动机具体的运行运行情况来进行涡轮载荷的计算, 事先要确定发动机的转速值, 在已知转速的基础上, 对轮盘外缘的温度载荷进行分析和确定。

2.3 计算结果和分析

通过对机械载荷耦合作用下响应进行分析, 可以了解到, 温度的响应相较于转速的响应来说, 具有一定的滞后性。因此, 对于机械载荷耦合作用下的响应分析可以采用有限元分析法来进行, 而有限元分析法中也具有多种不同的方法, 其中最适用于载荷耦合响应分析的为有限元瞬态分析法。这种分析方法的历时时间为1145s, 而分析结果则主要体现在最后一个状态上。

就有限元分析法可以了解到, 航空发电机涡轮的最大应力主要集中在涡轮盘的涡轮的中心位置, 在该位置应力值最大可达到1050MPa, 而就其运转中所产生的温度的平均值来说, 则平均温度会在300℃左右。航空发动机涡轮所采用的材料主要为镍基高温合金材料, 这种材料具有较高的抗高温性, 同时还具有较高的强度, 能够承受涡轮运转时所产生的高温以及应力, 涡轮运转过程中, 所产生的最大值在材料可承受的范围之内, 最大应力值可达到涡轮可承受的极限, 在极限范围内, 涡轮的中心处, 是处于一种塑性的状态, 因此, 在对涡轮中心位置的应力场进行分析时, 可以采用弹塑性应力场分析法, 但是这种分析方法相较于其他的分析方法来说, 具有保守性的特点。

3 有限元分析方法的应用

本文对航空发动机涡轮盘的热力边界条件进行了分析和简化, 计算了涡轮盘的换热系数, 并用所获得的换热系数, 进行了该涡轮盘的机械载荷耦合响应的有限元分析。该方法也适用于其它涡轮盘的热力边界条件分析;这里采用有限元瞬态分析方法, 采用第三类热力边界条件来处理温度载荷, 利用机械耦合算法, 机械载荷, 计算分析了涡轮盘在机械载荷耦合作用下的稳态响应, 同时该方法也可得到瞬态下的温度场和应力场结果。

结束语

综上所述, 本文主要针对航空发动机涡轮的航空机械载荷耦合响应进行了分析, 通过建立模型, 展开了对发动机涡轮运行系数的计算, 依据计算所得的系数, 分析出了航空机械载荷耦合响应的具体情况, 本文中主要采用的分析方法为有限元分析法, 而这种方法的应用范围较广, 能够有效的计算出具体的航空机械载荷耦合响应系数, 从而对航空发动机涡轮的运行形成合理的维护。

摘要：针对航空发动机涡轮航空机械载荷的耦合作用进行分析的过程中, 还具有一定的难度, 目前针对航空发动机涡轮航空机械载荷耦合响应分析的研究论文相对较少, 这就使得航空机械载荷耦合的作用无法得到有效的发挥, 从而严重影响到航空发动机涡轮的运行质量。本文就主要针对航空发动机涡轮航空机械载荷耦合响应进行了简要的分析和研究, 仅供同行交流和借鉴。

关键词：航空机械,发动机涡轮,载荷耦合相应

参考文献

[1]栾永先.先进涡轮盘结构强度对比分析[J].航空发动机, 2013 (3) .

[2]胡殿印, 裴月, 王荣桥, 李其汉.涡轮盘结构概率设计体系的研究[J].航空学报, 2008 (5) .

[3]王筱冬.涡轮盘应力场及破裂转速的有限元分析[J].榆林学院学报, 2009 (4) .

耦合响应 第4篇

储液结构广泛应用于人们日常生活的各个方面,作为一种重要的构筑物,在日常的生产工作中发挥着重要的作用[1]。储液结构涉及到石油、化工、水利、天然气以及核电等各种不同的行业,主要用于储存石油、天然气以及核反应燃料等易燃、易爆和具有放射性的物质,具有极大危险性,一旦泄露将会造成火灾或者极大的环境污染等二次灾害,如:在1964年日本新潟发生地震[2],储液罐以及电力设施遭到破坏,储油罐区起火,造成极大经济损失;2011年日本福岛核事故中[3],储存乏燃料组件的乏燃料水池受损导致降温系统丧失功能,造成重大核泄漏事故。因此储液结构的安全受到极大的关注和重视。我国处于世界两大地震构造系的交汇部位[4],属于地震多发国家。从而储液结构的抗震问题更是受到国内学者的重点关注。陈贵清和杨雪梅[5]对超大型储液罐考虑流固耦合进行了振动分析,用ANSYS建立储液罐—液体有限元模型,对储液罐流固耦合系统的振动频率以及破坏现象做了分析;刘云贺、王克成和陈厚群对储液池的抗震问题进行了探讨[6],探讨了液面波动、水池刚度以及宽深比变化对于储液池地震响应的影响;贾善坡,赵友清和许成祥[7]对储液容器内的液体晃动问题进行了动力学建模与模拟,以ABAQUS为平台,编写分析程序对液体晃动的特征频率和模态进行了计算分析;王翠翠,雷新弋[2]对立式的锚固储液罐基于有限元软件ADINA进行了地震响应分析。

对于储液结构的动力响应分析,涉及到流固耦合的强非线性问题,主要的分析方法是实验方法和数值分析方法。本文用大型有限元软件ADINA对某矩形储液池结构在不同情况下进行动力响应分析,期望能为以后的类似问题提供参考。

1 模型概况及ADINA建模

1.1 模型概况

本文中储液水池初始外廓长×宽×高为23 m×13.5 m×23.5 m,侧壁厚度为2 m,底板厚度为1.5 m。水池材料为C35混凝土,池内液体为水。逐渐改变储液高度和水池内部储液尺寸以研究在不同情况下,水池和内部液体的动力响应变化。

1.2 ADINA建模

应用大型有限元软件ADINA建立矩形储液水池的几何模型。水池C35混凝土采用3D-Solid单元模拟,池内水体采用3D-Fluid单元模拟,为考虑水体液面的晃动效应,水体材料采用ADI-NA软件中的势流体,并将水体上表面设置为势流体的自由液面来考虑重力对于液面波动的影响。储液结构体系的ADINA三维有限元模型如图1所示。

1.3 输入地震荷载

本文地震输入选用典型的EL-Centro波南北向时程记录,输入该地震波的前50 s时程记录,设置300个时间步长。输入地震荷载加速度时程如图2所示。

2 ADINA有限元分析结果

2.1 液体晃动频率与理论值对比分析

在保证水池几何尺寸不变的前提下,对水深分别为17 m,16 m,14 m,12 m和10 m的储液体系进行模态分析,可得出各个水深下水体晃动的一阶频率。在文献[8]中同济大学的李遇春等人基于Housner的质量弹簧系统提出了渡槽中液体晃动的频率理论计算公式(1),本文依据该公式针对本文中的储液水池结构计算水池中液体晃动的理论解以便和ADINA有限元模型的结果进行对比分析。经计算得出各水深情况下液体晃动的一阶频率如表1所示。

由表1中数据可以看出,ADINA有限元模态分析结果与文献[8]给出的公式计算结果基本吻合。两种计算值之间存在一点误差,这是因为在文献[8]中推导该计算公式时有两个经验参数,该经验参数的取值将会影响到晃动频率的计算结果。由此可见,有限元分析软件ADINA用于流固耦合问题的模拟分析是可行的,且准确度很高。

2.2 储液深度对液体晃动频率的影响

保持储液水池几何尺寸为初始几何尺寸23 m×13.5 m×23.5 m,储液池内液体深度分别取17 m,16 m,14 m,12 m和10 m,由ADINA计算得到各个储液深度下液体晃动的一、二、三和四阶频率和模态形状。由计算结果可知,在所取的几个储液深度下,液体晃动一、二、三和四阶模态振型一致,在此仅给出水深17 m情况下的一、二、三和四阶模态如图3~图6所示。各个储液深度下,各阶频率值如表2所示,各阶频率随储液深度的变化曲线如图7所示。

由图7可以看出随着储液深度的变化,液体第一阶晃动频率有细微变化,但很不明显,液体第二、三和四阶晃动频率随着储液深度的变化基本呈一条直线,没有变化,说明储液深度对于液体晃动频率的影响很小,基本可以忽略不计。

2.3 储液池几何尺寸对液体晃动频率的影响

保持储液池池壁厚度、底板厚度以及储液深度不变的情况下,改变储液池内部几何尺寸,研究储液池内部几何尺寸的变化对于液体晃动频率的影响。分别取内部几何尺寸为:19 m×9.5 m,21 m×10.5 m,23 m×11.5 m,25 m×12.5 m,27 m×13.5 m和29 m×14.5 m,计算得到液体晃动的一、二、三和四阶频率如表3所示,晃动频率随几何尺寸变化如图8所示。

由图8中变化曲线可以看出,水体晃动频率受到水池内部几何尺寸的影响较大,频率变化很明显,随着几何尺寸的增大,液体各阶晃动频率逐渐减小,相应的周期随之变大。与储液深度对液体晃动频率的影响对比可知,水池几何尺寸对于液体晃动频率的影响较储液深度的影响大,在进行储液结构体系设计时应加以考虑。

2.4 地震荷载作用下不同储液深度液体晃动响应

保持储液水池尺寸为初始几何尺寸23 m×13.5 m×23.5 m,分别取储液深度为:6 m,8 m,10 m,12 m和14 m。在模型Y向输入地震荷载,进行储液水池流固耦合的时程分析计算。选取储液池其中一个角部位置液面一点,对不同储液深度情况下的液面该点晃动时程进行对比分析,如图9所示。由图9中对比曲线可以看出:在不同储液深度时,液面同一位置点的晃动规律一致;在储液深度为6 m时,液体晃动波高最大;随着储液深度的增大,晃动波高减小;在储液深度达到10 m后储液深度的变化对于晃动波高影响不大,说明在保证几何尺寸以及输入地震荷载一定的情况下,此时储液深度的改变对于液体晃动波高的影响较小。

3 结语

本文对矩形储液水池在不同储液深度和不同几何尺寸情况下进行了模态以及地震响应分析,得出如下几个结论:

1)大型有限元软件ADINA可用于流固耦合分析,并且分析精度较高;

2)在保证储液水池几何尺寸一定的情况下,储液深度对于液体晃动频率的影响不大;

3)储液水池的几何尺寸对于液体晃动频率的影响较储液深度的影响大,且随着几何尺寸的增大,晃动频率逐渐减小;

4)在地震荷载作用下,不同储液深度的液体晃动规律一致,储液深度从一定值(本文为6 m)开始,随着深度的增加,液体晃动最大波高增大,当达到一定值(本文为10 m)时,储液深度的变化不再对液体晃动波高产生较大影响。

摘要：采用大型有限元软件ADINA,对某矩形储液结构体系进行了分析,通过对比液体晃动频率计算值与理论值,验证了ADINA对于流固耦合问题分析的实用性,并对不同储液深度和几何尺寸情况下的储液结构体系进行了模态和地震响应分析,得出了一些有意义的结论。

关键词：矩形储液水池,流固耦合,地震响应,模态分析

参考文献

[1]胡明祎.贮液结构地震反应数值分析及应用方法研究[J].中国地震局工程力学研究所2012届博士论文摘要(Ⅱ),2012(4):47-48.

[2]王翠翠,雷昕弋.考虑液固耦合储液罐非线性地震反应分析[J].防灾科技学院学报,2011,13(1):19-22.

[3]陈海英,刘远远,张春明,等.福岛乏燃料水池事故探讨[J].核安全,2012(2):76-78.

[4]李国强,李杰.地震工程学导论[M].北京:地震出版社,1992:1-348.

[5]陈贵清,杨雪梅.超大型储液罐流固耦合振动分析[J].唐山学院学报,2011,24(3):5-8.

[6]刘云贺,王克成,陈厚群.储液池的抗震问题探讨[J].地震工程与工程振动,2005,25(1):149-154.

[7]贾善坡,赵友清,许成祥.储液容器内液体晃动问题的动力学建模与模拟[J].郑州大学学报(工学版),2013,34(3):76-80.

耦合响应 第5篇

地面热模拟试验作为研究结构热问题的有效方法, 目前已成为高超声速飞行器结构设计、可靠性分析及产品性能检验的重要手段之一［3］。但是由于飞行器的飞行环境十分复杂, 单纯的地面试验是不可能全真模拟这种环境条件, 只能在简化的基础上近似模拟。所以, 用计算机数值模拟就变得相当重要。关于气动加热的最常见的数值模拟方法是求解N-S方程 ( CFD) , CFD的突出优点是精度较高, 但计算量太大, 需要大规模的计算机硬、软件资源, 耗时费力; 尤其是求解气动加热问题, 在目前条件下, 还不适合于工程设计和工程设计单位使用, 特别是方案设计阶段的初步估算［4］。

传统的热流工程计算方法是设定一个初始壁温用以计算冷壁热流 ( 冷壁热流: 将壁温设置为环境温度, 以等壁温边界计算得到的热流分布［5］) , 然后再用冷壁热流计算结构的温度场、热应力和热模态, 忽略了气动加热使壁温和表面性态变化后温度梯度、黏度和摩擦系数、附面层厚度等流场参数的变化［6］, 给热流、温度场等计算结果都带来较大的误差。

为了配合热结构试验的研究, 促进结构热试验技术向更高水平发展, 对试验件加热情况进行理论分析与数值仿真至关重要。仿真技术与结构热试验技术的密切结合, 已经成为解决结构热问题的发展方向。

高超声速飞行器结构热分析技术首先要考虑气动加热, 气动热流对结构传热则体现出实时耦合性。同时, 气动加热还是一个持续过程, 持续加热牵涉到各传热环节之间的耦合作用［7］。持续受热的特征决定了飞行器在某一时刻的受热情况不仅与这一时刻的飞行条件有关, 还与这之前飞行器自身热环境所经历的时间历程紧密相关。

用自主研发了飞行器气动加热与结构热响应耦合分析程序 ( CASRAS) ［8］, 利用它可以求解飞行器表面的温度分布, 以及在气动热载荷作用下飞行器的结构响应。结构分析程序采用了国际上航空航天领域大型商用软件Nastran。在得到飞行器结构温度分布之后, 利用Nastran可以进行结构的热应力分析、热模态分析、热稳定性分析等。

1 理论描述

1. 1 热响应分析

所谓热响应分析就是计算结构的热应力、热模态、热屈曲。当各时刻结构的温度场确定后, 增加适当的约束条件, 应用NASTRAN程序进行热响应数值分析。这样, 问题归结为确定结构的温度场。

1. 2 导热微分方程

传热过程的基本变量就是温度, 它是物体中的几何位置以及时间的函数。

在分析工程问题时, 经常要了解结构内部的温度分布情况, 例如高速飞行器的温度分布、发动机的工作温度、流体温度分布等。物体内部的温度分布取决于物体内部的热量交换, 以及物体与外部介质之间的热量交换, 一般认为是与时间相关的。物体内部的热交换采用以下的热传导方程 ( Fourier方程) 来描述［9］

式 ( 1) 中 ρ 为密度 ( kg/m3) ; c为比热[J/ ( kg·K) ];T为温度 ( K) ; t为时间 ( s) ; λx, λy, λz为导热系数[W/ ( m·K) ];为内热源 ( W/m3) 。

对于各向同性材料, 不同方向上的导热系数相同, 同时把式 ( 1) 简化为一维问题并忽略结构内部热生成, 则有

式 ( 2) 是典型的抛物型方程, 抛物型方程中的因变量与时间有关, 属于步进问题, 数学上称为初值问题, 其求解区域是一个开区间。计算时从已知的初值出发, 逐步向前推进, 依次获得适合于给定边界条件的解。在这类问题中, 特征线是与步进方向垂直的, 其依赖区与影响区以特征线为分界线, 如图1 所示。

步进问题的依赖区与影响区以特征线为界面截然分开是具有明确的物理意义的: 对于瞬态导热问题, 某一瞬时物体中的温度分布取决于该瞬时以前的情况及边界条件, 而与该瞬时以后将要发生的情况无关。步进问题的这种特点, 使得可以能够对气动加热工程方法与结构温度场进行时域上的耦合计算, 即从初始温度出发, 在时域上层层推进, 一直计算到所需时刻为止。

1. 3 气动加热工程计算方法描述

气动加热工程计算方法包括平板热流密度参考焓法、平板湍流热密度法、球头热流密度法、驻点热流密度法以及来流流动参数与正激波后流动参数的关联方法［10］, 这些方法需要大气参数表 ( 包含: 飞行高度、大气静温、大气静压、大气密度、声速) , 干空气物理参数表 ( 包含: 大气温度、大气比热、普朗特数、大气传导率) 及轨道参数表 ( 包含: 飞行时间、飞行高度、飞行马赫数) , 同时, 针对计算结构的特点, 工程计算方法中还包含有特征长度、锥角、后掠角几何参数等等。对于比较复杂的结构, 在耦合计算时采用了多种方法的混合。公式3 为某导弹遵循的气动方程

式 ( 3) 中: qc为冷壁热流 ( W/m2) ; hw为壁焓 ( J/kg) ; hc为恢复焓 ( J/kg) ; Tw为壁面温度 ( K) ; σ为波尔兹曼常数[W/ ( m2·K4) ]; ε 为黑度系数。

在耦合计算中, qr即为热壁热流, 即施加给结构表面的热流。可以看出, 式 ( 3) 中的壁面温度Tw及热壁热流是相互制约的, 按照步进问题的求解方法, 二者可以进行时域上的耦合计算。

2 耦合计算实施

2. 1 总体思路

耦合性是指飞行器各个传热环节之间存在的一种相互影响、相互制约关系, 在数学上表现为各环节求解时的联解性。对每一种计算模型, 按照气动加热工程计算方法计算初始温度下的热壁热流, 把热壁热流作为载荷施加到结构表面, 之后把整个求解时间分成n段, 通过n个计算任务分段完成。每个计算任务自动输入前一时间步的热流密度, 按照热流载荷进行下一步计算, 计算出新的壁面温度, 根据该计算任务所对应的时间段的轨道参数, 利用气动加热工程算法, 计算得到新的热流密度。最后, 所有分段计算结果按照弹道飞行时间顺序进行排序并以Patran格式进行自动转换, 形成完整的计算结果 ( 因为Patran每次计算的初始时间皆从零开始) 。该方法无需对壁温进行近似处理, 但整个分析是进行了多任务处理, 即需要多次Patran计算。在任务分段足够细致情况下, 可实现每个求解步真正意义上的气动/热耦合迭代分析, 这样, 可以得到较高精度的耦合分析结果。

2. 2 程序实施步骤

a) 按照公式 ( 3) 计算初始热壁热流, 此时Tw取当前环境温度;

b) 在Patran平台建立有限元模型;

c) 以集中热流方式对模型施加热流载荷;

d) 形成BDF卡片数据, 进行首次计算;

e) 由控制程序自动提取当前计算结果作为下一步的初始温度;

f) 提取壁面温度, 计算下一时刻热流密度;

g) 由控制程序自动形成新的Nastran卡片数据, 进行Nastran计算;

h) 转到执行步骤e) 直至整个时间推算完毕;

i) 控制程序对计算结果进行整理。

从上述处理过程不难发现, 耦合性的处理就是根据各环节传热情况实时调整 ( 修正) 各环节求解的边界条件, 每一环节在时间推进求解时都是在时变的边界条件下完成的, 处理是围绕结构壁面的温度响应计算进行的, 用每一刻得到的结构物面温度计算下一时刻的热流密度, 之后再计算新的物面温度, 直至过程结束。

2. 3 程序界面

应用PCL ( patran command language) 编程, 使得所有操作都能够在Patran平台上完成。程序界面见图2。

2. 4 气动热流载荷实现

首先按照图3 的热流分区在初温下计算气动热流 ( 由于气动加热公式中包含特征点到驻点的距离参数, 所以对有限元模型按照平行驻点前缘线的方法对模型进行分区, 即计算一系列热流值。) , 之后在Patran平台点击Loads/BCs, 出现窗体, 按标号次序输入各分区的热流值。分区原则是飞行器驻点前缘为一个热流区、一般部位以前缘为参考做一系列平行线进行多个分区, 在每一分区内选取一个有限元节点作为特征点, 特征点的作用是作为判据提取该点的温度值, 同时计算每一分区到前缘的垂直距离作为特征长度, 以雷诺数为判据计算出特征点的热流密度 ( 图4) 。

2. 5 Nastran分段计算

控制程序把Nastran计算模式限定为DOS单任务方式, 待Nastran计算结束后, 读取Nastran的结果文件 ( Jobname. F06) , 把结构表面温度作为下一步热壁密度计算的输入项, 把结构温度作为下一时间步计算的初始温度, 自动形成新的Nastran卡片数据 ( Jobname. BDF) , 图5 给出控制程序主流程。

2. 6 温度结果读取

提取结构文件, 采用PCL编程, 按照Patran模板格式 ( Pthermal_nod_T. RES_TMPL) 形成数据。瞬态分析的结果按秒时间数形成多个文件, 通过循环调用载入Patran数据库, 在Patran平台的结果显示 ( results) 进行云图显示和曲线绘制。

3 工程应用

3. 1 有翼导弹温度场分析

3. 1. 1 有限元模型和相关参数

图6 是某导弹的基本结构尺寸, 翼为梯形平板翼。图中导弹头部为尖拱弧形, 弧形的圆心坐标为:X = 2. 86D, Y = - 7. 792D, Z = 0. 0, 半径为R =757. 851 23 mm。坐标系原点位于头部, 导弹轴向为X方向, 弹体直径D = 90 mm。

针对导弹的特点, 采用了三种气流热流工程计算方法, 包括球头驻点、飞行器一般部位和机翼前缘。图7 给出示意图。

3. 1. 2 计算结果及分析

图8 是气动加热时间为96 s时, 在同一时刻, 弹翼温度相对于弹头降低20 K。

3. 2 导弹设备舱环境温度计算

3.2.1有限元模型和相关参数

飞行器仪器舱内设备正常的工作环境温度一般不能超过100 ℃, 当舱体外表面受到气动加热时, 舱壁温度急剧升高, 将会使舱内温度越限, 造成元器件性能恶化甚至失效。

本算例为某导弹在其外部增加一个钢盒, 钢盒外部承受气动热载荷, 钢盒内部包含有电子设备, 图9 给出有限元模型图。蒙皮和隔热材料选用壳单元, 空气选用体单元, 电子设备舱体等效为含有内部热生成的体单元。在空气和电子设备之间建立了接触边界, 接触热阻系数通过典型试验获得。

3. 2. 2 计算结果及分析

计算结果表明: 电子设备舱最高温度为71 ℃, 结构表面最高温度为341 ℃。由于接触热阻的影响, 电子设备舱表面温度明显低于其相邻空气的温度。图10 给出了电子设备舱温度曲线与试验曲线比较, 计算结果与试验测试结果非常接近, 相对误差在7%以内。

3. 3 陶瓷基TPS隔热性能计算

3.3.1有限元模型和相关参数

陶瓷基复合材料盖板式热防护系统是一种具有隔热与承载一体化功能的新型热防护结构, 主要由两块C/SiC盖板搭接而成, 在搭接形成的空腔填充密封材料, 盖板与底面板之间采用高温合金支架与紧固件连接, 紧固件与支架间安装密封垫, 盖板与底面板中间填充隔热材料。图11 给出TPS结构示意图。

由于不同材料的热沉不同, 以体单元的形式建立精细有限元模型。载荷为TPS面板同时承受均布气动压力和气动热载荷。

3. 3. 2 计算结果及分析

计算结果表明: TPS最高温度为446K, 发生在紧固件连接处 ( 图12) 。由于在紧固件与支架及基板之间全部增加了隔热垫片, 使得基板的温度得到了很好的控制。

3. 4 某型号温度计算

为了能够对热试验进行快速预估, 计算模型选取了结构的一小部分 ( 相当于小片试验件) , 这样, 仅仅需要对试件的厚度进行参数化建模, 几何模型就可以通过PCL编程自动建立。

3. 4. 1 物理模型的归一化

所谓“归一化处理”就是针对热试验当前控制点 ( control points) 的区域进行模型参数化, 仅考虑试验件厚度 ( test-piece thickness) 及厚度方向的多重材料性能, 这样便可实现有限元模型的完全自动化。图13 给出有限元模型。

3. 4. 2 计算结果分析

图14 给出试验与分析结果的比较曲线, 可以看出, 两条曲线在总体上是完全可以比拟的, 需要说明的是: ①由于试验的初始温度在200 K左右, 而计算的初温为100 K, 这是最高温度差的主要原因; ②弹道参数的尾端热流急剧下降, 试验无法进行主动降温 ( 试验过程中热流控制的方法是: 当控制热流计算值小于零时, 关闭加热器, 采用自然冷却) 。而数值分析则是严格按照给定曲线进行耦合计算的。从而使得两条曲线在最后有较大差异。

4 结论

( 1) 本文对高超声速飞行器气动热计算进行了研究。采用相似性理论及参考焓等方法确定飞行器表面的气动热流, 把气动加热工程算法公式与NAS-TRAN程序在时域上进行了耦合, 开发了CASRAS程序。

( 2) 通过理论值与试验值的比较, 表明采用CASRAS程序对飞行器内部电子设备预估, TPS的隔热性能评估具有较高的精度, 能够很好地满足工程设计的需要。同时, CASRAS程序针对小片试验提出的快速耦合计算可以为热试验加热功率的选用、试验方案的制定提供技术支持。

摘要：在Window环境下, 以PATRAN为平台, 利用二次开发语言PCL (patran command language) 把气动加热工程计算方法 与NASTRAN进行实时耦合, 通过客户化、增设特定的命令和窗体等技术开发了飞行器气动加热/结构热响应耦合分析程序CASRAS (couple up aeroheating&structural response analysis program system) 。该程序可对高超声速飞行器热防护系统进行快速隔热性能评估, 对高超声速飞行器内部电子设备环境温度进行快速预估。通过物理模型归一化处理, 能够对热试验方案的制定提供技术支持, 从而将试验与理论分析有机的结合在一起, 对热结构设计和结构热试验有很高的实用价值。

关键词：高超声速飞行器,气动加热,工程算法,MSC/PATRAN/NASTRAN,结构温度场,耦合

参考文献

[1] Peter A G.Computational fluid dynamics technology for hypersonic applications.AIAA, 2003:3259

[2] Mack A.Aerothermodynamic behaviour of a generic nosecap model including thermomechanical structural effects.Aerospace Science and Technology, 2007;11 (5) :386—395

[3] 蒋军亮, 成竹, 杨志斌.飞行器燃油系统气动热试验与数值分析研究.强度与环境, 2011;38 (4) :25—31

[4] 杨炳渊, 史晓鸣, 梁强.高超声速有翼导弹多场耦合动力学的研究和进展 (上) .强度与环境, 2008;35 (5) :55—63

[5] 赵晓利, 孙振旭, 安亦然, 等.高超声速气动热的耦合计算方法研究.科学技术与工程, 2010;10 (22) :5450—5455, 5461

[6] 杨恺, 高效伟.高超声速飞行器关键部位气动热计算.计算力学学报, 2012;29 (1) :13—18

[7] 程克明, 吕英伟.飞行器持续气动加热的耦合性分析.南京航空航天大学学报, 2000;32 (2) :150—155

[8] 杨志斌, 成竹, 蒋军亮.飞行器气动加热/结构热响应耦合分析程序 (CASRAS) .中国计算机软件著作权, 2012SR082798, 2012

[9] 王勖成, 邵敏.有限单元法基本原理与数值方法.北京:清华大学出版社, 1988