CN113192566A - 反应堆严重事故下熔池瞬态相变模拟方法 - Google Patents
反应堆严重事故下熔池瞬态相变模拟方法 Download PDFInfo
- Publication number
- CN113192566A CN113192566A CN202110484522.3A CN202110484522A CN113192566A CN 113192566 A CN113192566 A CN 113192566A CN 202110484522 A CN202110484522 A CN 202110484522A CN 113192566 A CN113192566 A CN 113192566A
- Authority
- CN
- China
- Prior art keywords
- equation
- iron
- melt
- concentration
- molten pool
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 74
- 230000008859 change Effects 0.000 title claims abstract description 39
- 230000001052 transient effect Effects 0.000 title claims abstract description 35
- 238000004088 simulation Methods 0.000 title claims description 13
- 230000008569 process Effects 0.000 claims abstract description 27
- 239000000126 substance Substances 0.000 claims abstract description 23
- 239000006185 dispersion Substances 0.000 claims abstract description 15
- 239000012530 fluid Substances 0.000 claims abstract description 12
- 238000009499 grossing Methods 0.000 claims abstract description 10
- 238000001228 spectrum Methods 0.000 claims abstract description 9
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 claims description 114
- 229910052742 iron Inorganic materials 0.000 claims description 44
- 239000000155 melt Substances 0.000 claims description 35
- 230000003595 spectral effect Effects 0.000 claims description 32
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 239000011159 matrix material Substances 0.000 claims description 13
- 229910052770 Uranium Inorganic materials 0.000 claims description 12
- QCWXUUIWCKQGHC-UHFFFAOYSA-N Zirconium Chemical compound [Zr] QCWXUUIWCKQGHC-UHFFFAOYSA-N 0.000 claims description 12
- QVGXLLKOCUKJST-UHFFFAOYSA-N atomic oxygen Chemical compound [O] QVGXLLKOCUKJST-UHFFFAOYSA-N 0.000 claims description 12
- 229910052760 oxygen Inorganic materials 0.000 claims description 12
- 239000001301 oxygen Substances 0.000 claims description 12
- JFALSRSLKYAFGM-UHFFFAOYSA-N uranium(0) Chemical compound [U] JFALSRSLKYAFGM-UHFFFAOYSA-N 0.000 claims description 12
- 229910052726 zirconium Inorganic materials 0.000 claims description 12
- 229910052751 metal Inorganic materials 0.000 claims description 7
- 239000002184 metal Substances 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000012886 linear function Methods 0.000 claims description 6
- 238000012885 constant function Methods 0.000 claims description 5
- 229910002059 quaternary alloy Inorganic materials 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 4
- 238000011161 development Methods 0.000 claims description 3
- 238000009792 diffusion process Methods 0.000 claims description 3
- 230000010355 oscillation Effects 0.000 claims description 3
- 101100173586 Schizosaccharomyces pombe (strain 972 / ATCC 24843) fft2 gene Proteins 0.000 claims description 2
- 239000007789 gas Substances 0.000 claims description 2
- 230000005484 gravity Effects 0.000 claims description 2
- 230000005501 phase interface Effects 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims 1
- 238000011160 research Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005290 field theory Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 230000003321 amplification Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000000116 mitigating effect Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000005191 phase separation Methods 0.000 description 1
- 238000004886 process control Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/10—Analysis or design of chemical reactions, syntheses or processes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E30/00—Energy generation of nuclear origin
- Y02E30/30—Nuclear fission reactors
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computing Systems (AREA)
- General Physics & Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Health & Medical Sciences (AREA)
- Mathematical Optimization (AREA)
- Health & Medical Sciences (AREA)
- Chemical Kinetics & Catalysis (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Algebra (AREA)
- Fluid Mechanics (AREA)
- Analytical Chemistry (AREA)
- Mathematical Analysis (AREA)
- Crystallography & Structural Chemistry (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Feedback Control In General (AREA)
Abstract
反应堆严重事故下熔池瞬态相变模拟方法,步骤如下：采用伪二元化简化的熔融物系统，利用耦合了流体力学方程纳维‑斯托克斯方程和相场方程卡恩‑希利亚德方程的方程组作为分层熔池瞬态演变过程中瑞利‑泰勒不稳定现象的控制方程；推导计算求解卡恩‑希利亚德方程所需要的自由能关系，化学势关系，化学迁移率和梯度自由能系数；利用涡量‑流函数方法对纳维‑斯托克斯‑卡恩‑希利亚德方程组进行数值离散，利用伪谱方法求解方程组，求解之前对方程组中所有的非光滑项光滑化；利用熔池瞬态相变模拟软件执行数值运算，模拟出完整的熔池瞬态相变过程并可视化。
Description
技术领域
本发明属于反应堆严重事故现象计算领域，具体涉及反应堆严重事故下熔池瞬态相变模拟方法。
背景技术
严重事故下，堆芯熔融，熔融物穿过下栅格板在堆芯下腔室重定位，出现反应堆熔池。熔融物包容策略被认为是一项有效的事故缓解策略，在该策略下，熔池衰变热的转移通过冷却外壁面实现。熔池主要由铀，氧，锆，铁四种元素组成，由于熔池元素间存在相变间隙，熔池被认为会出现相分离行为，并在热化学相互作用下出现双层或者三层的构型，由于实验条件的严苛性，多层熔池的瞬态演化数值模拟对IVR的安全性评估至关重要。
最近几年，国内在严重事故研究方面开始投入较多的人力和物力，并取得了很多研究成果，但是与核工业发达国家相比，我国对严重事故的研究相对不足，尤其是严重事故下熔池的瞬态相变模拟仍然是一片空白。目前我国尚没有具有自主知识产权的反应堆严重事故下熔池瞬态相变模拟软件。而开发反应堆严重事故下熔池瞬态相变模拟软件对我国核电安全分析有着重要的意义，它能通过对熔池内在相变机理的探究对IVR的安全性分析提供重要参考。
发明内容
为解决上述问题，本发明的目的在于提供一种反应堆严重事故下熔池瞬态相变模拟方法，本发明的原理基于材料学的相场模型，计算流体力学基本方程和熔体元素的热力学相关理论，结合伪谱方法这一高精度的数值解法，利用已开发的熔池瞬态相变模拟软件MPTPC_I执行数值运算，完成了反应堆严重事故下熔池瞬态相变模拟。
为了达到上述目的，本发明采用如下技术方案：
反应堆严重事故下熔池瞬态相变模拟方法，包括如下步骤：
步骤一：熔池瞬态相变过程控制方程开发：采用伪二元化简化的熔融物系统，利用铁的元素浓度指示铀、氧、锆、铁四元系统的相场分布并利用耦合了流体力学方程纳维-斯托克斯方程和相场方程卡恩-希利亚德方程的方程组作为分层熔池瞬态演变过程中瑞利-泰勒不稳定现象的控制方程，其中纳维-斯托克斯方程基于布辛涅司克近似；
步骤二：计算纳维-斯托克斯-卡恩-希利亚德方程组求解所需的热工水力参数。基于四元系统的热力学数据，计算化学迁移率M和梯度自由能系数κ；通过伪二元化简化的熔融物系统中铁浓度c与自由能f的关系f(c)计算自由能f，通过双势阱型的化学势关系μ(c)计算化学势μ；对分段线性的伪二元化简化的熔融物系统的密度和浓度关系p(c)光滑化处理；
步骤三：利用涡量-流函数方法对纳维-斯托克斯-卡恩-希利亚德方程组进行数值离散与求解，纳维-斯托克斯-卡恩-希利亚德方程组采用半隐式离散格式，时间离散采用一阶离散；对离散方程求解时，线性项直接投影到傅里叶空间计算，非线性项则先在物理空间中计算后利用快速傅里叶变换算法投影到傅里叶空间；每个时间步在傅里叶空间中计算出流体的涡量和铁浓度的谱系数；利用快速傅里叶逆变换算法将傅里叶空间中的变量的谱系数重新投影回物理空间；对每一个时间步计算得到的c绘图使得熔池的相变可视化；
步骤四：利用熔池瞬态相变模拟软件执行数值运算，模拟出完整的熔池瞬态相变过程并可视化：首先在每个时间步将铁浓度和涡量投影到傅里叶空间，在傅里叶空间中求解下一个时间步的铁浓度和涡量的谱系数，将计算出的铁浓度和涡量的谱系数从傅里叶空间中投影回物理空间得到此时铁的浓度和涡量，并绘制铁浓度的空间分布指示熔池的相变过程。
本发明与现有技术对比，具有如下优点：
1、采用了相场方法中的卡恩-希利亚德方程描述熔池。采用相场方法的最重要的优势在于数值模拟的过程中不需要对微结构的演化进行界面追踪；
2、采用伪二元化简化的熔融物系统。当前相场模型多着眼于双组分模型，多组分尤其是四元以上的相场模型并不成熟，本发明采用伪二元化简化的熔融物系统使得复杂熔池相变的模拟难度被大大简化，并有更加充分的理论支撑；
3、将材料学的相场理论与流体力学基本方程结合起来，并应用于目前国内研究刚刚起步的核反应堆严重事故领域尤其是在国际上仍是关注重点的分层熔池瞬态演变领域，从材料学机理的角度为严重事故安全性的评价提供新的思路。
4、采用的伪谱法这一数值解法有着很高的精度。伪谱法是传统有限差分法的极限情况，它的核心优势是运算指数收敛，这使得它十分适合于反应堆熔池这种条件极端复杂,模拟精度要求很高的多元系统的瞬态演变模拟。
5、采用涡量-流函数方法求解纳维-斯托克斯方程。涡量-流函数方法是伪谱法解流体力学方程的首选方法，因为它的求解过程中只涉及涡量的运算，且涡量，流函数，速度分量之间的相互转化在伪谱法的支撑下是十分简单快捷的。
6、采用了一种相对简单的光滑化方法对密度函数进行光滑化，有效地解决了吉布斯振荡问题。
附图说明
图1为相场理论的扩散界面示意图。
图2为吉布斯现象示意图。
图3a和图3b分别为光滑化处理前后的函数示意图。
图4为非线性项的谱方法运算示意图。
具体实施方式
步骤一：熔池瞬态相变过程控制方程开发
步骤二：计算纳维-斯托克斯-卡恩-希利亚德方程组求解所需的热工水力参数
步骤三：利用涡量-流函数方法对纳维-斯托克斯-卡恩-希利亚德方程组进行数值离散与求解
步骤四：利用已开发的熔池瞬态相变模拟软件执行数值运算，模拟出完整的熔池瞬态相变过程并可视化
下面进行详细说明：
采用伪二元化简化的熔融物系统，其中一元是铁，另一元是铀、氧、锆，将这三种元素等效成一种元素，利用铁元素的浓度指示由铀，氧，锆，铁四种元素组成的熔融物的相场分布，如图1所示，在这个伪二元化简化的熔融物系统中可以采用c和1-c分别表示两种组成的浓度，此时只需要一个变量即铁元素的浓度c就可以指示熔融物的相场分布，利用耦合了流体力学方程纳维-斯托克斯方程和相场方程卡恩-希利亚德方程的方程组作为分层熔池瞬态演变过程中瑞利-泰勒不稳定现象的控制方程，其中纳维-斯托克斯方程基于布辛涅司克近似，具体的控制方程如下：
其中
式1，2，3，4共同组成了纳维-斯托克斯-卡恩-希利亚德方程组，式中，c为铁元素的摩尔数占比，
为熔融物的速度，M为化学迁移率，f是自由能。式中μ为熔融物的化学势ρ为熔融物的实际密度，p为流体的压力，ρ1为布辛涅司克近似下熔融物的背景密度，η为熔融物的粘度，
表示重力，κ是梯度自由能系数。
步骤二：计算纳维-斯托克斯-卡恩-希利亚德方程组求解所需的热工水力参数。基于四元系统的热力学数据，计算化学迁移率M和梯度自由能系数κ；通过伪二元化简化的熔融物系统中铁浓度c与自由能f的关系f(c)计算自由能f，通过双势阱型的化学势关系μ(c)计算化学势μ；对分段线性的伪二元化简化的熔融物系统的密度和浓度关系ρ(c)光滑化处理；
化学迁移率M由以下公式计算：
式中wint为不同相之间的界面厚度。由于实际的相界面的厚度是介观的尺度，所以在实际模拟中它的值是一个由用户指定的较小的数，这里取L/16，L＝0.176m，本方法把[0,L]×[0,4L]的矩形空间作为计算空间，可以把这种处理看作是保持界面张力不变的前提下对界面进行了一定的放大，这种处理对宏观尺度的数值模拟产生的影响可忽略，wref为实际界面厚度，Vm是摩尔体积，R为摩尔理想气体常数，T为开尔文温度，DU，DZr，DO，DFe分别为铀、锆、氧、铁这四种元素的自扩散系数；
梯度自由能系数κ由以下公式计算：
式中γ为界面张力，cox和cmet分别为相平衡时富氧化物相和富金属相中对应的铁元素浓度。
自由能f由以下公式计算：
f＝A(c-cox)2(c-cmte)2 (8)
其中
化学势μ由以下公式计算：
ρ(c)通常由熔融物的热力学数据决定，通常它是过以下四个点的分段线性函数：
P1(c,ρ)＝(0,ρ_ox_int)
P2(c,ρ)＝(cox,ρ_ox_eq)
P3(c,ρ)＝(cmet,ρ_met_eq)
P4(c,ρ)＝(1,ρ_met_int)
其中ρ_ox_int和ρ_met_int分别为初始富氧化物相和初始富金属相的密度，ρ_ox_eq和ρ_met_eq分别为相平衡时的富氧化物相和富金属相的密度。
下面介绍光滑化的方法：
如图3a所示，由于分段线性函数不是光滑函数，在谱方法中这种非光滑性会导致吉布斯现象，它主要体现在间断处会出现吉布斯振荡，如图2所示，因此可以通过想办法把非光滑函数光滑化来解决这个问题。由于分段线性函数它的导函数是分段常数函数，因此采用tanh()型函数将间断的分段常数函数连接，使得其变成一个连续函数，对该连续函数的积分就是所需的光滑化后的连续函数,如图3b所示。
步骤三：利用涡量-流函数方法对纳维-斯托克斯-卡恩-希利亚德方程组进行数值离散与求解。纳维-斯托克斯-卡恩-希利亚德方程组采用半隐式离散格式，时间离散采用一阶离散；
对卡恩-希利亚德方程进行如下数值离散：
式中：cn和cn+1分别为第n和第n+1时间步长下的铁元素的浓度，τ为时间步长的倒数，
为第n时间步长下的自由能对铁浓度的导数，即化学势；
对纳维-斯托克斯方程进行如下数值离散：
式中：
和
分别表示第n和第n+1时间步长下的熔融物的速度；
对离散方程求解时，线性项直接投影到傅里叶空间计算，非线性项则先在物理空间中计算后利用快速傅里叶变换算法投影到傅里叶空间；线性项指的是变量及它的各阶导数的线性组合，变量在物理空间和傅里叶空间中的相互投影由快速傅里叶变换和快速傅里叶逆变换实现，图4是一个利用伪谱方法进行非线性变量求导的例子，它直观地描述了变量在不同空间中的投影过程。在傅里叶空间中，式(11)被转换成：
其中，
k2＝-(kx)2-(ky)2 (14)
式中：
和
分别表示第n和第n+1时间步长下的铁元素的浓度在傅里叶空间中的谱系数，
表示第n时间步长下的自由能对铁浓度的导数在傅里叶空间中的谱系数，k2对应拉普拉斯算子Δ在傅里叶空间中的傅里叶系数，kx和ky表示x方向和y方向的傅里叶系数，
Nx为x方向的网格数，Ny为y方向的网格数，Nx和Ny都是偶数；
每个时间步在傅里叶空间中计算出流体的涡量和铁浓度的谱系数；基于涡量-流函数基本关系，将式(12)转换为以涡量为变量的方程：
式中：wn和wn+1分别表示第n和第n+1时间步长下的熔融物的涡量，ρn表示第n时间步长下的熔融物的实际密度，在傅里叶空间中，式(15)被转换成：
其中
式中：
和
分别表示第n和第n+1时间步长下的熔融物的涡量在傅里叶空间中的谱系数，cx为铁元素的浓度在x方向的偏导数，cy为铁元素的浓度在y方向的偏导数，wx为熔融物的涡量在x方向的偏导数，wy为熔融物的涡量在y方向的偏导数，cxxy表示铁元素的浓度先对x方向求二阶偏微分后对y方向求一阶偏微分，cyyx表示铁元素的浓度先对y方向求二阶偏微分后对x方向求一阶偏微分，cxxx表示铁元素的浓度对x方向求三阶偏微分，cyyy表示铁元素的浓度对y方向求三阶偏微分；若已知每个时间步下熔融物的涡量和铁元素的浓度的谱系数，求解式(13)和式(16)即可得到下一个时间步下傅里叶空间中对应的值。
利用快速傅里叶逆变换算法将计算获得的谱系数投影回物理空间；
物理空间中的变量和傅里叶空间中变量的谱系数之间的相互投影的方法如下：
物理空间中的变量通过快速傅里叶变换算法即fft算法投影到傅里叶空间中得到变量的谱系数，傅里叶空间中的变量的谱系数通过快速傅里叶逆变换算法即ifft算法重新投影回物理空间得到变量本身，fft算法和ifft算法已经被集成到了MATLAB中并作为函数供用户调用，本发明是在MATLAB平台上完成数值运算的，如果变量c作为一个矩阵，其行序号和列序号分别对应铁浓度分布的x和y方向的空间网格序号，矩阵上的元素表示该位置的铁的浓度，那么fft2(c)表示对该变量矩阵执行二维的快速傅里叶变换，得到的结果为其在傅里叶空间的投影，其结果矩阵的行数和列数与kx和ky的向量元素个数相对应，矩阵上的值为对应kx和ky的变量的谱系数，同理ifft2
表示对变量c的谱系数矩阵执行二维的快速傅里叶逆变换。
绘制每一个时间步的铁的浓度的空间分布使得熔池的相变可视化；
上述公式的推导基于fft，ifft和涡量-流函数的基本原理，即：
Δψ＝-w (26)
u＝ψy (27)
v＝-ψx (28)
式中：ψ表示熔融物的流函数，
表示铁元素的浓度在傅里叶空间中的谱系数；
步骤四：利用已开发的熔池瞬态相变模拟软件执行数值运算，模拟出完整的熔池瞬态相变过程并可视化。首先在每个时间步将铁浓度和涡量投影到傅里叶空间，在傅里叶空间中求解下一个时间步的铁浓度和涡量的谱系数，将计算出的铁浓度和涡量的谱系数从傅里叶空间中投影回物理空间得到此时铁的浓度和涡量，并绘制铁浓度的空间分布指示熔池的相变过程。
Claims (4)
1.反应堆严重事故下熔池瞬态相变模拟方法，其特征在于：包括如下步骤：
步骤一：熔池瞬态相变过程控制方程开发：采用伪二元化简化的熔融物系统，利用铁元素的浓度指示由铀、氧、锆、铁四种元素组成的熔融物的相场分布并利用耦合了流体力学方程即纳维-斯托克斯方程和相场方程即卡恩-希利亚德方程的方程组即纳维-斯托克斯-卡恩-希利亚德方程组作为分层熔池瞬态演变过程中瑞利-泰勒不稳定现象的控制方程，其中纳维-斯托克斯方程基于布辛涅司克近似；
步骤二：计算纳维-斯托克斯-卡恩-希利亚德方程组求解所需的热工水力参数：利用铀、氧、锆、铁四元系统的热力学数据推导计算求解卡恩-希利亚德方程所需要的自由能、化学势、化学迁移率和梯度自由能系数；
步骤三：利用涡量-流函数方法对纳维-斯托克斯-卡恩-希利亚德方程组进行数值离散与求解：利用伪谱方法将方程组中的梯度运算转换为傅里叶空间中的线性运算，并完成方程组的求解，求解之前对方程组中所有的非光滑项光滑化；
步骤四：利用熔池瞬态相变模拟软件执行数值运算，模拟出完整的熔池瞬态相变过程并可视化：具体过程如下：
一、首先在每个时间步将铁浓度和涡量投影到傅里叶空间；
二、在傅里叶空间中求解下一个时间步的铁浓度和涡量的谱系数。
三、将计算出的铁浓度和涡量的谱系数从傅里叶空间中投影回物理空间得到此时铁的浓度和涡量，并绘制铁浓度的空间分布指示熔池的相变过程。
2.根据权利要求1所述的反应堆严重事故下熔池瞬态相变模拟方法，其特征在于：步骤一中所述的熔池瞬态相变过程控制方程开发的过程如下：
一、采用伪二元化简化的熔融物系统；
伪二元化简化的熔融物系统指的是假设在熔融物相变的过程中氧化物相向金属相传递的铀、锆、氧元素始终保持一定的比例，在该假设下，原本含有铀、氧、锆、铁四种元素的熔融物被简化为是由铁元素和另一种混合元素这两种组分组成，这种混合元素对应原本的铀、氧、锆三种元素，此时铁的浓度能指示熔融物的相变情况，因此采用的是伪二元化简化的熔融物系统，其中一元是铁，另一元是铀、氧、锆，将这三种元素等效成一种元素；
二、利用基于布辛涅司克近似的纳维-斯托克斯-卡恩-希利亚德方程组描述熔池瞬态相变过程中发生的瑞利-泰勒不稳定现象；
具体的控制方程如下：
其中
式1，2，3，4共同组成了纳维-斯托克斯-卡恩-希利亚德方程组，式中：c为铁元素的摩尔数占比，
为熔融物的速度，M为化学迁移率，f是自由能,μ为熔融物的化学势，ρ为熔融物的实际密度，p为流体的压力，ρ1为布辛涅司克近似下熔融物的背景密度，η为熔融物的粘度，
表示重力,κ是梯度自由能系数。
3.根据权利要求1所述的反应堆严重事故下熔池瞬态相变模拟方法，其特征在于：步骤二中所述的计算纳维-斯托克斯-卡恩-希利亚德方程组求解所需的热工水力参数所述的过程如下：
一、基于四元系统的热力学数据，计算化学迁移率M和梯度自由能系数κ；
化学迁移率M由以下公式计算：
式中：wint为相变过程中不同相之间的界面厚度，由于实际的相界面的厚度处于介观的尺度，所以在实际的宏观模拟中wint的值是一个由用户指定的较小的数，该较小的数取L/16，L＝0.176m，把[0,L]×[0,4L]的矩形空间作为计算空间，综上把这种处理看作是在保持界面张力不变的前提下对界面进行了一定的放大，这种处理对宏观尺度的数值模拟产生的影响忽略，wref为实际界面厚度,Vm是摩尔体积，R为摩尔理想气体常数，T为开尔文温度，DU，DZr，DO，DFe分别为铀、锆、氧、铁这四种元素的自扩散系数；
梯度自由能系数κ由以下公式计算：
式中：γ为界面张力，cox和cmet分别为相平衡时富氧化物相和富金属相中对应的铁元素浓度；
二、通过伪二元化简化的熔融物系统中铁浓度c与自由能f的关系f(c)计算自由能f，通过双势阱型的化学势关系μ(c)计算化学势μ；
自由能f由以下公式计算：
f＝A(c-cox)2(c-cmte)2 (8)
其中
化学势μ由以下公式计算：
三、对分段线性的伪二元化简化的熔融物系统的密度和浓度关系ρ(c)光滑化处理；
ρ(c)通常由熔融物的热力学数据决定，通常它是过以下四个点的分段线性函数：
P1(c,ρ)＝(0,ρ_ox_int)
P2(c,ρ)＝(cox,ρ_ox_eq)
P3(c,ρ)＝(cmet,ρ_met_eq)
P4(c,ρ)＝(1,ρ_met_int)
其中ρ_ox_int和ρ_met_int分别为初始富氧化物相和初始富金属相的密度，ρ_ox_eq和ρ_met_eq分别为相平衡时富氧化物相和富金属相的密度；
光滑化处理的方法如下：
由于分段线性函数不是光滑函数，在谱方法中这种非光滑性会导致吉布斯现象，主要体现在间断处会出现吉布斯振荡，因此通过把非光滑函数光滑化来解决这个问题；由于分段线性函数的导函数是分段常数函数，因此采用tanh()型函数将间断的分段常数函数连接，使得其变成一个连续函数，对该连续函数的积分就是所需的光滑化后的连续函数。
4.根据权利要求1所述的反应堆严重事故下熔池瞬态相变模拟方法，其特征在于：步骤三中所述的利用涡量-流函数方法对纳维-斯托克斯-卡恩-希利亚德方程组进行数值离散与求解的过程如下：
一、纳维-斯托克斯-卡恩-希利亚德方程组采用半隐式离散格式，时间离散采用一阶离散；
对卡恩-希利亚德方程进行如下数值离散：
式中：cn和cn+1分别为第n和第n+1时间步长下的铁元素的浓度，τ为时间步长的倒数，
为第n时间步长下的自由能对铁浓度的导数，即化学势；
对纳维-斯托克斯方程进行如下数值离散：
式中：
和
分别表示第n和第n+1时间步长下的熔融物的速度；
对离散方程求解时，线性项直接投影到傅里叶空间计算，非线性项则先在物理空间中计算后利用快速傅里叶变换算法投影到傅里叶空间；
线性项指的是变量及其各阶导数的线性组合，物理空间中的变量和傅里叶空间中变量的谱系数的相互投影由快速傅里叶变换和快速傅里叶逆变换实现，在傅里叶空间中，式(11)被转换成：
其中，
k2＝-(kx)2-(ky)2 (14)
式中：
和
分别表示第n和第n+1时间步长下的铁元素的浓度在傅里叶空间中的谱系数，
表示第n时间步长下的自由能对铁浓度的导数在傅里叶空间中的谱系数，k2对应拉普拉斯算子Δ在傅里叶空间中的傅里叶系数，kx和ky表示x方向和y方向的傅里叶系数，
Nx为x方向的网格数，Ny为y方向的网格数，Nx和Ny都是偶数；
二、每个时间步在傅里叶空间中计算出流体的涡量和铁浓度的谱系数；
基于涡量-流函数基本关系，将式(12)转换为以涡量为变量的方程：
式中：wn和wn+1分别表示第n和第n+1时间步长下的熔融物的涡量，ρn表示第n时间步长下的熔融物的实际密度，在傅里叶空间中，式(15)被转换成：
其中
式中：
和
分别表示第n和第n+1时间步长下的熔融物的涡量在傅里叶空间中的谱系数，cx为铁元素的浓度在x方向的偏导数，cy为铁元素的浓度在y方向的偏导数，wx为熔融物的涡量在x方向的偏导数，wy为熔融物的涡量在y方向的偏导数，cxxy表示铁元素的浓度先对x方向求二阶偏微分后对y方向求一阶偏微分，cyyx表示铁元素的浓度先对y方向求二阶偏微分后对x方向求一阶偏微分，cxxx表示铁元素的浓度对x方向求三阶偏微分，cyyy表示铁元素的浓度对y方向求三阶偏微分；若已知每个时间步下熔融物的涡量和铁元素的浓度的谱系数，求解式(13)和式(16)即能够得到下一个时间步下傅里叶空间中对应的值；
三、利用快速傅里叶逆变换算法将计算获得的谱系数投影回物理空间；
物理空间中的变量和傅里叶空间中变量的谱系数之间的相互投影的方法如下：
物理空间中的变量通过快速傅里叶变换算法即fft算法投影到傅里叶空间中得到变量的谱系数，傅里叶空间中的变量的谱系数通过快速傅里叶逆变换算法即ifft算法重新投影回物理空间得到变量本身，fft算法和ifft算法已经被集成到了MATLAB中并作为函数供用户调用，如果变量c作为一个矩阵，其行序号和列序号分别对应铁浓度分布的x和y方向的空间网格序号，矩阵上的元素表示该位置的铁的浓度，那么fft2(c)表示对该变量矩阵执行二维的快速傅里叶变换，得到的结果为其在傅里叶空间的投影，其结果矩阵的行数和列数与kx和ky的向量元素个数相对应，矩阵上的值为对应kx和ky的变量的谱系数，同理
表示对变量c的谱系数矩阵执行二维的快速傅里叶逆变换。
四、绘制每一个时间步的铁的浓度的空间分布使得熔池的相变可视化。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110484522.3A CN113192566B (zh) | 2021-04-30 | 2021-04-30 | 反应堆严重事故下熔池瞬态相变模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110484522.3A CN113192566B (zh) | 2021-04-30 | 2021-04-30 | 反应堆严重事故下熔池瞬态相变模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113192566A true CN113192566A (zh) | 2021-07-30 |
| CN113192566B CN113192566B (zh) | 2022-12-09 |
Family
ID=76983757
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110484522.3A Active CN113192566B (zh) | 2021-04-30 | 2021-04-30 | 反应堆严重事故下熔池瞬态相变模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113192566B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113946954A (zh) * | 2021-10-14 | 2022-01-18 | 中国核动力研究设计院 | 核反应堆压力容器下腔室熔融池瞬态结构获取方法及装置 |
| CN115064222A (zh) * | 2022-06-17 | 2022-09-16 | 上海玫克生储能科技有限公司 | 一种锂电池工况预测方法和系统 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407623A (zh) * | 2016-11-15 | 2017-02-15 | 南京航空航天大学 | 基于相场法的瞬态下焊接过程微观组织演变模拟方法 |
| CN110457791A (zh) * | 2019-07-26 | 2019-11-15 | 嘉兴学院 | 基于特征线的龙格-库塔时间离散的流体力学有限元算法 |
| CN112102894A (zh) * | 2020-09-04 | 2020-12-18 | 西安交通大学 | 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法 |
-
2021
- 2021-04-30 CN CN202110484522.3A patent/CN113192566B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407623A (zh) * | 2016-11-15 | 2017-02-15 | 南京航空航天大学 | 基于相场法的瞬态下焊接过程微观组织演变模拟方法 |
| CN110457791A (zh) * | 2019-07-26 | 2019-11-15 | 嘉兴学院 | 基于特征线的龙格-库塔时间离散的流体力学有限元算法 |
| CN112102894A (zh) * | 2020-09-04 | 2020-12-18 | 西安交通大学 | 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法 |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113946954A (zh) * | 2021-10-14 | 2022-01-18 | 中国核动力研究设计院 | 核反应堆压力容器下腔室熔融池瞬态结构获取方法及装置 |
| CN115064222A (zh) * | 2022-06-17 | 2022-09-16 | 上海玫克生储能科技有限公司 | 一种锂电池工况预测方法和系统 |
| CN115064222B (zh) * | 2022-06-17 | 2024-07-12 | 上海玫克生储能科技有限公司 | 一种锂电池工况预测方法和系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113192566B (zh) | 2022-12-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113192566B (zh) | 反应堆严重事故下熔池瞬态相变模拟方法 | |
| Roelofs et al. | Status and perspective of turbulence heat transfer modelling for the industrial application of liquid metal flows | |
| Dinh et al. | Turbulence modelling for large volumetrically heated liquid pools | |
| Su et al. | Numerical analysis of steam condensation over a vertical surface in presence of air | |
| Xing et al. | Closure to “discussion of ‘factors of safety for Richardson Extrapolation’”(2011, ASME J. Fluids Eng., 133, p. 115501) | |
| Yadigaroglu | Computational Fluid Dynamics for nuclear applications: from CFD to multi-scale CMFD | |
| Lakehal et al. | Direct numerical simulation of turbulent heat transfer across a mobile, sheared gas-liquid interface | |
| Agrawal et al. | Hydrogen distribution in nuclear reactor containment during accidents and associated heat and mass transfer issues—a review | |
| CN105247623B (zh) | 核反应堆中模拟流体流动及计算燃料组件机械形变的方法 | |
| Grahn et al. | Coupling of the 3D neutron kinetic core model DYN3D with the CFD software ANSYS-CFX | |
| CN106844853B (zh) | 结合阻力和能量分布包含格架搅混效应的子通道分析方法 | |
| Goluskin | Internally heated convection beneath a poor conductor | |
| Wang et al. | Validation of a methodology for thermal stratification analysis in sodium‐cooled fast reactors | |
| CN114091310B (zh) | 反应堆严重事故中包壳行为多尺度多物理场耦合分析方法 | |
| Tsige-Tamirat et al. | A review of models for the sodium boiling phenomena in sodium-cooled fast reactor subassemblies | |
| Borgohain et al. | Natural circulation studies in a LBE loop for a wide range of temperature | |
| Mazères et al. | Contribution to modeling of hydrogen effect on oxygen diffusion in Zy-4 alloy during high temperature steam oxidation | |
| Zou et al. | Validation of Pronghorn with the SANA Experiments | |
| Yadigaroglu | Thermal-hydraulic issues and achievements–A historical perspective | |
| Valette | Analysis of subchannel and rod bundle PSBT experiments with CATHARE 3 | |
| Taylor et al. | Mesh-independent large-eddy simulation with anisotropic adaptive mesh refinement for hydrogen deflagration prediction in closed vessels | |
| Mergui et al. | Transient double diffusive convection in a vertical enclosure with asymmetrical boundary conditions | |
| Zhuang et al. | Nonuniform three-layer models to predict transient flows in buoyancy-driven natural ventilation with a localized heat source | |
| CN105247622A (zh) | 用于模拟核反应堆中容器内流体流动的方法以及用于计算核反应堆堆芯组件的机械形变的方法，以及相关的计算机程序产品 | |
| Kim et al. | Analysis for Factors Affecting Molten Corium Concrete Ablation and Coolability |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |