引言

在核电管路系统中，T型管结构主支管流体的温差会使得管道内流体在掺混时出现剧烈的温度波动，引发管道产生热应力，导致疲劳破坏的发生。改变主支管流体的动量比可以使得管道内流体的掺混流型发生变化，进而改变管道内温度的波动区域，从而改善管道的热疲劳情况。Hosseini等^[1]通过实验探究了在T型管主管上游有无弯管的情况下，不同主支管流体动量比M_R对管道内流体流动的影响，并根据动量比大小将流型分为壁面流(M_R > 4)、再附着流(1.35 < M_R < 4)、偏转流(0.35 < M_R < 1.35)和冲击流(M_R < 0.35)。单祖华等^[2]通过数值模拟的方式，将重整化群(RNG)k-ε模型、Speziale-Sarkar-Gatski(SSG)雷诺应力模型和大涡模拟(LES)模型的计算结果与实验结果进行对比，证明了LES模型在冷热流体掺混区域能更好地表征流体的实际流动状态。卢涛等^[3]利用LES模型探究冷热流体掺混过程中流体速度波动与流体动量比变化之间的规律，获得了不同流型下的速度波动强度及波动位置。

对于流固耦合问题及热疲劳评估问题，国内外学者也进行了深入的研究。胡丽娜等^[4]为探究在启、停堆工况下管道的热应力状态，在耦合计算时对三通管道施加突变的温度载荷，得到管道热疲劳寿命与支管入射角度之间的关系式。Hannink等^[5]通过依次导入不同时刻下的管道温度场，计算出多个时刻下的管道热应力状态，并对计算结果进行了分析。Kumar等^[6]研究了热分层引起的热应力及其对管道疲劳寿命的影响。Guo等^[7]对水平T型管进行模拟试验和分析，设计了一个监测管路壁温的程序，能够通过测量外壁温度来监测水平T型管冷热流体混合时的壁温波动。Garrido等^[8]通过生成随机的管内壁面温度，使之与流场计算得到的随机温度载荷类似，进而探究在非稳态载荷下T型管的热应力波动。

在目前的研究中，评估不同动量比对管道热疲劳影响的工作较少，而T型管内流体的搅混流型对管道的热疲劳又有直接影响。本文通过LES模型对T型管内流体的搅混状态进行瞬态模拟，在保证使用同一套计算网格模型的前提下，通过CFD-post软件导出固体温度载荷及内壁面压力载荷，基于workbench平台进行流固耦合计算，并评估管道的热疲劳情况，为研究T型管内动量比对管道热疲劳的影响提供了一定的数据支撑。

1 数学模型 1.1 LES模型

本文将k-ε湍流模型计算得到的流场稳态解作为初始场，利用LES模型进行流场的瞬态计算。大涡湍流模型利用亚格子(SGS)尺度模型，以更好地捕捉瞬态流场中大尺度涡与小尺度涡之间的关系，SGS模型通常采用涡-黏度模型表示

$ \tau_{i j}-\frac{\tau_{k k} \delta_{i j}}{3}=-2 \mu_{\mathrm{t}} \bar{S}_{i j} $

(1)

式中，τ为亚格子尺度应力，δ_ij为克罗内克符号，μ_t为SGS湍流黏度，$\bar{S}_{i j} $为求解的应变率张量。

亚格子雷诺应力τ_ij可定义为

$ \tau_{i j}=\rho \overline{u_{i} u_{j}}-\rho \bar{u}_{i} \bar{u}_{j} $

(2)

式中，ρ为流体密度，u为流体流动速度。

在亚格子Smagorinsky-Lily模型中，SGS湍流黏度μ_t可以定义为

$ \mu_{\mathrm{t}}=\rho L_{S}^{2}|\bar{S}|^{2} $

(3)

式中，L_S表示亚格子尺度的混合长度，计算式为

$ L_{S}=\min \left(k d, C_{\mathrm{S}}V^{\frac{1}{3}}\right) $

(4)

其中，k为冯卡门常数，取0.41，d为第一层网格与壁面间距离，C_S为Smagorinsky常数，取0.1，V为计算单元体积。

1.2 应力计算理论模型

本文的应力计算基于ANSYS软件的有限元(FEM)分析理论进行，理论模型为

$ \boldsymbol{\sigma}=\boldsymbol{D} \times\left(\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{0}\right) $

(5)

式中，D为弹性矩阵，ε为节点形变，可定义为

$ \boldsymbol{\varepsilon}=\boldsymbol{B} \cdot \boldsymbol{\delta}^{\mathrm{e}} $

(6)

式中，B为应变矩阵，δ^e为节点位移，可定义为

$ \boldsymbol{K} \cdot \boldsymbol{\delta}^{\mathrm{e}}=\boldsymbol{Q}_{\mathrm{T}} $

