CN114023478B - 一维铅基系统瞬态安全分析方法 - Google Patents

Abstract

一种一维铅基系统瞬态安全分析方法，其主要步骤如下：1、根据输入文件，确定铅基系统的组织结构和组成台架各个组件的组件参数，设定计算时间和时间步长。2、系统初始化计算，得到零时刻的初始值。3、采用吉尔算法，计算当前t时刻的流场、温度场和压力长。4、根据步骤3的计算结果，采用吉尔算法进行下一时刻的计算，直到设定的总计算时间，计算停止。5、若计算达到规定时间，则终止计算。6、根据步骤1～5得到的计算结果，采用准稳态方法计算台架内氧浓度、铁浓度和氧化层浓度的分布。本发明可以分析铅铋快堆或铅基系统的稳态热工水力特性和事故工况下的瞬态相应特性，为推进铅铋快堆的设计和安全特性分析提供建议和指导。

Description

一维铅基系统瞬态安全分析方法

技术领域

本发明设计先进反应堆设计领域，具体涉及到一种铅铋反应堆及铅铋试验系统的瞬态安全分析方法。

背景技术

随着核能技术的进步，第四代反应堆的概念应运而生。铅铋反应堆作为第四代反应堆的一种，其所有具有的自然循环能力、高中子经济型等特点已经引起了学界的广泛重视。在设计和研究铅铋反应堆系统时，应当对铅铋反应堆的瞬态安全性引起足够重视。对于压水堆等传统反应堆，已经出现了一些系统计算程序，如RELAP,TRAC等。但由于铅铋合金的流动传热特性与水等常规流体之间存在着巨大的差异，因此这些程序并不能直接用于铅铋反应堆的建模。

在铅铋反应堆内，有一些独特的热工水力现象，传统的系统计算程序无法捕捉到并解释这些热工水力现象。因此，有必要针对铅铋反应堆内独特的热工水力现象，开发出一款专用于铅铋反应堆即铅基试验系统的瞬态安全分析软件。

发明内容

本发明公开一种一维铅基系统的瞬态安全分析方法。该方法针对铅铋反应堆及铅基试验系统的结构特点和铅铋反应堆内的特殊的热工水力现象，进行建模和计算分析。通过瞬态计算，分析反应堆在一些典型的事故工况下的响应特性，以此为依据评估安全特性，为铅笔反应堆的设计和安全特性分析提供建议和指导。

为了实现上述目的，本发明采取了以下技术方案予以实施：

一种一维铅基系统瞬态安全分析方法，步骤如下：

步骤1：根据输入文件，确定铅基系统的结构。

步骤2：计算初始值。

步骤3：计算瞬态值。

步骤4：利用吉尔算法进行计算。

步骤5：输出瞬态结果。

步骤6：根据计算得到的热工水力分析结果，进行氧控参数的初步计算.

步骤7：利用准稳态方法，计算氧控分析数据。

步骤8：输出氧控分析结果。

和现有的技术相比，本发明具有以下优点：

1.实现了铅基系统的模块化设计和模拟计算；

2.实现了铅基系统的一维氧化腐蚀分析，能够针对热工水力瞬态计算得到的流场计算结果计算得到氧化层的变化趋势。

3.与传统的系统分析程序相比，程序专为铅铋反应堆和铅基实验系统而设计，能够模拟在铅铋反应堆中特有的热工水力现象，计算可靠性更高。

附图说明

图1为本发明的计算方法N-S图。

图2为划分的管道堆芯组件的控制体。

图3a为堆芯组件的截面示意图，图3b为堆芯控制体的简图。

图4为铅池模型示意图。

图5为划分的换热器控制体。

图6为氧控装置的示意图。

具体实施方式

下面结合附图和具体实施方式对本发明方法作进一步详细说明：

如图1所示，本发明一种空间核电推进系统球床反应堆物理热工耦合分析方法，步骤如下：

步骤1：用户根据所模拟的铅铋实验台架的组织结构和组件参数，模拟计算的总时间和时间步长，模拟计算的关键热工参数和瞬态情况，编写并配置输入文件。根据输入文件的信息，铅铋实验台架被视为数个组件的有机结合，将每个组件按照需求沿径向划分为若干个控制体。

