考虑气动加热瞬态效应的复合材料壁板超音速LCO特性研究

刘雅婷, 段静波,2, 徐步青,2, 高艺航

(1. 石家庄铁道大学 工程力学系,石家庄 050043; 2. 石家庄铁道大学 河北省工程力学基础学科研究中心,石家庄 050043; 3. 国防科技大学 空天科学学院,长沙 410073; 4. 北京航天系统工程研究所,北京 100076)

超音速飞行器在高速飞行过程中会产生大量的热量。严重的气动加热会显著改变飞行器蒙皮结构温度场,导致恶劣的热环境,从而可能导致飞行器力学性能、气动弹性大变形和失稳边界发生变化,进一步降低飞行器的安全性[1]。为了解决这些问题,许多研究者对飞行器在热环境下的气动弹性行为进行了研究。通常,文献研究中主要采用稳态温度场,忽略了气动加热产生的瞬态温度场对力学性能和结构稳定性的影响。自吴志刚等[2]提出分层求解的方法之后,研究者们广泛采用这种单向耦合的方法来研究瞬态气动加热对气动弹性行为的影响[3]。Tran等[4]研究了F-16二维机翼在气动加热和气动力作用下的温度分布和时域响应。Mcnamara等[5]将气流与结构之间的传热结合起来,讨论了气动加热对结构稳定性的影响。Miller等[6]提出了一种流固热松耦合模型,并证明了压力、位移插值和热流密度的估算保持了壁板运动参数和温度的二阶精度。Huo等[7]采用基于活塞理论和薄激波层的时适应松耦合方法研究了C/SiC复合材料板的气热弹性行为。Chen等[8]采用时间自适应高超声速热-流-固松耦合模型用于热气动弹性响应分析,重点研究了低展弦比机翼的自振和热模态变化。Khalafi等[9]利用三阶活塞理论和广义微分求积法分析了功能梯度板的热气动弹性特性。Ye等[10]研究了高超声速机翼热气动弹性变形对机翼扭转角和下表面压力的变化规律的影响。Kamali等[11]采用有限元求解器计算热解和结构解模拟了热气动弹性耦合,并将计算结果与试验数据进行了对比。

随着高超声速飞行器的出现,气动和结构的耦合作用愈加显著,气动弹性变形对瞬态气动加热的影响不容忽视[12]。Thornton等[13]采用有限元法分析了不锈钢薄板的流-热-固耦合作用。结果表明,气动弹性变形会影响加热速率和温度分布。Li等[14]采用热气动弹性双向耦合方案,研究了功能梯度曲板不同拱高对颤振特性的影响。Wan等[15]采用双向耦合的静态热气动弹性分析方法,得到了高超声速机翼的热流密度、气动力、热场和变形。研究发现,由温度引起的变形对气动力分布有很大影响。Ye等[16]开发了松耦合热气动弹性分析系统,采用单向耦合和双向耦合计算方法研究了热气动弹性变形对高超声速进气道性能和流场的影响。Huang等[17]利用模态降阶法方法和紧耦合模式进行了热气动弹性的计算。

基于此,本文建立了考虑气动加热瞬态温度效应的复合材料壁板热气动弹性模型。在模型验证后,重点讨论了气动加热瞬态温度场中斜激波、气动压力、传热系数、初始干扰力等关键因素对复合材料壁板非线性热气动弹性LCO响应的影响规律,为超音速飞行器复合材料壁板的工程设计提供技术支持。

1 热气动弹性建模

以超音速气流中的复合材料层合板为研究对象,如图1所示。壁板沿长度x方向为a,沿长度y方向为b,壁板的厚度为h。建模过程中,考虑超音速气流经过楔形体产生的斜激波、气动加热产生的壁板内瞬态温度场,以及壁板结构的大挠度变形等。

图1 复合材料壁板示意图Fig.1 Schematic diagram of the composite panel

1.1 激波参数和气动力

超音速气流经过楔形体会被压缩,在楔形体前缘产生斜激波,因此,复合材料壁板边界层气体参数会发生变化。壁板表面激波后气体参数与来流气体参数关系如下[18]

(1a)

(1b)

(1c)

(1d)

式中:β为斜激波角;ρ∞,V∞,Ma∞和T∞分别为来流密度、来流速度、来流马赫数和温度;ρe,Ve和Te为激波后气流参数。

当来流马赫数Ma∞>1.6时,壁板表面沿x方向的气动力分布根据一阶活塞理论获得[19]

(2)

(3)

式中,D110为铺层角为0°时壁板弯曲刚度矩阵的第一个元素。

1.2 壁板的本构关系

复合材料壁板在高速气流中可能会产生非线性变形,因此,基于Von karman大变形假设,壁板的位移-应变关系可以表示为