(7)

式中，K为刚度矩阵，Q_T为载荷。

ε₀为由温度变化引起的变形，可定义为

$ \boldsymbol{\varepsilon}_{0}=\alpha T\left[\begin{array}{llllll} 1 & 1 & 1 & 0 & 0 & 0 \end{array}\right]^{\mathrm{T}} $

(8)

式中，α为材料的线膨胀系数，T为温度变化量。

2 计算模型

T型管的管道结构如图 1所示。主、支管管径分别为D_m与D_b，主、支管流体流速分别为V_m与V_b，主管流体动量为M_m，支管流体动量为M_b，主、支管流体动量比M_R计算式如式(11)所示，其中，ρ_m、ρ_b分别为主、支管流体的密度。

$ M_{\mathrm{m}}=\rho_{\mathrm{m}} \cdot D_{\mathrm{m}} \cdot D_{\mathrm{b}} \cdot V_{\mathrm{m}}^{2} $

(9)

$ M_{\mathrm{b}}=\rho_{\mathrm{b}} \cdot \frac{{\rm{ \mathsf{ π} }} D_{\mathrm{b}}^{2}}{4} \cdot V_{\mathrm{b}}^{2} $

(10)

$ M_{\mathrm{R}}=\frac{M_{\mathrm{m}}}{M_{\mathrm{b}}}=\frac{4 \rho_{\mathrm{m}} D_{\mathrm{m}} V_{\mathrm{m}}^{2}}{{\rm{ \mathsf{ π} }} \rho_{\mathrm{b}} D_{\mathrm{b}} V_{\mathrm{b}}^{2}} $

(11)

本文以液态水为工作介质，主管内流体为热流体，支管内流体为冷流体。T_m与T_b分别为主、支管流体温度，为定量分析管道内流体温度的波动情况，对流体温度进行无量纲归一化处理，流体无量纲温度T_i^*定义为

$ T_{i}^{*}=\frac{T_{i}-T_{\text {cold }}}{T_{\text {hot }}-T_{\text {cold }}} $

(12)

式中，T_i为流体温度，T_hot与T_cold分别为热流体与冷流体温度。

为表征流体时均温度大小，定义无量纲时均温度T_mean^*为

$ T_{\text {mean }}^{*}=\frac{1}{n} \sum\limits_{i=1}^{n} T_{i}^{*} $

(13)

流体温度的波动情况由无量纲均方根温度T_rms^*表示，其定义式为

$ T_{\mathrm{rms}}^{*}=\sqrt{\frac{1}{n} \sum\limits_{i=1}^{n}\left(T_{i}^{*}-T_{\text {mean }}^{*}\right)^{2}} $

(14)

在本文中，主支管流体温差保持在100 ℃不变，通过改变主支管流体管径比或速度比的方式来改变流体的动量比，具体的计算工况如表 1所示。表中，case1、case2及case3工况保持支管入口质量流量一致，改变主支管管径比与速度比，以实现主支管流体动量比的改变，使得T型管内流体的搅混流型由偏转流型转变为壁面流型；case4工况则通过进一步增大主支管流体速度比，来研究动量比对管道热疲劳的影响。

管道内监测点布置方案如图 2所示，共布置如图所示的P1~P9共9个监测面，面与面之间间隔100 mm，每个监测面上有L1~L8共8个监测方向，每个方向布有内壁面监测点以及据内壁面1 mm处的流体监测点。同时，管道主支管交汇处内壁面与主管相交的相贯线为PA监测线，其上布置PA1~PA4共4个应力监测点，与支管相交的相贯线为PB监测线，其上布置PB1~PB4共4个应力监测点，用以监测主支管交汇区热应力的波动情况。为便于后续分析，以监测点P3-L5-I/F为例，其表示在P3监测面L5方向上，内壁面/近壁面流体的监测点。

为保证计算的准确性，对管道模型进行结构化网格划分。为满足大涡模拟的计算要求，划分网格时控制边界层第一层的网格高度，保证y+值小于5，同时进行网格无关性验证。网格划分方案如表 2所示。

图 3为在case4工况条件下，P3-L5-I监测点的网格无关性验证结果。如图所示，相较于网格数量最少的mesh1网格的数值计算结果，网格数量更多的mesh2与mesh3网格的数值计算结果相差更小。因此mesh2网格能在满足计算精度的同时提高计算效率，故选择该网格作为本文数值模拟的网格划分方案。

3 数值计算及结果分析 3.1 CFD计算

本文利用Fluent软件对T型管道进行数值模拟计算，通过标准k-ε湍流模型计算得到收敛的稳态结果，以该结果作为LES模型的初始场进行瞬态计算。求解选用couple耦合求解器及二阶迎风格式，主支管入口均为速度入口(velocity inlet)，出口为自由出流(outflow)，管道外壁面为绝热壁面，内壁面为耦合无滑移边界。计算总时长为20 s，时间步长为0.005 s，数据监测频率为200 Hz。由于瞬态计算结果在10 s后逐渐趋于稳定，故取15~20 s时间内的数据进行对比分析。