步骤2：根据铅铋实验台架的构成结构和初始条件，进行铅铋实验台架的稳态初始化计算，从而得到全台架的每个控制体的冷却剂流量和温度，作为零时刻的初始值。

堆芯的热功率：

P_t＝P_n (1)

第i个控制体的冷却剂质量流量：

M_i＝M_in (2)

第i个控制体的冷却剂压力：

P_i＝P_i+1+ΔP_fi+ΔP_ai (3)

第i个控制体的冷却剂温度：

式(1)-(4)中：

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_p,g——冷却剂比容/J·kg^-1·K^-1

ΔV_i——第i个控制体的体积/m³

q”’_c——堆芯的功率密度/W·m^-3

步骤3：建立铅铋实验台架的堆芯中子物理模型、冷却剂的热工水力模型和铅铋冷却剂的流动传热模型和热工水力模型，其中铅铋冷却剂的热工水力模型包括换热器、堆芯、管道、铅池的热工水力模型。各种组件的模型如下:

堆芯中子物理模型，采用点堆动力学方程求解堆芯的裂变功率，点堆动力学方程考虑六组缓发中子和反应性反馈；

式(5)-(6)中：

P(t)——t时刻堆芯裂变功率/W

t——时间/s

ρ(t)——t时刻总反应性/$

β_eff——缓发中子的总份额

Λ——中子代时间/s

λ_i——第i组缓发中子的衰变常数/s^-1

C_i(t)——t时刻第i组缓发中子先驱核的浓度/m^-3

β_i——第i组缓发中子的份额

停堆后，反应堆的功率由反应堆余热功率提供：

式(7)中：

P_decay——t时刻堆芯衰变功率

P₀——停堆前堆芯裂变功率

A、τ、α——常系数

用户也根据裂变产物的衰变模型进行衰变余热的计算：

式(8)中：

P_decay——t时刻堆芯衰变功率

N——衰变产物种类数量

——第i种裂变产物的半衰期/s

u_i——第i种裂变产物的衰变功率/W

——第i种裂变产物的份额

n_fts——堆芯裂变功率/W

冷却剂的热工水力模型：认为冷却剂在堆芯内的流动是不可压缩的径向一维流动。控制方程如下：

式(9)-(11)中：

ρ——流体密度/kg·m^-3

z——控制体高度/m

W——流体流量/kg·s^-1

A——控制体流道面积/m²

P——压强/Pa

f——摩擦阻力系数

D_e——等效水力直径/m

H——流体焓/J·kg^-1

U——流体速度/m·s^-1

管道模型：在计算过程中，假设流体是不可压的，因此对于整根管道而言，流体的流量是不变的。管道模型的结构如图2所示，在计算过程中，管道组件被划分成了数个控制体，并包括一个壁面边界条件(图中以阴影部分表示)，图中T_w,i是第i个节点的壁面温度，W_i是第i个控制体内的流量，H_i是第i个节点的流体焓值，W_input,W_output分别是管道组件的入口流量和出口流量。

管道温度的变化率：

管道流量的变化率：

式(12),(13)中：

ρ——流体密度/kg·m^-3

V——控制体体积/m³

H_i——第i个控制体的焓/J·kg^-1

W——流体流量/kg·s^-1

Q——热边界处的热量/W

L_ctrl——控制体长度/m

A——控制体流通面积/m²

P_in——管道入口压强/Pa

P_out——管道出口压强/Pa

ΔP_i——各个控制体的压力差/Pa

式中，热边界中指定的热量可以是由用户指定的定热流边界条件，或者由点堆模型计算得出的反应堆功率。堆芯的节点模型如图3a和图3b所示分别为堆芯组件的截面示意图和堆芯控制体的简图。

堆芯模型：

式(14),(15)中：

ρ——流体密度/kg·m^-3

V——控制体体积/m³

C_p,i——第i个控制体内流体的比热容/J·kg^-1

T_i——第i个控制体的温度/℃

n——燃料棒数量

D_cs——燃料棒包壳直径/m

l_i——第i个控制体的长度/m

h——换热系数/J·K^-1

ρ_w——包壳密度/kg·m^-3