(5)

式中:u0,v0和w0分别为中性层在x,y和z方向的位移;ε0为中面的总面内应变向量;κ为板的曲率;γ为剪切应变;θx和θy分别为中性面绕x轴和y轴的转角。

同时,考虑瞬态非均匀温度场对复合材料壁板热应力和力学性能的影响,结构的应力-应变关系定义为

(6)

根据层合板的层合效应,n层对称复合材料壁板的内力可以由下式进行计算

(7a)

Fs=Asγ

(7b)

(7c)

式中:N,M和Fs为壁板的膜力、弯扭矩和剪力;A,D,As分别为刚度矩阵、弯曲刚度矩阵和剪切刚度矩阵;NΔT和MΔT为温度变化产生的膜力和弯矩。

1.3 气动加热模型

考虑壁板上表面对外界的气体辐射效应,根据牛顿冷却公式,可以得到气动加热产生的复合材料壁板表面的热流

(8)

式中:α为热传递系数;Tw为壁温;Taw为绝热壁温可以表示为

(9)

α=ρ*VeSt*cp

(10)

式中:Ve为边界层气流速度;cp为定压比热容;ρ*为参考温度下的边界层气流密度;St*为参考温度下的斯坦顿数。参考温度可以由下式进行计算

(11)

1.4 壁板热传导模型

由于气动加热产生的热流会在气流方向发生变化,因此,考虑壁板热传导沿厚度方向(z方向)和气流方向(x方向)的温度差异,采用二维热传导模型分析,其控制方程为

(12)

式中:ρ和c分别为材料密度和比热容;kx和kz为材料在x和z方向的热传递系数。

根据壁板上表面与空气的对流换热边界条件可得

(13)

式中:k为壁板与空气之间的对流换热系数;q为气动加热的热流密度,可以由式(8)求得。

将边界条件代入控制方程,采用有限差分法计算瞬态温度场。沿x方向和z方向进行网格划分,z方向网格数与壁板层数一致。具体差分格式可以表示为

(14)

式中,i和j分别为x和z方向的节点编号。

1.5 气动弹性模型

复合材料壁板结构变形产生的虚应变能δUM、温度变化产生的虚应变能δUT、壁板的虚动能δT,气动力和阻尼力产生的外力虚功δW,具体如下

(15a)

(15c)

δW=∬AΔpδwdA

(15d)

式中,ρ为材料密度。

然后,利用四边形四节点单元对复合材料层合板进行离散。单元的位移向量可表示为

u=Nuδe

(16)

式中,u,δe和Nu分别为各单元的未知位移向量、节点位移向量和形函数矩阵。具体表达式为

u=[u,v,w,θx,θy]

(17a)

δe=[u1,v1,w1,θx1,θy1…u4,v4,w4,θx4,θy4]

(17b)

(17c)

根据哈密顿变分原理

(18)

将式(16)代入式(18),通过变分可以得到复合材料壁板的热气动弹性分析的非线性运动微分方程

(19)

式中:M为质量矩阵;Ca为气动阻尼矩阵;K0,Knon,Ka和KT分别为线性刚度矩阵、非线性刚度矩阵、气动力矩阵和热应变刚度矩阵;fΔT为温度引起的壁板总体载荷列阵。对应矩阵的具体表达式见附录A。

在对式(19)求解时,首先进行模态降阶,缩减后的运动方程可以表示为

(20)

[F({q}i+1)]{q}i+1-{f}i+1=0

(21)

其中,

(22b)

最后,由于运动方程中具有非线性项,采用牛顿迭代法对式(21)进行求解。

2 气动加热和气动弹性振动耦合计算

由于计算中考虑气动加热的瞬态效应,因此,气动加热和气动弹性振动需在时间增量步内分别求解,再进行瞬态温度响应和气动弹性振动响应的数据交换,依此进行气动加热与气动弹性振动耦合计算。具体而言,首先,根据t时刻温度得到材料参数,进而计算t～t+Δt时刻壁板的气动弹性响应。其次,根据壁板的变形量可以获得壁板表面t+Δt时刻的气动力分布。同时,可以计算t+Δt时刻与壁板变形有关的气动参数。在气动加热分析中,利用更新后的气动参数可以得到t+Δt时刻的温度场。接着,用同样的方法可以得到下一时间步壁板的气动弹性响应,直至最后。计算流程图如图2所示。

图2 热气动弹性耦合计算流程图Fig.2 Flow chart of aerothermoelastic coupling calculation

3 算例验证与结果分析

3.1 算例验证

为了验证现有算法和程序,首先采用本文的程序计算壁板的瞬态温度场,并与Culler等研究进行对比,验证壁板中心点温度随时间的变化。壁板几何参数和材料物性参数均与Culler等研究中的一致,计算结果如图3所示。可以看到,本文计算结果与文献结果吻合很好。这说明本文程序关于壁板瞬态温场计算的程序是正确的。