图 4为T型管在不同动量比工况下支管入口截面的温度分布云图。由图可以看出，随着主支管流体动量比逐渐增大，冷热流体混合区域逐渐向支管方向发展。在4例工况中，case1与case2工况管道内的搅混流型为偏转流，两种工况管道壁面的温度几乎没有波动，case3与case4为壁面流工况，其中case3工况管道温度波动范围最大。

为表征冷热流体掺混区域流体的流动状态，比较不同工况下P3截面位置的主管速度矢量分布图，结果如图 5所示。由图 5可以看出，在搅混流型为偏转流条件下，流体在主管内的掺混区域随着动量比的增大，从case1工况的L6方向发展到case2工况的L5方向；而在搅混流型为壁面流的case3与case4工况条件下，流体掺混发生在主管内部的L4方向上，该现象是由在主支管交汇处支管上游弯管产生的二次流造成的。对比近壁面流体速度矢量的大小可以发现，壁面流工况下近壁面流体的速度及掺混区域均明显大于偏转流工况。

由于支管上游弯管的影响，在弯管产生的二次流的作用下，掺混区域内流体从L6方向向L5方向旋转，故以L5方向为例，比较不同动量比工况下，近壁面流体的无量纲时均温度与无量纲均方根温度的分布曲线，如图 6、图 7所示。

图 6为不同动量比工况下，L5方向近壁面流体在各监测点上的无量纲时均温度大小。如图 6所示，偏转流工况(case1与case2)下在近壁面位置的无量纲时均温度均近似为1，即冷热流体掺混主要发生在流体中心，近壁面流体温度保持不变；而壁面流工况(case3和case4)下近壁面流体的温度沿主管轴向先下降再升高，在P3截面位置最低，其中case3工况P3截面流体监测点的无量纲时均温度最低，为0.42。

图 7为不同动量比工况下，L5方向近壁面流体在各监测点上的无量纲均方根温度大小。如图 7所示，对于偏转流工况，case1与case2工况近壁面流体的无量纲均方根温度近似为0，即偏转流工况近壁面位置流体的温度波动很小。而在壁面流工况中，case3与case4工况近壁面流体的无量纲均方根温度明显增大，且case3工况温度的波动程度要明显大于case4工况，其主要原因为，在壁面流工况中，随着动量比逐渐增大，冷热流体掺混区域会逐渐靠近主管上壁面，使得近壁面流体温度的波动区域减小。无量纲均方根温度最大位置位于case3工况P2截面的流体监测点处，为0.18。

本文中，case1、case2及case3工况主支管流体温差及流量均未发生改变，而case4工况相对于case3工况仅改变流体的速度比，根据式(11)中对动量比的定义，影响这3例工况动量比变化的主要因素为主支管流体的速度比，随着速度比的增大，主支管流体间动量比增大，冷热流体的掺混区域逐渐靠近主管上壁面。

3.2 CFD-FEM流固耦合计算

在T型管结构中，最易发生疲劳破坏的位置主要出现在主支管相贯区以及冷热流体掺混区域。在本文中，冷热流体发生掺混的位置位于T型管相贯处，故相贯区的PA1~PA4、PB1~PB4这8个监测点是应力计算时要重点关注的位置。

根据Fluent数值模拟的计算结果，将15~20 s时间内管道的瞬时温度载荷与内壁面瞬时压力载荷导出，通过workbench平台加载到固体计算域中，固体计算域的约束条件及所施加的载荷如表 3所示，管道材料选用316L不锈钢。

本文采用第三强度理论将管道温度载荷与压力载荷带来的多轴应力转化为单轴应力，第三强度理论的当量应力σ定义为

$ \sigma=\sigma_{1}-\sigma_{3} $

(15)

式中，σ₁为第一主应力，σ₃为第三主应力。

为定量分析T型管相贯区的应力大小及应力波动情况，计算相贯区内8个应力监测点在15~20 s内的时均应力σ_mean^*与均方根应力σ_rms^*，其计算公式如式(16)、(17)所示。

$ \sigma_{\text {mean }}^{*}=\frac{1}{N} \sum\limits_{i=1}^{N} \sigma_{i} $

(16)

$ \sigma_{\mathrm{rms}}^{*}=\sqrt{\frac{1}{N} \sum\limits_{i=1}^{N}\left(\sigma_{i}-\sigma_{\text {mean }}^{*}\right)^{2}} $

(17)

式中，N为总的计算时间步数，σ_i为每个时间步计算得到的监测点应力。