C_p,w——包壳比热容/J·kg^-1

V_i,w——第i个控制体的包壳体积/m³

T_w,i——第i个控制体的包壳温度/℃

Q_t,i——第i个控制体的燃料功率/W

外部热量Q_t可设置为恒定热流密度或是由点堆方程计算导出。

在堆芯通道内的流动传热关系式和摩擦阻力关系式：

铅池模型：铅池模型由一个入口、出口和一个控制体组成，如图4所示，图4显示了铅池模型的简图。铅池内待求解的量包括了铅池的液位高度L_p和铅池内液体体积m_p，它们的变化取决于入口管道和出口管道的流量W_in,W_out，同时也考虑了因为铅池内流体膨胀导致的膨胀流量W_ex。

式(19)-(20)中：

L_p——液位高度/m

m_p——铅池内质量/kg

ν——铅池内液体比体积/m³·kg^-1

T_pool——铅池内液体温度/℃

W_in——入口流量/kg·s^-1

W_out——出口流量/kg·s^-1

W_ex——膨胀流量/kg·s^-1

膨胀流量的计算公式如下：

式(21)中：

V_i——铅池第i个控制体的体积/m³

ρ_i——铅池第i个控制体的密度/kg·m^-3

换热器模型：换热器模型如图5所示。一二次侧之间的换热量计算如下：

Q＝kAΔT (22)

其中k是导热系数，根据换热器的结构，可以算出来一次侧和二次侧的导热系数：

式(23)-(24)中：

k_p——一次侧导热系数/W·m^-2·K^-1

k_s——二次侧导热系数/W·m^-2·K^-1

h_p——一次侧换热系数/W·m^-2·K^-1

h_s——二次侧换热系数/W·m^-2·K^-1

R_pp——一次侧管壁热阻/m²·K·W^-1

R_ps——二次侧管壁热阻/m²·K·W^-1

R_mp——一次侧夹层热阻/m²·K·W^-1

R_ms——二次侧夹层热阻/m²·K·W^-1

为了方便计算，将中间夹层分开，并分别计算了各个部分的热阻：

式(25)-(29)中：

κ_tube——管壁材料导热系数W·m^-1·K^-1

D_ip——一次侧管壁内径/m

D_op——一次侧管壁外径/m

D_is——一次侧夹层内径/m

D_os——一次侧夹层外径/m

D_av——平均直径/m

那么一次侧和二次侧之间传输的热量为：

式(30)中：

Q_p——根据一次侧计算的传热量/W

Q_s——根据二次侧计算的传热量/W

k_p——一次侧传热系数/W·m^-2·K^-1

k_s——二次侧传热系数/W·m^-2·K^-1

A_p——一次侧传热面积/m²

A_s——二次侧传热面积/m²

T_p——一次侧温度/℃

T_w——壁面温度/℃

T_s——二次侧温度壁温的计算公式可以由能量守恒导出：

式(31)中：

ρ_wi——第i个控制体的管壁密度/kg·m^-3

V_wi——第i个控制体的管壁体积/m³

C_pwi——第i个控制体的管壁比热容/J·kg^-1·K^-1

T_wi——第i个控制体的管壁温度/℃

Q_pi——第i个控制体的一次侧热量/W

Q_si——第i个控制体的二次侧热量/W

步骤4：根据步骤3的计算结果，利用吉尔算法进行下一时刻的计算，直到达到设定的总计算时间，热工水力计算停止。

步骤5：将步骤1～4得到的计算值输出到输出文件中。该输出文件可以用各类数据分析程序进行分析，如origin,matlab等。

步骤6：将步骤1～4得到的计算值进行进一步整理，得到当前时刻的流场信息，根据该信息进行氧化腐蚀行为的计算。首先计算初始时刻的流体内氧浓度，铁浓度和氧化层浓度。

氧浓度控制方程：

铁浓度控制方程：

式中(32)-(33)：

C_O——氧元素浓度/wt％

C_Fe——铁元素浓度/wt％

j_O——氧的输运通量/m·wt％·s^-1

——氧浓度源/wt％·s^-1

j_Fe——铁通量/m·wt％·s^-1

U——管道截面周长/m

A——管道截面积/m²

