矩形管道载氚热磁流体输运特性的数值模拟研究.pdf

June2023NuclearFusionand Plasma Physics2023年6 月Vol.43,No.2第43卷第2 期子体物理核聚变与等离文章编号：0 2 54-6 0 8 6(2 0 2 3)0 2-0 16 1-0 8D0I:10.16568/j.0254-6086.202302007矩形管道载氙热磁流体输运特性的数值模拟研究韩佳佳，汪卫华，浦文婧1，江海燕3，储德林，史博1(1.陆军炮兵防空兵学院，合肥2 30 0 31；2.安徽大学，合肥2 30 0 39；3.合肥工业大学，合肥2 30 0 0 9)摘要：聚变堆液态金属包层矩形管道中的氙输运过程与磁流体动力学(MHD)流动传热过程耦合在一起，形成了复杂的载氙热磁流体输运特性。基于开发的MHD流动与传热数值模拟程序对矩形管道中液态金属MHD流动传热特性及其氙输运的影响进行了数值模拟。该程序首先求解了动量守恒方程，并与理论解进行了对比验证，然后与能量守恒方程耦合求解，得到了温度影响下矩形管道中的液态金属流场分布，在此基础上对强磁场高核热梯度影响下的氙浓度分布进行了数值模拟，得到了氙浓度在管道中的分布特性。结果显示，液态金属在矩形管道中的流动传热对氙输运过程产生了显著影响。关键词：聚变堆；MHD；氙输运；数值方法中图分类号：TL334文献标识码：A1引言氙是聚变堆的重要燃料，聚变反应中所需的氙主要来源于高能中子与包层中的氙增殖剂进行的反应。按照氙增殖剂形态的不同，可以把包层分为固态包层和液态包层。液态包层主要以液态金属锂或者含锂的液态金属合金(如PbLi等)作为氙增殖材料。氙在液态金属中的输运问题不仅关系到聚变燃料的来源，而且与聚变堆稳定运行的安全性密切相关。通过研究液态金属中的氙输运特性，不仅可以为液态包层氙提取系统优化设计提供参数支持，也能够为氙输运过程中的安全性分析提供参考。在聚变堆高温强磁场环境中，流动的液态金属在强磁场作用下产生了感应电流，磁场与感应电流的相互作用产生了阻碍流体运动的洛伦兹力，形成了复杂的磁流体动力学(MHD)现象2 ,MHD现象改变了液态包层中的流场分布，同时也影响到了氙的输运特性。谢波等人对锂铅合金中的释氙行为进行了理论分析4 与实验研究5。曾正魁等人利用SuperMC软件系统分析了DFLL包层中的氙增殖性能，并进行了实验验证6 。贺青云对液态包层复杂几何构件下的MHD流动、传热和氙输运问题进行了研究，分析了氙溶解度与扩散系数等关键变量对氙分布的影响7 。由于矩形管道是液态金属包层中的一种典型结构8 ，王俊等人利用COMSOL软件研究了强磁场下方管中锂铅流体内的氙输运及其分布9。目前的研究主要集中于物性参数不随温度发生变化的理想情况下的氙输运特性。但在聚变堆环境下，高能中子进入包层管道时能够产生高功率非均匀分布的核热，形成了包层管道中的非均匀温度场，管道中液态金属的物性参数也随之发生变化，影响到了流场分布3，改变了包层管道中不同位置处的氙浓度，从而影响了氙的输运特性。本文基于自主开发的聚变堆液态金属包层收稿日期：2 0 2 1-0 2-2 6；修订日期：2 0 2 2-11-0 4基金项目：国家重点研发计划(2 0 18 YFE0310400)；国家自然科学基金(1197 50 2 2)；安徽省自然科学基金(2 0 0 8 0 8 5QA37)作者简介：韩佳佳(198 5-)，男，山东东营人，博士，从事核聚变包层技术研究。162第43 卷核聚变与等离子体物理MHD流动和传热数值模拟程序MHPl101，考虑温度对粘性系数的影响，将动量方程与能量方程耦合求解，得到矩形管道中的速度分布与温度分布，并根据质量传输方程求解得到氙在矩形管道中的输运特性，可为液态金属包层的优化设计提供参考。2计算模型2.1控制方程本文对不可压缩液态金属在矩形管道中的流动传热与氙输运特性进行了数值模拟。不考虑液态金属重力的影响，对于完全发展的磁流体管流，用于计算速度和温度分布的无量纲方程为：Vu=0(1)Qu+u.Vu=-Vp+=V(Vu)+N(JB)(2)atReV20=V.(uxB)(3)aT1+V(uT)V?T+EcNJ?+Q(4)tPe其中，u、p 和T分别为无量纲化后的液态金属速度、压力和温度；J、B、和分别为无量纲化后的感应电流、磁感应强度、电势和能量源项；Re、N、Pe 和Ec分别表征液态金属流动传热特性的雷诺数、相互作用参数、贝克莱数和埃克特数。无量纲参数与液态金属的物性参数有关，而物性参数会随着温度的改变而发生变化。以液态金属锂铅的粘性系数为例,其动力粘性系数与温度的关系为 :11640u=1.8710-exp(5)RT根据式(5)，在聚变堆包层工况条件下，如果液态金属的温度从50 0 K变化到6 0 0 K，动力粘性系数将下降37.3%。与粘性系数不变时的理想状态相比，粘性系数的变化改变了管道中的流场分布和温度分布，影响了液态金属中氙的输运特性。液态金属中的氙由高能中子和锂反应生成，其输运行为主要包括氙的生成、对流和扩散。反应生成的氙离子在很短的范围内停止运动并在液态金属增殖剂中传输，因此氙的输运过程主要受到对流和扩散的影响，描述两者作用强弱之比的无量纲数施密特数(Sc）定义为：Sc=L(6)pD式中，p为液态金属密度；D为氙在液态金属中的扩散系数，由于实验装置和曲线拟合模型的差异，不同实验中测得的扩散系数随温度的变化关系不同12 ,本文的计算模型采用的是Fauvet的测量结果13。此时 PbLi 的施密特数为194，对流作用远大于扩散作用，可以根据质量传输方程得到氙在矩形管道中的浓度分布从质量守恒定律出发，建立以下的质量传输方程：2Cr+V(ucr)=DVc+ST(7)t式中，C为液态金属中的氙浓度；S为中子在液态金属中的产氙率。在研究MHD流动传热特性对氙输运影响时，将质量传输方程进行无量纲化，质量传输方程的无量纲形式为：0CL+V(ucr)=1V2CT+ST(8)atPeL式中，PeL为质量传输方程无量纲贝克莱数，其定义为：Pe,=UoL(9)D电磁学与流体力学基本方程(1)(4)及质量传输方程(8)共同组成了矩形管道中求解氙输运特性的控制方程组。2.2定解条件为求得控制方程的数值解，需要给定由初始条件和边界条件组成的定解条件。2.2.1速度定解条件矩形管道入口处，速度边界条件设置为：u=1,V=W=0(10)式中，u、V、w 分别为速度u沿x、y、z 轴的分量。对于充分发展流动，速度在壁面处应用无滑移边界条件：=V=W=O(11)管道出口处为Neumann边界条件：dudvdw=0(12)dxdxdx(14)163韩佳佳等：矩形管道载氙运特性的数值模拟研究第2 期2.2.2压力定解条件压力在管道入口与壁面处设置为Neumann边界条件：0(13)dn管道出口处的边界条件设置为：p=02.2.3电势定解条件控制方程组中电磁学方程中的电势在入口、出口与壁面处的边界条件设置为Neumann边界条件：=0(15)an2.2.4温度定解条件由于聚变堆包层被连接固定于真空室中，真空室处于低温恒温器的内部，是密封的环形不锈钢容器15。真空室与包层之间通过辐射传递的热量远小于聚变反应产生的中子在包层中沉积的核热，因此在数值模拟中没有考虑包层与真空室之间的热量交换，管道壁面上设置为绝热边界条件：dTdT=0(16)dydz根据液态包层中常用的液态金属锂铅的实际工况，矩形管道入口温度值设置为：Tim.=753K(17)出口温度为充分发展Neumann边界条件：aT=0(18)an2.2.5氙浓度定解条件计算模型中管道入口处的氙浓度设置为：C=0(19)假设氙在锂铅流体与管道壁面之间建立了化学平衡，那么壁面上的氙浓度与平衡分压系数相关9,即:C,=KCL(20)其中，C,与CL分别为壁面两侧锂铅流体和结构材料中的氙浓度；K为壁面上的平衡分压系数16 ，其值为8.510 4。氙浓度的出口边界条件设置为Neumann边界条件：aCL=0(21)dn3数值方法矩形管道中液态金属的速度和温度分布采用耦合方法求解。其中，动量守恒方程通过相容守恒格式计算电流和洛伦兹力，基于不可压缩流动投影算法17 数值求解速度分布。能量守恒方程采用有限体积法数值求解3。将动量守恒方程和能量守恒方程耦合求解后可以得到矩形管道中的速度场分布，然后通过求解质量传输方程(8)得到管道中氙浓度分布。数值模拟算法的流程如图1所示。流场初始化计算液态金属物性参数计算NS方程得到速度分布计算能量方程得到温度分布否计算收敛是计算质量传输方程得到氙浓度分布否计算收敛是输出计算结果图1数值模拟程序算法流程图控制方程组采用非均匀结构化网格进行离散，其中速度u、压力p、电势、温度和氙浓度储存在网格中心。为保证电流在计算过程中守恒，将电流密度J储存在网格面元中心，并对边界层网格加密，在垂直于磁场的哈特曼层内划分为七个网格，平行于磁场的侧层内划分为九个网格。在动量守恒方程与能量守恒方程的数值求解的基础上10 ，采用有限体积法对氙质量传输方程进行离散：164第43卷核聚变与等离子体物理acdv+V.（(uc)dv=vc.dv+S.dv(22)离散后的方程采用隐式方法进行求解，其中时间项离散为：.0CAxAyAzT(23)tAt扩散项离散为：v2c.dv=CTi+1.j.k人ij.k-1.jkAyAzXi+1.j.kXi.j.kXi.j,k-Xi-1,j,kCTJkCCj+1,kTi.j-1.kAxAz(24)XXi.j.kj+1,ki,j,kXi,j-1,kCTiJkCj.k+1i,j.ki.j.k-1+AxAyXj.k+1,j，kik式中，对流项采用QUICK格式进行离散，采用延迟修正方法将网格界面上的氙浓度表示为：cH=c+(cH-cl)*(25)式中，cH与c分别为离散项的高阶格式与低阶格式；*表示上一时间步的迭代值，低阶格式采用一阶迎风格式。采用延迟修正后的对流项与源项分别离散为：JV(ucr)dv=(AyAz(26)S.dv=StAxAyAz+S*(27)式中，Sr为产氙率；S*为由于在QUICK格式中采用延迟修正方法而引入的源项因为粘性系数在数值计算的迭代过程中发生变化，动量守恒方程与能量守恒方程的耦合求解过程非常复杂。为避免计算结果的发散，在耦合求解过程中采用了松弛技术。对于离散后的质量传输代数方程，使用交替方向隐式迭代法(ADI方法)，分别沿x轴、y轴和z轴方向采用三对角阵算法(TDMA方法)逐行扫描迭代求解。为了提高迭代收敛速度，在每次迭代完成后使用块修正技术进行加速4结果与讨论为了验证数值模拟程序，建立了矩形管道中充分发展的MHD液态金属流动模型，流动区域的无量纲尺寸设置为6 2 2。使用Fortran90语言编写数值模拟程序代码，在加载Intel Visual FortranCompilerXE2011的VisualStudio2010平台上编译运行。为加快计算过程运行速度，数值模拟程序利用个人计算机多核共享存储结构的特征，采用多线程程序设计技术OpenMP，在个人计算机平台上实现了8 线程并行计算。为分析网格数量与程序计算结果之间的关联性，数值模拟程序分别将流体区域划分为636363(网格1)、6 96 96 9(网格2)和7 57 575(网格3)三个网格进行网格无关性验证。在矩形管道中间位置x=3.0，z=0 处水平方向上的速度与温度计算结果如图2 所示。从图2 中可以看出，不同网格数量情况下计算结果的变化很小，说明该数值模拟程序具有良好的网格无关性。计算结果如图2所示。2网格！a1.8网格网格31.61.41.210.80.60.40.20-0.8 0.6-0.4-0.200.40.6 0.81000网格！950网格900网格38508007507006506005505000.8-0.6-0.4-0.200.40.60.8图2网格无关性验证结果a速度计算结果；b温度计算结果。4.1温度影响下的速度分布在不考虑温度对液态金属物性参数影响时，165韩佳佳等：矩形管道载氙热磁流体输运特性的数值模拟研究第2 期MHD流动计算中的 Shercliff算例与Hunt算例有解析解，通常将其作为MHD数值模拟程序的基准算例2 0 。本文中的数值模拟程序分别计算了哈特曼数Ha=300，R e=10、H a=50 0，R e=10 与 Ha=1000,Re=500三种情况下的Shercliff与Hunt算例中的速度分布情况。图3显示了Ha=500时的速度分布，从图3中可以看出，数值解和解析解非常接近。不同基准算例中的流量与压力梯度的计算结果列于表1中，其中，压力梯度的最大误差为6%，流量的最大误差为0.1%。计算结果表明，本文中的数值模拟程序可以准确地模拟管道中的液态金属MHD流动分布。1.2a10.8一解析解数值解0.60.40.20-1-0.500.5J9b876一解析解5数值解4320-0.500.5图3Shercliff与Hunt算例计算结果a-Shercliff算例(Ha=500,Re=10);b-Hunt算例(Ha=1000,Re=500)。表1不同算例中的压力梯度与流量计算结果Shercliff算例Hunt算例哈特曼数5001000压力梯度解析解5.263.47压力梯度数值解4.983.36压力梯度计算误差5.32%3.17%入口流量4.004.00出口流量数值解3.9983.999流量计算误差0.05%0.025%考虑到温度对液态金属粘性系数的影响，在MHD数值模拟程序的基础上，本文将动量方程与能量方程耦合求解，得到了温度影响下Hunt算例中的速度场分布。图4为温度影响下速度场分布的数值模拟结果，可以看出，受到洛伦兹力的作用，在矩形管道侧壁附近出现了剪切流。剪切流的出现增强了液态金属在侧壁的换热能力18 ，降低了该区域的温度。由式(5)可知，温度的降低增加了侧壁附近液态金属的粘性系数，而高的粘性系数会减慢液态金属的流动速度。而在Hunt算例解析解的计算过程19 中，没有考虑到温度变化的影响，因此与本文中的数值模拟结果存在差异。a10.5u03.50N02.9212.34-0.521.76二11.180.530.6000.540.02110解析解b数值解86421-0.500.5图4温度影响下Hunt算例中的速度场分布a-速度分布数值模拟；b数值解与解析解的对比。Ha=1000,Re=500。4.2矩形管道中的氙浓度分布本文中计算的氙输运模型为：矩形管道尺寸(长第43卷166核聚变与等离子体物理宽高)为0.3m0.1m0.1m，管道壁厚0.0 1m，侧壁与垂直磁场方向的哈特曼壁均为导电壁面，管道内为液态锂铅流体，入口锂铅流速uo=0.01ms-。管道位于磁感应强度B=5T的外部磁场中，磁场方向沿z轴，管道截面如图5所示。0.2m0.01m侧壁N液态锂铅Y轴哈特曼壁图5矩形管道截面图液态锂铅的物性参数取锂铅材料数据库中的数值，管道内部的产氙率ST一般通过中子学计算或实验结果得到，本文中的产氙率计算模型采用数值拟合方法得到的双冷锂铅包层(DCLL)中的产氛率14 为：St=2.0610-+5.7010exp(0.0061(28)+2.2510-exp(-0.08276这里y为管道内部与第一壁(FW)的距离。根据4.1节中液态金属在矩形管道中的流动传热特性，再利用MHD数值模拟程序数值求解质量输运方程计算得到了管道中的氙浓度分布。图6 为管道出口处的液态金属温度分布、速度分布与氙浓度分布的数值模拟结果。abu/10-2ms-T/K1.517901.267841.017780.767720.51766Z0.26Z7600.017542.5Cd一产氛率氙浓度21.5c/10molm.s55.470.541.6127.74_8.05-0.02500.0250.0513.870图6矩形管道出口处速度、温度与氙浓度分布图速度分布；b温度分布；c氙产生率与氙浓度分布的对比；d氙浓度分布。a图6 a为管道出口处的速度分布，图6 b为管道出口处的温度分布。在侧壁附近出现剪切流，增强了换热能力，因此管道侧壁的温度小于管道中心的温度。同时，侧壁附近的剪切流也增强了氙在液态167第2 期韩佳佳等：矩形管道载氙热磁流体输运特性的数值模拟研究金属中的输运作用，使得侧壁附近的氙浓度小于该处的氙产生率，利用均值归一化后管道出口处的氙产生率与氙浓度分布如图6 c所示。此外，管道的形状对液态金属的换热和氙的输运也会产生影响，如图6 b所示，矩形管道弯曲处的几何形状为直角，不利于热量交换，因此这一区域的温度要高于管道中心温度。该区域同样不利于氙的扩散与对流，导致管道壁四个直角弯曲处的氙浓度远高于管道中心。氙在管道出口处的浓度分布如图6 d所示，在模型工况下，数值计算的结果显示，矩形管道出口处的氙浓度最大值为55.47 10 molm3.s-l，是这一区域氙浓度平均值的5.7 3倍。过高的氙浓度会显著增加管道处的氙压，增加了氙泄漏的风险，是液态金属包层设计中需要关注的重要问题5总结基于开发的MHD数值模拟程序，研究了矩形管道中液态金属的流动传热特性及其对氙输运的影响。通过将动量守恒方程与能量守恒方程耦合求解，获得了温度影响下矩形管道中液态金属的速度分布，并采用与标准算例对比的方法验证了数值模拟结果的可靠性。在获得管道中速度分布的基础上，通过数值求解质量传输方程得到了四周均为导电壁的矩形管道中的氙浓度分布。结果显示，氙浓度分布与管道中的流场分布呈现出显著的负相关特性，并在管道壁面相交处产生了氙的累积效应，使得氙浓度显著增大，为开展聚变堆液态金属包层中氙输运问题的研究提供了参考。参考文献：1 吴宜灿，郁杰，胡丽琴，等聚变堆安全特性评价研究.核科学与工程,2 0 16,36(6):8 0 2-8 10.2Morgan L,Pasley J.Tritium breeding control withinliquid metal blankets.Fusion Engineering and Design,2018,33(3):107-112.3 Han J J,Zuo C J,Wang W H.Numerical investigation ofthe MHD flow under the influence of heat transfer inrectangular ducts J.IEEE Transacations on PlasmaScience,2019,47(2):1399-1404.4谢波，翁葵平锂铅合金中氙扩散行为的理论研究.分子科学学报,2 0 0 9,2 5(5):352-356.5 谢波，吴宜灿，陈晓军，等锂铅合金释氙实验研究.原子核物理评论,2 0 11,2 8(3):37 1-37 6.6曾正魁，蒋洁琼，王海霞，等双功能液态锂铅包层氙增殖性能分析.核科学与工程，2 0 2 0，40(2):218-226.7贺青云.液态包层复杂几何构件下的MHD流动、传热和氙输运研究D.合肥：中国科学技术大学,2 0 17.8 Abdou M,Morley N B,Smolentsev S,et al.Blanket/firstwall challenges and required R&D on the pathway toDEMO J.Fusion Engineering and Design,2015,100:2-43.9王俊，张国书，张龙，等.方管锂铅磁流体氙增殖剂中的氙输运模拟.核聚变与等离子体物理，2 0 18,38(3):254-259.10韩佳佳，汪卫华.聚变堆液态金属包层MHD流动和传热数值模拟程序开发与验证J.核科学与技术，2 0 19,7(3):83-90.11 Mas de les Valls E,Sedano L A,Batet L,et al.Lead-lithium eutectic material database for nuclear fusiontechnology J.Journal of Nuclear Materials,2008,376(3):353-357.12 Ricapito I.PbLi/T database:status of the knowledge R.ENEA CR Brasimone,FPN-FISING,2007.13 Fauvet P,Sannier J.Hydrogen behaviour in liquid17Li83Pb alloy J.Journal of Nuclear Materials,1988,155(1):516-519.14 Zhang H J.Modeling and analysis of tritium transport inmulti-region lead-lithium liquid metal blankets D.LosAngeles:University of California,2014.15袁保山，姜韶风，陆志鸿.托卡马克装置工程基础M北京：原子能出版社，2 0 11.38 1-38 7.16 Zhang H,Ying A,Abdou M,et al.Modeling tritiumtransport in PbLi breeder blankets under steady state J.Fusion Science&Technology,2011,60(2):814-818.17 Ni M J,Munipalli R,Morley N B,et al.A current densityconservative scheme for incompressible MHD flows at alow magnetic Reynolds number.Part I:On a rectangularcollocated grid system J.Journal of ComputationalPhysics,2007,277(1):174-204.18周涛聚变堆液态金属磁流体动力学模拟程序研发与应用研究D.北京：中国科学院研究生院,2 0 11.19 Hunt J C R.Magnetohydrodynamic flow in rectangularducts J.Journal of Fluid Mechanics,1965,21(4):577-590.20 Smolentsev S,Badia S,Bhattacharyay R,et al.Anapproach to verification and validation of MHD codes forfusion applications J.Fusion Engineering and Design,2015,100:65-72.168第43卷核聚变与等离子体物理Numerical simulation on transport characteristics of tritium loadedthermomagnetic fluid in rectangular ductHAN Jia-jia,WANG Wei-hua,PU Wen-jing,JIANG Hai-yan3,CHU De-lin,SHI Bo(1.PLA Army Academy of Artillery and Air Defense,Hefei 230031;2.Anhui University,Hefei 230039;3.Hefei University of Technology,Hefei 230009)Abstract:The transport process of tritium in the liquid metal blanket of fusion reactor is coupled with theflow and heat transfer process of MHD,resulting in the complex tritium transport characteristics in the strongmagnetic field.Based on the developed MHD flow and heat transfer numerical simulation program,the flow andheat transfer characteristics of liquid metal MHD in a rectangular duct and its influence on tritium transport werenumerically simulated.Firstly,the momentum conservation equation is solved and validated by the theoreticalsolution,and then coupled with the energy conservation equation,the liquid metal flow field distribution in therectangular duct under the influence of temperature is obtained.On this basis,the tritium concentrationdistribution under the influence of strong magnetic field and high nuclear thermal gradient is numericallysimulated,and the distribution characteristics of tritium concentration in the duct are obtained.The results showsthat the flow and heat transfer of liquid metal in rectangular duct have significant effect on tritium transport.Key words:Fusion reactor;MHD;Tritium transport;Numerical method