相贯区各监测点的时均应力大小如图 8所示。在相贯区与主管相交的PA监测线上，除位于上游的PA2监测点处时均应力不随T型管内流体搅混流型的变化而出现明显改变外，其余各监测点在壁面流工况下的时均应力均大于偏转流工况下的时均应力。其主要原因为壁面流工况下近壁面流体的温度波动大，使得其内壁面温度波动增大，时均应力大于偏转流工况，其中，case4工况PA3监测点的时均应力最大，为215.82 MPa。由于这4例工况冷热流体掺混区域主要集中在主管下游，位于上游的PA2监测点的温度波动较小，使得该位置的应力主要受主支管温度梯度大小的影响，本文中各工况温差均为100 ℃，温度梯度分布相近，故在PA2点的时均应力受流型的影响小。在相贯区与支管相交的PB监测线上，由于T型管内没有湍流穿透现象发生，支管温度波动小，其时均应力主要受温度梯度的影响，故各工况监测点的时均应力大小基本一致。

图 9为各工况相贯区监测点的均方根应力图。如图 9所示，case1与case2工况在PB监测线上的应力波动大于PA监测线，case3工况下PA监测线应力波动剧烈，其中PA1监测点的应力波动最大，为31.66 MPa；case4工况应力波动程度相比于其他各工况波动程度较小，仅在PA3位置应力波动达到6.77 MPa。即在case3工况下，相贯区靠近主管位置的应力波动最大，其原因是流体的温度波动主要发生在这一位置，从而使得壁面温度的波动更剧烈，应力波动也随之增大。

3.3 热疲劳分析

根据有限元分析方法得到各点应力的波动信号后，对应力信号进行预处理，保留信号的峰值与谷值，并对预处理信号进行雨流计数分析^[9]，统计在15~20 s时间段内，不同平均应力下各交变应力幅的循环次数，得到雨流矩阵N。

由于计算得到的应力循环载荷为非对称循环载荷，故本文利用Goodman平均应力修正曲线对应力载荷进行修正，得到其等效的对称循环应力载荷，如式(18)所示

$ \sigma_{\mathrm{a}}=\sigma_{-1}^{*}\left(1-\frac{\sigma_{\mathrm{m}}}{\sigma_{\mathrm{b}}}\right) $

(18)

式中，σ_a为交变应力幅，σ_m为平均应力，σ_b为材料强度极限，σ_-1^*为等效对称循环应力，本文中管道材料为316L不锈钢，其强度极限为520 MPa。

根据美国核管理委员会制订的NUREG/CR-6909文件^[10]中的奥氏体不锈钢疲劳寿命S-N曲线，求得在不同交变应力下管道的许用循环次数矩阵N_f。雨流矩阵与许用循环次数矩阵相除即可求得管道的疲劳损伤矩阵。管道监测点上的疲劳损伤矩阵D定义为

$ \mathit{\boldsymbol{D}} = \frac{\mathit{\boldsymbol{N}}}{{{\mathit{\boldsymbol{N}}_{\rm{f}}}}} $

(19)

根据线性累计疲劳损伤准则，将D中各元素相加可求得管道的累计疲劳损伤D，当D大于1时，即可认为发生了疲劳损伤，故疲劳寿命L可定义为

$ L=\frac{{\tau}}{D} $

(20)

式中，τ为总的计算时长，τ=5 s。

管道的疲劳寿命由最先发生疲劳位置的寿命决定，通过对比各工况下所有应力监测点的寿命，找出最先发生疲劳的位置，以判断不同动量比工况下管道的疲劳寿命。各工况管道的疲劳寿命如表 4所示。

在偏转流工况case1与case2中，管道的预测疲劳寿命分别为1.47×10⁹年和5.40×10¹⁵年，在正常使用时间内不会出现热疲劳现象，而在壁面流工况case3与case4中，管道的预测寿命分别为8.81×10^-4年及4.62×10^-2年。由此可知，降低主支管流体动量比，使搅混流型从壁面流转变为偏转流能有效改善管道的热疲劳现象，但在搅混流型相同时，通过降低动量比的方式并不能改善管道的热疲劳。

4 结论

(1) 在偏转流工况条件下，T型管内流体温度的波动区域相比于壁面流工况条件更接近主管中心，通过降低动量比，使T型管内流体搅混流型从壁面流发展为偏转流，能有效降低管道近壁面流体的温度波动程度。

(2) 相较于壁面流工况，采用偏转流工况能通过减小近壁面流体的温度波动程度来降低管壁热应力的波动程度，从而改善T型管道的热疲劳现象，提高管道的热疲劳寿命。

(3) 只有当动量比降低至T型管内流体的搅混流型发生变化，即由壁面流型发展为偏转流型后，管道的热疲劳现象才能得到改善，流型未发生改变时，调整动量比并不能改善管道的热疲劳现象。