[发明专利]核反应堆事故工况下耦合传热传质的包壳氧化分析方法

权利要求书

1.一种核反应堆事故工况下耦合传热传质的包壳氧化分析方法，其特征在于：包括如下步骤：

步骤1：建立包壳的几何结构：环形包壳由三个相构成，从内侧到外侧分别为β-Zr，α-Zr(O)和ZrO₂；包壳外侧与水或水蒸汽接触，包壳内侧与包壳和芯块间的气隙接触；氧原子由外侧进入包壳；

步骤2：建立计算模型的控制方程，具体包括如下内容：

1)氧原子扩散方程：

式中：F_O为氧原子通量密度，v为变形速度，C_O为氧原子浓度，D为氧原子扩散系数，r为空间坐标，t为时刻；

2)锆原子守恒方程：

式中：F_Zr为锆原子通量密度，C_Zr为锆原子浓度；

3)导热方程：

式中：Q为热流密度，ρ为密度，C_p为比热容，T为温度，λ为导热系数；

步骤3：建立相界面迁移方程，具体包括如下内容：

1)在β-Zr与α-Zr(O)的界面和α-Zr(O)与ZrO₂的界面，分别建立如下方程：

式中：p和p-1分别为相邻两相的编号，C_O,p为相界面处p相的氧原子浓度，C_O,p-1为相界面处p-1相的氧原子浓度，ξ为相界面位置，v_p为相界面处p相的变形速度，v_p-1为相界面处p-1相的变形速度，D_p为相界面处p相的氧原子扩散系数，D_p-1为相界面处p-1相的氧原子扩散系数，C_Zr,p为相界面处p相的锆原子浓度，C_Zr,p-1为相界面处p-1相的锆原子浓度；

2)在包壳和芯块间的气隙与β-Zr的界面和ZrO₂与水蒸汽的界面，分别建立如下方程：

式中：ξ₁为包壳和芯块间的气隙与β-Zr的界面位置，ξ₄为ZrO₂与水蒸汽的界面位置，v₁为包壳和芯块间的气隙与β-Zr界面处的变形速度，v₃为ZrO₂与水蒸汽界面处的变形速度；

步骤4：设置相界面处的边界条件，具体包括如下内容：

1)建立氧原子扩散方程的边界条件：

对于β-Zr与α-Zr(O)的界面和α-Zr(O)与ZrO₂的界面，在恒温条件下建立的边界条件为：式中：为相界面处p-1相的平衡氧原子浓度，为相界面处p相的平衡氧原子浓度；对于温度瞬变条件，边界条件修正为：式中：τ为弛豫时间；

对于包壳和芯块间的气隙与β-Zr的界面，建立的边界条件为：

对于ZrO₂与水蒸汽的界面，若蒸汽不足，建立的边界条件为：式中：F_O,0为已知的氧原子通量密度，随时间变化；若蒸汽充足，在恒温条件下建立的边界条件为：式中：为ZrO₂与水蒸汽界面处ZrO₂的平衡氧原子浓度；对于温度瞬变条件，边界条件修正为：

2)建立导热方程的边界条件：

对于β-Zr与α-Zr(O)的界面和α-Zr(O)与ZrO₂的界面，建立的边界条件为：T_p-1＝T_p，式中：λ_p-1为相界面处p-1相的导热系数，λ_p为相界面处p相的导热系数，T_p-1为相界面处p-1相的温度，T_p为相界面处p相的温度，为相界面处p-1相的温度梯度，为相界面处p相的温度梯度，Q_ch为相界面处的化学反应热流密度；

对于包壳和芯块间的气隙与β-Zr的界面和ZrO₂与水蒸汽的界面，根据不同的模拟条件有三种可能的边界条件，分别为：①第一类边界条件：T＝T₀，式中：T₀为已知的温度，随时间变化；②第二类边界条件：式中：Q₀为已知的热流密度，随时间变化；③第三类边界条件：式中：h_G为对流换热系数，T_G为包壳外部流体主流温度；

步骤5：离散步骤2建立的控制方程，基于步骤4建立的边界条件求解获得氧原子浓度场、锆原子浓度场、温度场、变形速度场；具体包括以下子步骤：

1)对初始时刻的包壳几何结构划分网格；

2)采用控制容积积分法分别离散氧原子扩散方程、锆原子守恒方程和导热方程，得到：

式中：a_O,P1、a_O,W、a_O,E、a_O,P2、a_Zr,w、a_Zr,e、a_Zr,b、a_T,P1、a_T,W、a_T,E、a_T,P2为离散方程中的系数，分别为t+Δt时刻离散节点P、W、E处的氧原子浓度，为t时刻离散节点P处的氧原子浓度，分别为t+Δt时刻网格界面w、e处的变形速度，分别为t+Δt时刻离散节点P、W、E处的温度，为t时刻离散节点P处的温度；

3)采用三对角阵算法分别求解导热方程和氧原子扩散方程的离散方程，获得t+Δt时刻的温度场和氧原子浓度场；

4)根据温度场和氧原子浓度场更新平衡氧原子浓度、氧原子扩散系数、导热系数、密度和比热容；

5)氧原子浓度和锆原子浓度的数量关系为：C_OM_O+C_ZrM_Zr＝ρ，式中M_O和M_Zr分别为氧原子和锆原子的摩尔质量；根据氧原子浓度和锆原子浓度的数量关系，计算得到t+Δt时刻的锆原子浓度场；

6)根据锆原子浓度场，计算锆原子守恒方程的离散方程中的系数a_Zr,w、a_Zr,e、a_Zr,b；联立离散方程组，求解t+Δt时刻的变形速度场；

7)若本次迭代计算得到的温度场、氧原子浓度场、变形速度场与上一次迭代计算结果间的相对误差小于设定值，则执行步骤6，否则再次执行子步骤2)～7)；

步骤6：基于步骤3建立的相界面迁移方程，计算各个相界面的迁移速度；采用梯形求积公式计算Δt时间内各个相界面的位移，获得t+Δt时刻各个相界面的位置，更新包壳的几何结构；

步骤7：若此时的计算时间大于设定的终止时间，则结束计算，否则对更新后的包壳几何结构划分网格，执行步骤5～7。