铣槽式再生冷却推力室采用等肋宽的薄壁结构,其喉部自内向外分别为内壁、冷却通道和外壁,冷却通道使用肋片进行强化传热。推力室中实际发生的热量传递包括:燃气对内壁的加热、冷却剂对冷却通道的对流换热。考虑热载荷对受力的单向耦合,采用稳态热平衡方程对温度场进行传热分析,热传导的控制方程为无热源的导热微分方程,即

式中:c_p为定压比热容; ρ为材料密度; T为结构温度; t为时间; ▽为哈密顿算子; λ为材料热导率。

发动机推力室在工作过程中同时受到温度载荷和机械载荷的作用。结构的总应变ε_total由热应变ε_th、弹性应变ε_ela、塑性应变ε_pla和蠕变应变ε_cr构成,即

ε_total=ε_th+ε_ela+ε_pla+ε_cr (2)

热应变与温度的变化成线性关系,即

ε_th=α_thΔTI (3)

式中:α_th为线膨胀系数; ΔT为相对参考温度的温度变化量; I为二阶单位张量。

在分离弹性应变和塑性应变时,使用von Mises屈服准则来判断材料是否达到屈服,屈服条件表达式为

式中:f为屈服函数; s为偏应力张量; a为背应力张量,即屈服面中心; σ_y为屈服应力。

达到屈服面后塑性应变产生,塑性应变沿着屈服面梯度方向增加,即

式中:dλ为塑性乘子,即等效塑性应变增量; σ为应力张量。

计算蠕变应变时,考虑陈化理论,即Norton定律。在蠕变过程中,时效、扩散和恢复等因素影响蠕变行为,其中时间是最大的影响因素。Norton模型的表达式为

式中:ε^·_cr为蠕变速率; A为材料蠕变常数; σ_eq为Mises等效应力; n为蠕变指数; Q_c为蠕变激活能; R为气体常数。蠕变本构主要影响结构的应变构成,进而影响弹性应变和裂纹扩展的弹性应变驱动力。本文暂未考虑蠕变对材料断裂韧性的影响。

定义相场序参量c表示材料的断裂情况,如果材料完好未损坏,则c为 1; 如果存在裂纹,则c为 0。因此,变量c可以看作是损伤力学模型中的损伤参数。断裂模型中总自由能密度由弹性应变能、塑性应变能和断裂界面能组成,表达式为

ψ=ψ^ela+(1-c)²ψ^pla+ψ^c (7)

假定弹性应变中只有拉伸分量驱动裂纹扩展,将弹性应变谱分解为

ε_ela=PΛP^T (8)

式中:Λ=diag(λ₁,λ₂,λ₃)是与弹性应变张量相似的对角阵,λ₁、λ₂、λ₃为弹性应变的特征值; P为特征矩阵。

根据谱分解,弹性应变的正负分量为

弹性应变能密度分解为驱动裂纹扩展的拉伸部分ψ^ela+和不可驱动裂纹扩展的压缩部分ψ^ela-,表达式为

根据Borden的塑性断裂相场模型^[22],塑性应变能密度表达式为

ψ^pla=∫σdε_pla (11)

裂纹的相场演化方程表示为

式中:η为动力学系数; l₀为界面宽度尺寸; G为临界断裂能密度; ψ^h为驱动裂纹扩展的历史变量,即

ψ^h=max_t[ψ^ela+(ε_ela)]+ψ^pla(ε_pla) (13)

推力室结构在工作过程中会经历周期性的热应力或热应变,单次循环载荷作用下,应力水平可能远低于使材料断裂的临界应力,但是仍会经历往复循环的塑性变形。因此,需要考虑循环载荷下的相场疲劳断裂模型。通过引入疲劳退化函数,对材料的断裂韧性进行修正,使其在循环载荷下不断衰减,模拟结构在疲劳载荷下的断裂失效。采用的退化函数为

式中:α_T为疲劳阈值;(-overα)(t)为疲劳累积函数,用来计算历史变量α=‖ε_ela+ε_pla‖的累积,累积函数表达式为

h(α(·overα))为符号函数,定义为

疲劳断裂的相场演化方程表示为

不同阶段下结构的温度分布如图2所示。在热试车阶段,结构的升温比较明显,最高温度出现在内壁下表面,G点温度最高可达881.17 K。

图2 各阶段温度分布
Fig.2 Temperature distribution at different stages