图3 壁板中心点瞬态温度场的变化Fig.3 Variation of transient temperature field at the center point of the composite panel

其次,通过稳态温度场中四边简支矩形壁板的颤振验证本文关于气动弹性模型计算程序的正确性。复合材料壁板的几何尺寸和材料参数与欧阳小穗等研究中的相同。采用欧阳小穗等研究中的两种铺层结构。将本程序得到的结果与欧阳小穗等的研究结果进行比较,如图4所示。这说明本文关于壁板气动弹性计算程序是正确的。

图4 LCO振幅随无量纲动压的变化Fig.4 Variation of LCO amplitude with dimensionless dynamic pressure

3.2 气动加热瞬态温度场中壁板时域响应特性

以四边固定的矩形复合材料层合壁板为研究对象,研究其在考虑气动加热瞬态温度效应下的颤振响应特性。壁板尺寸及材料参数,如表1所示。考虑壁板上表面气动加热以及壁板上表面对外界的气体辐射效应,气动加热产生的热流为q,壁板左右两侧为绝热边界条件,下表面为恒温边界条件,T=300 K,如图5所示。

表1 壁板尺寸和材料参数Tab.1 Panel dimensions and material parameters

图5 热边界条件示意图Fig.5 Schematic diagram of thermal boundary conditions

为了便于比较,本文同时给出了瞬态温度场和稳态温度场中复合材料壁板的颤振时域响应和温度响应,如图6所示。由图6(a)中可以看出,复合材料壁板在考虑气动加热瞬态温度效应下的LCO幅值与稳态温度下的计算结果基本一致。这说明通过稳态温度场近似预测瞬态热环境下壁板颤振特性是合理的。然而,从壁板整个时间响应历程来看,受温度场瞬态特性的影响,复合材料壁板的气动弹性振动时域响应不同于稳态温度场情形,主要体现在振动的早期阶段。从图6(a)可以看到,稳态温度场中的复合材料壁板振动响应很快就会进入LCO状态,时间过程相对较短。但是考虑了气动加热产生的温度场瞬态效应后,复合材料壁板在进入LCO状态前有一个相对较长的振动历程,而且壁板在进入LCO状态前其振动位移是先收敛后增大,见图6(a)所示。主要原因是气动加热下壁板内温度场存在一个温度逐渐升高过程,如图6(b)所示。值得注意的是,计算得到的瞬态温度场最终将是随时间持续微幅振荡变化的,这主要是由壁板最终处于LCO状态导致的。此外,图6(c)和图6(d)分别给出了复合材料壁板上一个周期内不同时刻的位移振型图和温度场分布图。可以看到,壁板的周期变形影响其表面的气动力分布,进而影响气动加热热流,最终导致壁板温度场在长度方向以图6(d)所示形态微幅波动。

图6 复合材料壁板气动弹性振动时域响应与瞬态温度响应(λ=166)Fig.6 Aeroelastic vibration time-domain and transient temperature responses of the composite panel(λ=166)

3.3 主要参数影响分析

进一步讨论斜激波、传热系数、初始扰动力等因素对考虑气动加热瞬态温度效应下壁板气动弹性LCO响应和瞬态温度响应的影响。

首先研究斜激波对复合材料壁板在瞬态温度场中LCO响应的影响。如图7所示给出了楔形体半锥角为10°时复合材料壁板的振动时域响应和瞬态温度响应。从图7(a)中可以发现,斜激波会使壁板的LCO幅值显著增加,也就是说,斜激波的产生会大大降低壁板的临界颤振动压,降低壁板的颤振稳定性边界。图7(b)给出了斜激波作用下壁板的瞬态温度场。从图中可以看出,斜激波的产生也会使壁板的温度场显著增大,这也是导致壁板临界颤振动压变化的原因所在。而且,从图8中不同动压下瞬态温度场响应可以看出,动压越大,斜激波对壁板瞬态温度场影响越明显,进而对壁板的颤振影响也会越大。

图7 激波前后复合材料壁板气动弹性振动时域响应与瞬态温度响应(λ=157)Fig.7 Aeroelastic vibration time-domain and transient temperature responses of the composite panel ahead of and behind shock wave(λ=157)

图8 不同动压下激波前后复合材料壁板瞬态温度响应Fig.8 Transient temperature responses of the composite panel ahead of and behind shock waves at different dynamic pressures