u——流体速度/m·s^-1

L_ctrl——控制体长度/m

氧的输运通量：

式中(34)-(36)：

——氧化层溶解速率

j_ox——氧化层生长速率/m·wt％·s^-1

——Fe溶解至LBE后的耗氧速率j_Fe——铁的溶解速率/m·wt％·s^-1

M_O——氧的分子量

M_Fe——铁的分子量

α——铁在氧化层中的比率

M_Cr——铬的分子量

最低氧浓度限值：

铁元素的扩散模型：

C_Fe,s＝C_O ^-1.3310^{11.35-12844/T} (38)

氧化模型：

式(37)-(39)中：

C_O,min——最小氧浓度/wt％

T_f——流体温度/℃

C_Fe,s——LBE内铁的最大浓度/wt％

δ——氧化层厚度/cm

K_P1——氧化层增长速率/cm²·s^-1

K_P2——氧化层溶解速率/cm²·s^-1

316L不锈钢在纯气体环境下的氧化层增长速率：

T91不锈钢在纯气体环境下的氧化层增长速率：

式(40)-(41)中：

——氧在LBE中的分压

气体氧分压和液态LBE中氧浓度的转化关系式：

氧化层消耗速率：

氧控装置如图6所示，氧元素从壁面(图中A_cs处)输运而来，并进入管道，而管道中流体图中A为氧控装置横截面积，A_cs为氧控装置的表面积。氧控装置的供氧模型：

壁面厚度变化：

式(42)-(45)中：

ρ_LBE——LBE密度/kg·m^-3

ρ_ox——氧化层密度/kg·m^-3

f_ox ^Fe——铁元素在氧化层中所占的质量分数

C^out——氧控装置出口氧浓度/wt％

Cⁱⁿ——氧控装置入口氧浓度/wt％

L_ctrl——控制体长度/m

K——质量输运系数

δ——氧化层长度/m

C^s——氧化层表面氧浓度/wt％

A_cs——氧控装置表面积/m²

K_p——氧化层总生长速率/m·wt％·s^-1

步骤7：将步骤6得到初始时刻的量按照准稳态方法进行计算，得到每一时刻的流体内的氧浓度、铁浓度和氧化层浓度。

步骤8：将步骤7得到的计算值输出到输出文件当中。该输出文件可以用各类数据分析程序进行分析，如origin,matlab等。

Claims (1)

1.一种一维铅基系统瞬态安全分析方法，其特征在于：步骤如下：

步骤1：用户根据所模拟的铅铋实验台架的组织结构和组件参数，模拟计算的总时间和时间步长，模拟计算的关键热工参数和瞬态情况，编写并配置输入文件；根据输入文件的信息，铅铋实验台架被视为数个组件的有机结合，将每个组件按照需求沿径向划分为若干个控制体；

步骤2：根据铅铋实验台架的构成结构和初始条件，进行铅铋实验台架的稳态初始化计算，从而得到全台架的每个控制体的冷却剂流量和温度，作为零时刻的初始值；

堆芯的热功率：

P_t＝P_n (1)

第i个控制体的冷却剂质量流量：

M_i＝M_in (2)

第i个控制体的冷却剂压力：

P_i＝P_i+1+ΔP_fi+ΔP_ai (3)

第i个控制体的冷却剂温度：

式(1)-(4)中：

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_p,g——冷却剂比容/J·kg^-1·K^-1

ΔV_i——第i个控制体的体积/m³

q″′_c——堆芯的功率密度/W·m^-3

步骤3：建立铅铋实验台架的堆芯中子物理模型、冷却剂的热工水力模型和铅铋冷却剂的流动传热模型，其中铅铋冷却剂的流动传热模型包括换热器、堆芯、管道、铅池的热工水力模型；各种组件的模型如下:

堆芯中子物理模型，采用点堆动力学方程求解堆芯的裂变功率，点堆动力学方程考虑六组缓发中子和反应性反馈；

式(5)、(6)中：

P(t)——t时刻堆芯裂变功率/W

t——时间/s

ρ(t)——t时刻总反应性/$

β_eff——缓发中子的总份额

Λ——中子代时间/s

λ_i——第i组缓发中子的衰变常数/s^-1

