详情

原文传递利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法

专利名称：	利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法
摘要：	一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，主要步骤如下：1、针对具体问题进行粒子建模，设定粒子初始布置；2、根据各粒子位置划分网格；3、在时间步长Δt之后，使用有限容积法求解能量方程；4、粒子法处理显式计算动量方程中的粘性项和源项，估算粒子速度及位置；5、计算粒子数密度，求解压力Poisson方程；6、计算速度的修正值，并计算各粒子在下一时层的速度及位置；7、输出计算结果，并根据求解的具体问题输出壁面热流密度和熔融池内物质的相态。本方法具有粒子法精确捕捉自由液面及模拟相变问题的优点，又具有较高的计算精度，能够通过此方法来研究高瑞利数熔融池换热特性，为核电厂反应堆严重事故安全特性研究提供重要依据。
专利类型：	发明专利
国家地区组织代码：	陕西;61
申请人：	西安交通大学
发明人：	田文喜;李勇霖;陈荣华;苏光辉;秋穗正
专利状态：	有效
申请日期：	2019-05-13T00:00:00+0800
发布日期：	2019-07-23T00:00:00+0800
申请号：	CN201910394189.X
公开号：	CN110044959A
代理机构：	西安智大知识产权代理事务所
代理人：	何会侠
分类号：	G01N25/20(2006.01);G;G01;G01N;G01N25
申请人地址：	710049 陕西省西安市碑林区咸宁西路28号
主权项：	1.一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，其特征在于：步骤如下：步骤1：对高瑞利数熔融池进行粒子建模，不同种类的粒子代表不同的物质，每种粒子根据所表示的物质具有其对应的质量、密度和熔化焓，使用粒子i来表示某一粒子，则其质量、密度、熔化焓对应标记分别为mi、ρi、设定粒子的初始布置，包括粒子的位置、速度、压力和温度，对应标记分别为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0；步骤2：根据各粒子的位置划分网格：在本步骤中，设置一个粒子搜索半径，对于每一个粒子i，其坐标位置为(xi，yi，zi)，与其间距lij≤1.2l0的所有粒子作为搜索的对象，式中lij为j粒子与i粒子的距离，j粒子的坐标标记为(xj，yj，zj)，lij使用公式(1)进行计算，并将搜索到的粒子使用j1，j2，……，jn符号进行标记，其中下标n表示搜索到的粒子总数；将j1，j2，……，jn粒子中的每三个粒子进行组合，将有种组合，由公式(2)进行计算，对每个组合中的粒子进行考察，在这些组合中，每个组合内的三个坐标点必须满足构成面的要求，取其中的一个组合进行解析，其中一种组合粒子标记为j1，j2，j3，粒子中心分别标记为O1，O2，O3，空间坐标对应为(x1，y1，z1)，(x2，y2，z2)，(x3，y3，z3)，连接O1、O2为直线O1O2，连接O1、O3为直线O1O3，连接O2、O3为直线O2O3，若要判断O1、O2、O3能否构成一个平面，只需要判断直线O1O2与直线O1O3之间存在一个不为0的夹角，则只需满足公式(3)，对所有的组合进行判断，将能够构成面的组合保留，不能构成面的组合删去，随后考虑i粒子与每个组合中三个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体，因为组合中的三个粒子能够构成平面，所以要使这四个粒子能够构成四面体，只需满足i粒子不在组合中三个粒子所构成的平面上即可，添加i粒子，粒子中心标记为O，空间坐标定义为(xi，yi，zi)，若要使O、O1、O2、O3四个点构成四面体，因为O1、O2、O3三个点已经满足构成面的要求，则只需要满足i粒子到平面O1O2O3的距离不为0，使用公式(4)进行判定，式中： ——i粒子到平面O1O2O3的距离m； a、b、c、d——平面O1O2O3的平面方程的四个系数，通过公式(5)获得， i粒子与j1、j2、j3粒子构成的平面O1O2O3的距离若满足公式(4)，则说明i粒子能够与三角形O1O2O3构成四面体OO1O2O3，四面体OO1O2O3的四个面分别为三角形OO1O2、三角形OO1O3、三角形OO2O3和三角形O1O2O3；使用海伦公式计算四个三角形的面积，四个面的面积对应标记为海伦公式如公式(6)所示，式中， S——三角形的面积； L1、L2、L3——三角形的三边长； p——三角形的半周长，即p＝(L1+L2+L3)/2；四面体OO1O2O3的体积通过公式(7)进行计算，式中， V——四面体体积，在四面体OO1O2O3中，标记为 S——三角形O1O2O3的面积，标记为使用公式(6)计算； h——O点与三角形所在平面O1O2O3的距离，记为使用公式(4)计算；在使用粒子位置划分网格时，还有几个参量需要求解，其一是四面体OO1O2O3中四个顶点O、O1、O2、O3分别对四面体总体积的分属控制体积，对应标记为VO、这个分属控制体积由公式(8)进行计算，实际上，四面体OO1O2O3的四个点的分属控制体积与其相对的三角形面积成反比，即有公式(9)所示的关系，其二是对四面体OO1O2O3中四个面的三角形进行面积划分，在三角形O1O2O3中，三角形O1O2O3的面积通过公式(6)求出；点M12为线段O1O2的中点，点M13为线段O1O3的中点，点M23为线段O2O3的中点。此时需要在三角形中寻找一点P，连接P、O1，连接P、O2，连接P、O3，连接P、M12，连接P、M13，连接P、M23，此时三角形O1O2O3被分成了6个小三角形，分别标记为三角形O1PM12，三角形O1PM13，三角形O2PM12，三角形O2PM23，三角形O3PM13，三角形O3PM23；因为点M13为线段O1O3的中点，所以三角形O1PM13与三角形O3PM13面积相等，标记为同样的，三角形O1PM12与三角形O2PM12面积相等，标记为三角形O2PM23与三角形O3PM23面积相等，标记为对于各个三角形的面积，有如公式(10)和公式(11)所示的几何关系，三角形的三个角的角度由公式(3)所示的余弦定理求出；由公式(10)和公式(11)求出三角形O1O2O3所分成的6个小三角形的面积；同理也能求得四面体OO1O2O3另外三个面所分成的共计18个小三角形的面积；在四面体OO1O2O3中寻找一点Q，对于顶点O的控制体积而言，连接QP1，QP2，QP3，QM1，QM2，QM3，QO，至此点O的控制体积分成了六个四面体，分别为四面体OP1M1Q，四面体OP1M2Q，四面体OP2M1Q，四面体OP2M3Q，四面体OP3M2Q和四面体OP3M3Q，其体积对应标记为和对于四面体的体积，有公式(12)所示的关系，设点Q距离平面O1O2O3、平面OO1O2，平面OO1O3，平面OO2O3的距离分别为H、H1、H2、H3，对于顶点O的控制体积列如公式(13)的方程，同理对于顶点O1、O2、O3也列出相似的方程，联立这4个方程及公式(8)能够求出H、H1、H2、H3的值；至此，在O、O1、O2、O3组成的四面体系统中，O1对O起影响的点O的控制体积为记为O2对O起影响的点O的控制体积为记为O3对O起影响的点O的控制体积为记为所以，在O、O1、O2、O3组成的四面体系统中，点O1对O的有效传热面积标记为点O2对O的有效传热面积标记为点O1对O的有效传热面积标记为通过公式(14)计算，式中， ——O点到O1点的距离； ——O点到O2点的距离； ——O点到O3点的距离；以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传热面积数据之后才能开始进行温度计算；步骤3：在温度计算中，能量守恒方程如公式(15)所示，式中， T——温度K； t——时间s κ——导热系数W/(m2·K)； ρ——密度kg/m3； cp——比定压热容J/(kg·K)； ST——温度的源项；公式(16)给出能量方程的Δt时间步长的积分形式，式中， t——时间/s； Δt——时间步长s； T——温度K； V——体积m3； CV——控制体积m3； κ——导热系数W/(m2·K)； ρ——密度kg/m3； cp——比定压热容J/(kg·K)； ST——温度的源项；对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(16)能够离散成公式(17)的形式，公式(17)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；式中： Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传热面积； κ——导热系数W/(m2·K)； rj——j粒子的坐标矢量； ri——i粒子的坐标矢量； ——j粒子在n时层的温度K； Tin——i粒子在n时层的温度K； ρi——i粒子的密度kg/m3； cpi——i粒子的比定压热容J/(kg·K)； Tin+1——i粒子在n+1时层的温度，即在n时层的时间步长Δt后的温度K；求解公式(17)即得到i粒子在时间步长Δt之后的温度Tin+1；步骤4：显式计算动量方程中的粘性项及源项，得到粒子速度及位置的估算值及ri；动量方程如公式(18)所示式中， u——速度m/s； t——时间s； p——压力Pa v——运动粘度m2/s； β——热膨胀系数m-1； g——重力加速度m2/s； Tm——背景温度，定义为环境温度K； τ——粒子应力项，表示不可解的亚粒子尺度的流动对可解尺度的贡献，这一项是用于求解高瑞利数熔融池内湍流的；显式计算动量方程的粘性项及源项，即计算项，其中粘性项的拉普拉斯算子由公式(19)离散求解，式中， d——计算的维度； n0——初始粒子数密度，通过计算，是核函数，使用公式(20)计算，只需将括号内的自变量修改为即可，为j粒子初始时刻的坐标矢量，为i粒子初始时刻的坐标矢量； λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(21)计算，式中， re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；速度和位置的估算值由公式(22)进行计算，式中， ——i粒子在上一时间步的速度矢量m/s； rin——i粒子在上一时间步的位置矢量m； Tin——i粒子在上一时间步的温度K；步骤5：根据估算得到的粒子位置，通过公式(23)计算此时的粒子数密度ni，只需要将公式中的\|rj-ri\|使用位置的估算值替代即可，由公式(24)所示的压力Poisson方程求解得到此时各粒子的压力值pi；公式(24)等式的左边使用如同公式(19)的形式离散，式中， ——i粒子在n+1时层的压力Pa； α、γ——静态系数；步骤6：由计算得到的压力值，根据公式(25)计算速度的修正值；式中， u′i——i粒子的速度修正值；再由计算得到的速度修正值，根据公式(26)，公式(27)分别求解得到粒子在下一时层的速度及位置； rin+1＝ri*+u′iΔt 公式(27) 步骤7：输出计算所得到的高瑞利数熔融池中粒子的种类、速度位置压力及温度值Tin+1；并对熔融池内部粒子信息进行处理，按照不同的粒子种类表示不同的物质即得到熔融池内的硬壳分布，而对于壁面热流密度根据公式(28)求解；式中， q——穿过壁面的热流密度/W/m2； κwall——壁面材料的导热系数/W/(m·K)； ——A、B粒子在n+1时层的温度K，A、B粒子分别为位于同一径向的内壁面粒子和外壁面粒子； ——A、B粒子在n+1时层的位置矢量m。
所属类别：	发明专利