图9～图11分别给出了复合材料壁板在长度方向和厚度方向的温度场分布。从图9中不难发现,壁板前缘受气动加热的影响最严重,温度梯度较大,而沿着壁板的长度方向温度梯度逐渐减缓。另外,瞬态温度场会随动压的增大而增大,但相比于激波前,激波后动压对瞬态温度场的影响变得更为明显。图10所示动压为λ=164时,复合材料壁板内不同位置的瞬态温度场随时间的变化情况,可以看出,壁板内瞬态温度场随着时间变化逐渐趋于等幅波动,但不同位置处的温度场波动幅值不同。图11进一步给出了壁板厚度方向的温度场分布。可以看到,壁板沿厚度方向会存在非均匀的温度梯度。而且,通过动压λ=155和λ=164两组温度分布结果对比,可以得出,随着动压增大,壁板厚度方向的温度变化明显增大,这将会影响壁板的临界颤振动压。因此,考虑壁板厚度方向的温度非均匀分布是十分必要的。

图9 不同动压下激波前后复合材料壁板上表面温度场分布Fig.9 Temperature field distribution on the upper surface of the composite panel ahead of and behind shock waves under different dynamic pressures

图10 激波后复合材料壁板不同位置的瞬态温度响应(λ=164)Fig.10 Transient temperature responses of the composite panel at different positions behind shock wave(λ=164)

图11 激波后复合材料壁板沿厚度方向的温度场分布Fig.11 Temperature field distribution along thickness of the composite panel behind shock wave

图12为传热系数对气动加热瞬态温度场中复合材料壁板LCO响应的影响。其中,传热系数k分别取值为20 W/m/K,40 W/m/K,60 W/m/K三种情况。从图中可以看出,随着传热系数从60 W/m/K降低到20 W/m/K,壁板的LCO幅值显著增大,这意味着壁板的临界颤振动压会明显下降。原因在于传热系数较大时,气动加热会在复合材料壁板内产生变化相对较小的瞬态温度场,如图12(b)所示,进而影响了壁板的临界颤振动压。

图12 不同传热系数下复合材料壁板气动弹性振动时域响应与温度响应(λ=180)Fig.12 Aeroelastic vibration time-domain and transient temperature responses of the composite panel at different heat transfer coefficients(λ=180)

图13给出了不同初始扰动力作用下考虑气动加热瞬态温度效应的复合材料壁板振动时域响应和瞬态温度响应。从图13(a)中可以看出,初始扰动力不会影响复合材料壁板的颤振临界动压,但是在时间历程上会影响复合材料壁板初始阶段的收敛振动,以及壁板进入LCO的时间。当初始扰动力较大时,壁板初始阶段的收敛振动幅值变化显著,且壁板进入LCO的时间较早。当初始扰动力较小时,壁板初始阶段的收敛振动幅值变化不明显,壁板进入LCO的时间相对较晚。从图13(b)中可以看到,不同初始扰动力作用下,壁板的瞬态温度场也会发生变化,这其中的原因主要是初始扰动影响了壁板的气动弹性振动响应,进而影响了壁板的瞬态温度场。

图13 不同初始扰动力下复合材料壁板气动弹性振动时域响应与瞬态温度响应(λ=157)Fig.13 Aeroelastic vibration time-domain and transient temperature responses of the composite panel at different initial disturbance forces(λ=157)

4 结 论

本文建立了气动加热瞬态温度效应的复合材料壁板气动弹性模型,给出了壁板的气动弹性时域响应和瞬态温度响应特性,并在考虑气动加热瞬态温度效应下,重点讨论分析了斜激波、来流动压、传热系数、初始干扰力等对复合材料壁板的气动弹性LCO响应的影响。主要结论如下:

(1)考虑了气动加热产生的温度场瞬态效应后,复合材料壁板的气动弹性振动时域响应早期阶段不同于稳态温度场情形,壁板在进入LCO状态前有一个相对较长的振动历程,且壁板振动位移是先收敛后增大,计算得到的温度场最终是随时间持续微幅振荡变化的。

(2)斜激波会使壁板的温度场升高,并且相比于激波前,激波后动压对瞬态温度场的影响更为明显,这进而导致复合材料壁板的LCO幅值显著增加,大大降低壁板的临界颤振动压,降低壁板的颤振稳定性边界。

(3)传热系数的变化对壁板的LCO振幅有显著的影响,传热系数减小,会导致壁板内瞬态温度场变化加剧,从而通过影响壁板刚度使壁板的LCO振动幅值增大,降低壁板的临界颤振动压。

(4)初始扰动力不会影响复合材料壁板的颤振临界动压,但是在时间历程上会影响复合材料壁板初始阶段的收敛振动,以及壁板进入LCO的时间,也会影响瞬态温度场振动发生的时间。

附录A

(A.1)

(A.2)

(A.3)

(A.4)

(A.5)

(A.6)

(A.7)

其中,

(A.8)

(A.9)

(A.12)