C_i(t)——t时刻第i组缓发中子先驱核的浓度/m^-3

β_i——第i组缓发中子的份额

停堆后，反应堆的功率由反应堆余热功率提供：

式(7)中：

P_decay——t时刻堆芯衰变功率

P₀——停堆前堆芯裂变功率

A、τ、α——常系数

用户也根据裂变产物的衰变模型进行衰变余热的计算：

式(8)中：

P_decay——t时刻堆芯衰变功率/W

N——衰变产物种类数量

——第i种裂变产物的半衰期/s

u_i——第i种裂变产物的衰变功率/W

——第i种裂变产物的份额

n_fts——堆芯裂变功率/W

冷却剂的热工水力模型：认为冷却剂在堆芯内的流动是不可压缩的径向一维流动，控制方程如下：

式(9)-(11)中：

ρ——流体密度/kg·m^-3

z——控制体高度/m

W——流体流量/kg·s^-1

A——控制体流道面积/m²

P——压强/Pa

f——摩擦阻力系数

D_e——等效水力直径/m

H——流体焓/J·kg^-1

U——流体速度/m·s^-1

管道模型：在计算过程中，假设流体是不可压的，因此对于整根管道而言，流体的流量是不变的；

管道温度的变化率：

管道流量的变化率：

式(12),(13)中：

ρ——流体密度/kg·m^-3

V——控制体体积/m³

H_i——第i个控制体的焓/J·kg^-1

W——流体流量/kg·s^-1

Q——热边界处的热量/W

L_ctrl——控制体长度/m

A——控制体流通面积/m²

P_in——管道入口压强/Pa

P_out——管道出口压强/Pa

ΔP_i——各个控制体的压力差/Pa

式中，热边界中指定的热量是由用户指定的定热流边界条件，或者由点堆模型计算得出的反应堆功率；

堆芯模型：

式(14),(15)中：

ρ——流体密度/kg·m^-3

V——控制体体积/m³

C_p,i——第i个控制体内流体的比热容/J·kg^-1

T_i——第i个控制体的温度/℃

n——燃料棒数量

D_cs——燃料棒包壳直径/m

l_i——第i个控制体的长度/m

h——换热系数/J·K^-1

ρ_w——包壳密度/kg·m^-3

C_p,w——包壳比热容/J·kg^-1

V_i,w——第i个控制体的包壳体积/m³

T_w,i——第i个控制体的包壳温度/℃

Q_t,i——第i个控制体的燃料功率/W

外部热量Q_t设置为恒定热流密度或是由点堆方程计算导出；

在堆芯通道内的流动传热关系式和摩擦阻力关系式：

铅池模型：铅池模型由一个入口、出口和一个控制体组成；铅池内待求解的量包括了铅池的液位高度L_p和铅池内液体体积m_p，它们的变化取决于入口管道和出口管道的流量W_in,W_out，同时也考虑了因为铅池内流体膨胀导致的膨胀流量W_ex；

式(19)、(20)中：

L_p——液位高度/m

m_p——铅池内质量/kg

ν——铅池内液体比体积/m³·kg^-1

T_pool——铅池内液体温度/℃

W_in——入口流量/kg·s^-1

W_out——出口流量/kg·s^-1

W_ex——膨胀流量/kg·s^-1

膨胀流量的计算公式如下：

式(21)中：

V_i——铅池第i个控制体的体积/m³

ρ_i——铅池第i个控制体的密度/kg·m^-3

换热器模型：一二次侧之间的换热量计算如下：

Q＝kAΔT (22)

其中k是导热系数，根据换热器的结构，计算出一次侧和二次侧的导热系数：

式(23)-(24)中：

k_p——一次侧导热系数/W·m^-2·K^-1

k_s——二次侧导热系数/W·m^-2·K^-1

h_p——一次侧换热系数/W·m^-2·K^-1

h_s——二次侧换热系数/W·m^-2·K^-1

R_pp——一次侧管壁热阻/m²·K·W^-1

R_ps——二次侧管壁热阻/m²·K·W^-1

R_mp——一次侧夹层热阻/m²·K·W^-1

R_ms——二次侧夹层热阻/m²·K·W^-1