选取结构上7个参考点的温度变化曲线,如图3所示。可以看出:自下而上,内壁温度最高,肋片次之,外壁温度最低。由于内壁与肋片的材料热导率高于外壁,外壁上参考点的升温速率也明显慢于其他两部分。外壁上B、C两点的温度基本一致。

图3 各参考点温度变化
Fig.3 Temperature evolution of reference points

推力室单次工作各阶段结构的应力/应变分布如图4和图5所示。在开机阶段,结构升温幅度较小,整体的压应力和压应变的数值较低。由于肋片膨胀受到内、外壁的限制,整体处于压应力水平。在热试车阶段,结构整体升温,受热膨胀,由于受到对称边界的限制,结构整体表现为受压。肋片膨胀受到外壁限制,C处会产生较高的压力。该阶段下,结构内壁下表面温度最高,G点受到左侧边界的限制,受到的压应变最大,最高达到1.34%。在后冷阶段,结构弹性变形恢复,内壁由于压载较大,降温后有残余的塑性变形。内壁下表面F点和G点分别表现为残余拉应变和残余压应变。F点受到拉伸载荷,残余应变约为1%。在关机阶段,结构的温度分布与后冷阶段相差不大,结构最终的受载情况和残余应变与后冷阶段基本一致。

图4 各阶段应力分布
Fig.4 Stress distribution at different stages

图5 各阶段应变分布
Fig.5 Strain distribution at different stages

根据应力和应变分布图,确定D、E、F、G这4点在工作过程中受到的局部应力-应变较大,为服役危险点。这4点的应力、应变曲线如图6所示。4个点在循环中都经历压应力变为拉应力,D、E、G点始终受到压应变,F点最终由压应变转变为拉应变。

对10个工作循环下的推力室喉部结构进行仿真。根据单次循环应力应变分析,结构的危险点位于内壁下表面的F和G点,循环仿真后重点关注F和G点的应力应变情况。F点和G点10次循环的应力应变曲线分别如图7和图8所示。随着循环数的增加,F、G两点表现为应力控制的非对称循环加载,并且每次循环结束残余应变值不断增加。其中,F点为拉伸应变,G点为压缩应变。F、G两点残余应变发展曲线如图9所示。

图6 单次工作中4个特征点的应力-应变曲线
Fig.6 Stress-strain curves of four points during single cycle

图7 F点10次循环应力-应变响应
Fig.7 Stress-strain response of point F in 10 cycles

图8 G点10次循环应力-应变响应
Fig.8 Stress-strain response of point G in 10 cycles

图9 F、G点10次循环残余应变的变化
Fig.9 Residual strain evolution at F and G in 10 cycles

由于单次循环结构承受的最大应力小于材料断裂破坏的临界应力,结构中的断裂序参量增长较小。循环加载并计算结构的断裂序参量,直至结构中出现序参量值达到0.95则停止计算。最终发现,91个循环后结构失效,计算停止,F点断裂序参量最先达到临界值c=0.95,表明结构的疲劳寿命为91次循环,失效点位于内壁下表面F点。

单次循环后的序参量分布如图 10所示。图 10表明,结构中内壁下表面中心点F序参量最大,在单次工作结束后,F点的损伤最大,但此时序参量的值还处于很低的水平,对结构的承载能力影响不大,不会使应力衰减,仍可以进行多次循环工作。结构在内壁下表面受到的温度载荷最高,并且F点由于对称边界受到约束,在温度恢复常温后,受到的拉应力最大,所以序参量演化最快。根据应力分布云图,在结构变截面处,即C、D点会出现应力集中。C点受到的温度载荷较低,所以温度载荷恢复后受到的载荷水平较低,没有出现损伤。D点受到的温度载荷较高,相场值为6×10^-7,但与F点相比仍较低,因此结构并未在变截面处出现裂纹。

10次循环后的序参量分布如图 11所示,结果与单次循环相似,仍在F点序参量值最高。

结构中的断裂序参量演化如图 12所示。随着循环数增加,下表面的断裂序参量值不断上升,并且EF边界明显快于内壁其他区域,从F点到E点从大到小基本呈线性分布。

图 10 单次工作后相场序参量分布
Fig.10 Distribution of phase-field order parameter after 1 cycle

图 11 10次工作后相场序参量分布
Fig.11 Distribution of phase-field order parameter after 10 cycles

图 12 相场序参量演化
Fig.12 Evolution of phase-field order parameter

结构失效破坏时的径向位移如图 13所示。疲劳破坏时,EF边界明显向下凹陷,内壁不断变薄,最终在F点处断裂失效。

图 13 结构破坏时的径向位移
Fig.13 Radial displacement at failure