为了方便计算，将中间夹层分开，并分别计算各个部分的热阻：

式(25)-(29)中：

κ_tube——管壁材料导热系数W·m^-1·K^-1

D_ip——一次侧管壁内径/m

D_op——一次侧管壁外径/m

D_is——一次侧夹层内径/m

D_os——一次侧夹层外径/m

D_av——平均直径/m

那么一次侧和二次侧之间传输的热量为：

式(30)中：

Q_p——根据一次侧计算的传热量/W

Q_s——根据二次侧计算的传热量/W

k_p——一次侧传热系数/W·m^-2·K^-1

k_s——二次侧传热系数/W·m^-2·K^-1

A_p——一次侧传热面积/m²

A_s——二次侧传热面积/m²

T_p——一次侧温度/℃

T_w——壁面温度/℃

T_s——二次侧温度

壁温的计算公式由能量守恒导出：

式(31)中：

ρ_wi——第i个控制体的管壁密度/kg·m^-3

V_wi——第i个控制体的管壁体积/m³

C_pwi——第i个控制体的管壁比热容/J·kg^-1·K^-1

T_wi——第i个控制体的管壁温度/℃

Q_pi——第i个控制体的一次侧热量/W

Q_si——第i个控制体的二次侧热量/W

步骤4：根据步骤3的计算结果，利用吉尔算法进行下一时刻的计算，直到达到设定的总计算时间，热工水力计算停止；

步骤5：将步骤1～4得到的计算值输出到输出文件中；该输出文件用各类数据分析程序进行分析；

步骤6：将步骤1～4得到的计算值进行进一步整理，得到当前时刻的流场信息，根据该信息进行氧化腐蚀行为的计算；首先计算初始时刻的流体内氧浓度，铁浓度和氧化层浓度；

氧浓度控制方程：

铁浓度控制方程：

式(32)-(33)中：

C_O——氧元素浓度/wt％

C_Fe——铁元素浓度/wt％

j_O——氧的输运通量/m·wt％·s^-1

——氧浓度源/wt％·s^-1

j_Fe——铁通量/m·wt％·s^-1

U——管道截面周长/m

A——管道截面积/m²

u——流体速度/m·s-1

L_ctrl——控制体长度/m

氧的输运通量：

式(34)-(36)中：

——氧化层溶解速率

j_ox——氧化层生长速率/m·wt％·s^-1

——Fe溶解至LBE后的耗氧速率

j_Fe——铁的溶解速率/m·wt％·s^-1

M_O——氧的分子量

M_Fe——铁的分子量

α——铁在氧化层中的比率

M_Cr——铬的分子量

最低氧浓度限值：

铁元素的扩散模型：

氧化模型：

式(37)-(39)中：

C_O,min——最小氧浓度/wt％

T_f——流体温度/℃

C_Fe,s——LBE内铁的最大浓度/wt％

δ——氧化层厚度/cm

K_P1——氧化层增长速率/cm²·s^-1

K_P2——氧化层溶解速率/cm²·s^-1

316L不锈钢在纯气体环境下的氧化层增长速率：

T91不锈钢在纯气体环境下的氧化层增长速率：

式(40)-(41)中：

——氧在LBE中的分压

气体氧分压和液态LBE中氧浓度的转化关系式：

氧化层消耗速率：

氧控装置的供氧模型：

壁面厚度变化：

式(42)-(45)中：

ρ_LBE——LBE密度/kg·m^-3

ρ_ox——氧化层密度/kg·m^-3

f_ox ^Fe——铁元素在氧化层中所占的质量分数

C^out——氧控装置出口氧浓度/wt％

Cⁱⁿ——氧控装置入口氧浓度/wt％

L_ctrl——控制体长度/m

K——质量输运系数

δ——氧化层长度/m

C^s——氧化层表面氧浓度/wt％

A_cs——氧控装置表面积/m²

K_p——氧化层总生长速

步骤7：将步骤6得到初始时刻的量按照准稳态方法进行计算，得到每一时刻的流体内的氧浓度、铁浓度和氧化层浓度；

步骤8：将步骤7得到的计算值输出到输出文件中，该输出文件能够用各类数据分析程序进行分析。