CN114510677B - 基于间断有限元的中子输运方程处理方法、计算机程序产品 - Google Patents
基于间断有限元的中子输运方程处理方法、计算机程序产品 Download PDFInfo
- Publication number
- CN114510677B CN114510677B CN202210062137.4A CN202210062137A CN114510677B CN 114510677 B CN114510677 B CN 114510677B CN 202210062137 A CN202210062137 A CN 202210062137A CN 114510677 B CN114510677 B CN 114510677B
- Authority
- CN
- China
- Prior art keywords
- finite element
- transport equation
- neutron transport
- unit
- function
- 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.)
- Active
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 17
- 238000004590 computer program Methods 0.000 title claims description 10
- 238000000034 method Methods 0.000 claims abstract description 53
- 238000006243 chemical reaction Methods 0.000 claims abstract description 21
- 239000011159 matrix material Substances 0.000 claims abstract description 19
- 239000000243 solution Substances 0.000 claims description 26
- 230000004907 flux Effects 0.000 claims description 15
- 238000013507 mapping Methods 0.000 claims description 11
- 238000012545 processing Methods 0.000 claims description 8
- 238000012360 testing method Methods 0.000 claims description 8
- 239000012086 standard solution Substances 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 31
- 238000013461 design Methods 0.000 abstract description 6
- 238000004458 analytical method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 230000010355 oscillation Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 125000003275 alpha amino acid group Chemical group 0.000 description 2
- 150000001875 compounds Chemical class 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000035755 proliferation Effects 0.000 description 2
- 235000014676 Phragmites communis Nutrition 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000004992 fission Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 239000012088 reference solution Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Monitoring And Testing Of Nuclear Reactors (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种基于间断有限元的中子输运方程处理方法、程序产品,处理方法包括构造区间B样条小波单元的形状函数和转换矩阵;在离散纵标中子输运方程基础上,推导常规间断有限元的离散纵标中子输运方程;在常规间断有限元的离散纵标中子输运方程基础上,根据构造的区间B样条小波单元的形状函数和转换矩阵,推导出区间B样条小波间断有限元的离散纵标中子输运方程;求解区间B样条间断有限元的离散纵标中子输运方程;根据求解结果指导核反应堆的设计。能够解决目前常规间断有限元求解中子输运方程时,若网格多则计算耗时,若网格少又计算精度不高、收敛速度慢的技术问题。
Description
技术领域
本发明涉及核工程中子输运方程求解数值模拟领域,具体涉及一种基于间断有限元的中子输运方程处理方法、计算机程序产品。
背景技术
中子输运方程是小型复杂反应堆设计及屏蔽计算的基础和关键,其解也是反应堆及屏蔽计算设计的重要参数。中子输运方程是一个含有7个自变量(包含3个空间位置自变量,2个角度自变量,1个能量自变量和1个时间自变量)的积分-微分方程。然而要精确的求解这个方程基本上是不可能的,即便在稳态情况下,求解中子输运方程依然是很困难的。随着计算机的发展,数值计算在工程领域得到了非常广泛的应用,同时在中子输运方程的求解中也得到了广泛的应用。
有限元方法是现代科学与工程计算领域中最重要的数值方法之一,在中子输运方程求解领域也得到了广泛的研究和应用。间断有限元方法是传统连续有限元方法的创新形式、改进和发展。间断有限元方法起源于一阶双曲问题和椭圆边值问题两类不同问题数值方法的研究,具有以下优点:(1)是物理守恒性质可在单元上满足;(2)数值解具有高分辨率;(3)适用于处理不同介质的物理问题(解在区域内部有间断性);(4)容易实现局部网格加密和各单元多项式的独立选取等优点。虽然间断有限元方法具有以上优点,但是,目前该方法在中子输运方程求解中仍存在一些问题:(1)低阶线性间断有限元相比于高阶间断有限元需要更多的计算网格;(2)高阶间断有限元相比于线性间断有限元在每个单元上需要花费更多的计算时间,同时,可能出现数值“龙格”振荡现象,影响计算的精度;(3)相比于连续有限元,间断有限元方法需要更多的自由度,使得计算效率降低。
小波分析是近年来发展起来的一种全新的数值分析方法,是与有限元方法联合起来进行数值求解的一种方法,与传统的拉格朗日基函数相比,小波函数具有多尺度和多分辨析的特点,能够提供多种基函数作为有限元的插值函数,由此构造的小波基单元可以根据实际需要任意改变分析尺度,使得在变化梯度小的求解域能够用大的分析尺度,在变化梯度大的求解域能够用小的分析尺度。提供了另一种除了单元网格加密和插值阶次升高的自适应求解算法,使得求解的数值稳定性好、运算速度快、求解精度高。目前,在研究的小波主要有Daubechies小波、第二代小波和区间B样条小波(BSWI),而BSWI被认为是所有小波函数中最好的一种,主要是因为,BSWI除了具有小波的特点之外,它还具有局部化性质,能克服求解边值问题出现的数值振荡;另外,BSWI具有显示的表达式,能够极大的提高数值计算效率。因此,基于BSWI这样的一些特点,在连续有限元方面开展了很多相关的研究。但在间断有限元和中子输运方程求解领域中,还缺少采用区间B样条小波函数开展研究和应用。
发明内容
为了克服目前构建小型复杂反应堆时,采用间断有限元方法处理中子输运方程,存在低阶线性间断有限元相比于高阶间断有限元需要更多的计算网格、高阶间断有限元相比于线性间断有限元在每个单元上需要花费更多的计算时间,且可能出现数值“龙格”振荡现象,影响计算的精度,以及相比于连续有限元,间断有限元方法需要更多的自由度,使得计算效率降低,而区间B样条小波函数虽相应具有很多优势,但缺少在间断有限元和中子输运方程求解中开展研究和应用的技术问题,提供一种基于间断有限元的中子输运方程处理方法、计算机程序产品,在网格较少时,计算精度能够大幅提高、收敛速度快,即使网格多,也能快速计算,耗时短。
为实现上述目的,本发明提供如下技术方案:
一种基于间断有限元的中子输运方程处理方法,其特殊之处在于,包括以下步骤:
S1,构造区间B样条小波单元的形状函数和转换矩阵;
S2,根据离散纵标中子输运方程,推导常规间断有限元的离散纵标中子输运方程;
S3,在常规间断有限元的离散纵标中子输运方程基础上,根据构造的区间B样条小波单元的形状函数和转换矩阵,推导出区间B样条小波间断有限元的离散纵标中子输运方程;
S4,求解区间B样条间断有限元的离散纵标中子输运方程;
S5,根据区间B样条间断有限元的离散纵标中子输运方程的求解结果,构建核反应堆。
进一步地,所述步骤S1具体为:
S1.1,根据边值问题维数,确定其相应维数方程;
S1.2,通过网格剖分,将所述维数方程的求解域划分为多个子域,并将每个子域均映射为标准求解域[0,1];
S1.3,以区间B样条尺度函数作为插值函数构造单元,并根据B样条尺度函数的阶数和尺度将单元分为多个部分,并以此确定单元节点总数;
S1.4,确定每个单元节点在标准求解域[0,1]中的映射值;
S1.5,确定未知场函数和剖分单元节点处的场函数值;
S1.6,根据经步骤S1.5得到的未知场函数和剖分单元节点处的场函数值,得到区间B样条小波单元的形状函数和转换矩阵。
进一步地,所述步骤S2具体为:
S2.1,根据稳态多群离散纵标中子输运方程,得到单群中子输运方程;
S2.2,根据Galerkin间断有限元理论,在单群中子输运方程的两侧同时乘以一个试函数,并在单元空间域内进行积分,得到常规间断有限元的离散纵标中子输运方程
进一步地,所述步骤S3具体为:
将形状函数的行形式写为列形式,并带入经步骤S2得到的常规间断有限元的离散纵标中子输运方程,得到区间B样条小波间断有限元的离散纵标中子输运方程
进一步地,所述步骤S4具体为:
S4.1,设定边界条件,并根据实际待解决问题,选择源项;
S4.2,根据相应的收敛准则,求得区间B样条间断有限元的离散纵标中子输运方程的解。
进一步地,所述边界条件包括真空边界条件和全反射边界条件;
所述源项包括固定源和临界源。
进一步地,步骤S1.1中,所述边值问题维数为一维,确定其一维方程为:
L(u(x))=f(x),Ω={xx∈[c,d]}
其中,L表示微分算子,Ω为求解域,u为关于x的场函数,x为变量,c为x的下限值,d为x的上限值;
步骤S1.2具体为,通过网格剖分,将求解域Ω划分为多个子域Ωi,i=1,2,...,将任一子域映射为标准求解域Ωs={ξ|ξ∈[0,1]},其中,ξ为标准单元中的变量;
步骤S1.3具体为,以区间B样条尺度函数作为插值函数构造单元,并根据B样条尺度函数的阶数m和尺度j将单元分为n个部分,
n=2j+m-2;
确定单元节点总数为n+1;
步骤S1.4具体为,确定各单元节点实际坐标值xh为:
xh∈[x1,xn+1],1≤h≤n+1;
定义转换公式:ξ=(x-x1)/le,0≤ξ≤1,其中,le为实际单元的长度
将各单元节点的实际坐标值xh代入转换公式,得到各单元节点实际坐标值的映射值ξh:
ξh=(xh-x1)/le,0≤ξh≤1,1≤h≤n+1;
步骤S1.5具体为,确定未知场函数为:
其中,u(ξ)为场函数;i为尺度函数;为小波差值系数列向量,/>为m阶j尺度的尺度函数行向量,/>为小波系数,/>为小波尺度函数;
确定剖分单元节点处的场函数值ue为;
ue=[u(ξ1)u(ξ2)···u(ξn+1)]T
其中,u(ξ1)为单元ξ1节点的未知场函数值,u(ξ2)为单元ξ2节点的未知场函数值,u(ξn+1)为单元ξn+1节点的未知场函数值;
步骤S1.6具体为,将各单元节点实际坐标值的映射值ξh代入未知场函数中的ξ,得到:
ue=Reae;
其中,Re=[ΦT(ξ1)ΦT(ξ2)···ΦT(ξn+1)]T,ΦT(ξ1)为Φ的转置矩阵;
再将ue代入步骤S1.5中的未知场函数,得到:
u(ξ)=Φ(Re)-1ue=Neue;
其中,Ne=Φ(Re)-1=ΦTe,Ne为区间B样条函数单元的形状函数,Te为区间B样条尺度函数单元的转换矩阵。
进一步地,步骤S2.1具体为:
确认稳态多群离散纵标中子输运方程为:
其中,r表示空间位置,k=1,···,M,g=1,···,G,M表示中子飞行方向总数,G表示能群总数;为第k个离散方向的第g群中子飞行方向,Σt g为第g群宏观中子总截面,为第k个离散方向的第g群中子角通量密度,Sg(r)为第g群的源项;
得到单群中子输运方程为:
其中,Ωk为第k个离散方向的中子飞行方向,Σt为宏观中子总截面,为第k个离散方向的中子角通量密度,S(r)为源项;
步骤S2.2具体为:根据Galerkin间断有限元理论,在单群中子输运方程的两侧同时乘以一个试函数v;
在单元空间域内进行积分,得到:
根据散度定理有:
其中,
Ve为单元空间域,Γe为单元空间域Ve的边界,n为边界单位外法向量,入射边界表示为Γe -,出射边界表示为Γe +;
得到:
对中子角通量密度和试函数v应用常规拉格朗日的形状函数进行展开,得到常规间断有限元的离散纵标中子输运方程:
进一步地,步骤S4中,所述边界条件为真空边界条件或反射边界条件;所述实际待解决问题包括固定源问题和临界源问题。
本发明还提供了一种计算机程序产品,包括计算机程序,其特殊之处在于,该程序被处理器执行时实现上述方法的步骤。
与现有技术相比,本发明的有益效果是:
1.本发明基于间断有限元的中子输运方程处理方法,与现有处理中子输运方程的方法相比,本发明的方法具有收敛速度快和计算精度高的特点。
2.本发明的处理方法中,在处理中子输运方程过程时,相同物理网格下,能够达到2阶或者更高阶常规间断有限元的精度。
3.本发明的处理方法中,在处理中子输运方程过程时,能够避免高阶常规间断有限元的“龙格”数值不稳定现象。
附图说明
图1为本发明基于间断有限元的中子输运方程处理方法实施例的流程示意图;
图2为本发明实施例一中的ISSA临界基准问题几何示意图;
图3为本发明实施例一与通过其他不同方法的计算结果示意图;
图4为同网格下本发明实施例一与采用常规间断有限元方法求解中子输运方程的计算精度对比图;
图5为本发明实施例二中的固定源问题几何示意图;
图6为本发明实施例二与其他不同方法的计算结果示意图;
图7为本发明实施例二与采用其他不同阶常规间断有限元方法求解中子输运方程的结果比较图;
图8为本发明中子域Ωi的离散示意图。
具体实施方式
下面将结合本发明的实施例和附图,对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例并非对本发明的限制。
本发明提供了一种计算中子输运方程的新型间断有限元处理方法,用于指导核反应堆的设计,主要内容包括:1、构造区间B样条小波单元的形状函数和转换矩阵;2、在离散纵标中子输运方程基础上,推导常规间断有限元的离散纵标中子输运方程;3、在常规间断有限元的离散纵标中子输运方程基础上,根据构造的区间B样条小波单元的形状函数和转换矩阵,推导出区间B样条小波间断有限元的离散纵标中子输运方程;4、求解区间B样条间断有限元的离散纵标中子输运方程。5、根据求解结果指导核反应堆的构建。
首先,根据***的边界条件,针对不同类型的问题,选择不同的源项,即固定源或者临界源,再根据不同的收敛准则,得到方程的解。
如图1,下面具体给出获得中子输运方程解的新型间断有限元方法:
第一步,构造区间B样条小波单元的形状函数和转换矩阵。
以一维边值问题为例进行说明。该一维方程为:
L(u(x))=f(x),Ω={x|x∈[c,d]} (1)
(1)式中,L表示微分算子,Ω为求解域,通过网格剖分将Ω分成许多子域Ωi,i=1,2,…对任一子域Ωi={x|x∈[a,b]},可以映射为标准的求解域Ωs={ξ|ξ∈[0,1]}。当采用阶数为m,尺度为j的区间B样条尺度函数(记为BSWImj)作为插值函数构造单元时,单元上的节点分布及相应的坐标如图8所示,Ωi(单元长度le=xn+1-x1)被分成n=2j+m-2个部分,单元节点总数为n+1。单元各节点实际坐标值为:
xh∈[x1,xn+1],1≤h≤n+1 (2)
定义转换公式:
ξ=(x-x1)/le,0≤ξ≤1 (3)
式(3)将x映射成标准求解区间[0,1],将(2)式代入(3),可得到每个节点xh的映射值ξh:
ξh=(xh-x1)/le,0≤ξh≤1,1≤h≤n+1 (4)
当采用BSWImj尺度函数作为插值函数时,未知场函数u(ξ)可表示为:
其中:
表示小波插值系数列向量;
表示m阶j尺度的尺度函数行向量。
根据式(4),可得剖分单元节点处的场函数值为:
ue=[u(ξ1)u(ξ2)···u(ξn+1)]T (6)
然后将式(4)代入式(5),可得:
ue=Reae (7)
式(7)中,矩阵Re为:
Re=[ΦT(ξ1)ΦT(ξ2)···ΦT(ξn+1)]T (8)
将(7)代入(5)中,可得
u(ξ)=Φ(Re)-1ue=Neue (9)
式(9)中,
Ne=Φ(Re)-1=ΦTe (10)
其中
Te=(Re)-1 (11)
式(10)中Ne表示BSWImj单元的形状函数,Te表示BSWImj单元的转换矩阵。
第二步,常规间断有限元的离散纵标中子输运方程推导,在离散纵标中子输运方程基础上,推导出常规间断有限元的离散纵标中子输运方程。
稳态多群离散纵标中子输运方程为:
其中,k表示中子飞行方向数,g表示能群数。Ωg为第g群中子飞行方向,Σt g为第g群宏观中子总截面,为第g群中子角通量密度,S(r)为方程的源项,包括散射源、固定源和临界源。对于各项同性源,可表示为:
式(13)中,为第g群的角中子通量密度,qg为第g群固定源,/>为第g'群到第g群的散射截面,υg为第g群的平均裂变中子数,χg为第g群的裂变中子份额,/>为***有效增值系数,/>为第g群的裂变截面,式(14)中,φg为第g群标通量,wk为离散纵标方程的求积权重。
常规间断有限元的离散纵标中子输运方程,由于单群中子输运方程解法和多群解法的一致性,本发明基于单群离散纵标中子输运方程进行弱解方程的推导,并且忽略能群标识g。基于式(12)的多群离散纵标中子输运方程,给出单群中子输运方程为:
根据Galerkin间断有限元理论,将方程(19)两边同时乘以试函数v(r),并在单元空间域内进行积分,得到方程为:
根据散度定理有:
其中,
其中,式(17)、式(18)中,Ve为单元空间域,Γe为空间域Ve的边界,n为边界单位外法向量,入射边界表示为Γe -,出射边界表示为Γe +,并将此2式代入方程(16)中,可得到:
为了进一步求解方程(19),需要对中子角通量密度和试函数v应用常规拉格朗日的形状函数进行展开,然后代入方程(19)中有:
方程(20)即是中子输运方程的弱形式,也即常规间断有限元的离散纵标中子输运方程。
第三步,推导区间B样条小波间断有限元的离散纵标中子输运方程,在第二步常规间断有限元的离散纵标中子输运方程基础上,根据第一步中构造的区间B样条小波单元的形状函数和转换矩阵,推导出区间B样条小波间断有限元的离散纵标中子输运方程。
为了将区间B样条小波形状函数应用于常规间断有限元中,首先需要将形状函数的行形式写成列形式,即:
N=(Ne)T (21)
然后将式(21)代入方程(20)中,可得到:
方程(22)即是基于区间B样条间断有限元的离散纵标中子输运方程。方程(22)和方程(20)最大的不同在于使用了不同类型的形状函数,前者使用了区间B样条小波函数所构造的形状函数,后者使用了常规拉格朗日型的形状函数,不同形状函数的使用导致方程的收敛性发生极大的变化。
第四步,求解第三步中推导的区间B样条间断有限元的离散纵标中子输运方程。首先需给出***的边界条件,然后针对不同类型的问题,选择不同的源项,即固定源或者临界源,最后根据不同的收敛准则,得到方程的解。
边界条件,对于中子输运方程的求解***,一般要么是真空边界条件要么全是反射边界条件,表示如下:
真空边界条件:
全反射边界条件:
其中,m和m'是关于方向对称的。
对于固定源问题,在求解过程中,满足逐点中子通量密度收敛准则即可,该准则如下所示:
其中,分别为第k-1次迭代中子角通量密度和第k次迭代中子角通量密度,EPF一般满足1×10-4。只需满足这样的一次迭代收敛过程,就达到求解的要求,这种迭代为内迭代。
对于临界源问题求解,除了要求满足逐点中子通量密度收敛准则外,还需满足特征值收敛准则,表示如下:
其中,keffk-1、keffk分别为第k-1次迭代***有效增殖系数和第k次迭代***有效增殖系数,EPK一般满足1×10-5,其中,
而式(27)中,Qf为裂变源,表示为:
Qf=υΣfφ(r) (28)
这种特征值的收敛,一般为外迭代,角中子通量密度的收敛为内迭代。即是,根据不同的问题选择不同的迭代策略,最终得到问题的精确解。
实施例一
如图2所示为ISSA临界基准问题,该问题题为单群问题,包括两个材料区,右边界为真空,左边界为全反射。
根据前述详细技术方案,首先,根据第一步到第三步给出的公式进行处理,然后,根据问题类型,选择源项计算公式(13)中的临界源,再次根据边界条件公式(23)和公式(24),最后根据临界问题的收敛准则公式(25)和公式(26),最终求解公式(22),得到问题的解。图3给出了不同程序的计算结果,表明达到同样的计算精度,本发明方法所用网格最少,且精度最高,计算时间最少,收敛速度快。图4给出了在相同网格条件下,本发明方法计算精度达到了2阶常规间断有限元的计算精度。其中,图3中,mesh为计算物理网格数,Sn-ARES为常规差分法,T-DGFEM为常规间断有限元方法,SPARK为区间B样条小波间断有限元方法,即本发明方法,Reference为参考解。图4中,p1-T-DGFEM为1阶常规间断有限元方法,p2-T-DGFEM为2阶常规间断有限元方法,p3-T-DGFEM为1阶常规间断有限元方法,SPARK为区间B样条小波间断有限元方法,即本发明方法。
实施例二
如图5为单群固定源问题(Reed cell problem),该问题是一个单群5区固定源问题,用来验证在复杂条件下的计算精度问题,该问题左边界为全反射边界条件,右边界为真空边界。
基于前述的具体方法,根据问题类型,选择源项计算公式(13)中的固定源,再根据边界条件公式(23)和公式(24),最后根据固定源问题的收敛准则公式(25),最终求解公式(22),得到该问题的解。图6给出了不同程序的计算结果,其中,达到相同精度,本发明方法SPARK只用了16个网格,而常规间断有限元方法T-DGFEM用了160个物理网格,Exact为精确解,本发明方法,即SPARK极大的减少了物理网格的计算量。图7给出了本发明方法与不同阶常规间断有限元计算结果比较,由图可以看到,在相同网格条件下,本发明方法的计算精度达到了3阶常规间断有限元的计算精度。
本发明中方法处理中子输运方程后,可以为复杂反应堆的设计及屏蔽设计提供一种较快速的处理方法。
本发明的中子输运方程处理方法也可作为计算机程序产品,包括计算机程序,程序被处理器执行时实现上述处理方法的步骤。
以上所述仅为本发明的实施例,并非对本发明保护范围的限制,凡是利用本发明说明书及附图内容所作的等效结构变换,或直接或间接运用在其他相关的技术领域,均包括在本发明的专利保护范围内。
Claims (8)
1.一种基于间断有限元的中子输运方程处理方法,其特征在于,包括以下步骤:
S1,构造区间B样条小波单元的形状函数和转换矩阵,具体为:
S1.1,根据边值问题维数,确定其相应维数方程;
S1.2,通过网格剖分,将所述维数方程的求解域划分为多个子域,并将每个子域均映射为标准求解域[0,1];
S1.3,以区间B样条尺度函数作为插值函数构造单元,并根据B样条尺度函数的阶数和尺度将单元分为多个部分,并以此确定单元节点总数;
S1.4,确定每个单元节点在标准求解域[0,1]中的映射值;
S1.5,确定未知场函数和剖分单元节点处的场函数值;
S1.6,根据经步骤S1.5得到的未知场函数和剖分单元节点处的场函数值,得到区间B样条小波单元的形状函数和转换矩阵;
S2,根据离散纵标中子输运方程,推导常规间断有限元的离散纵标中子输运方程,具体为:
S2.1,根据稳态多群离散纵标中子输运方程,得到单群中子输运方程;
S2.2,根据Galerkin间断有限元理论,在单群中子输运方程的两侧同时乘以一个试函数,并在单元空间域内进行积分,得到常规间断有限元的离散纵标中子输运方程;
S3,在常规间断有限元的离散纵标中子输运方程基础上,根据构造的区间B样条小波单元的形状函数和转换矩阵,推导出区间B样条小波间断有限元的离散纵标中子输运方程;
S4,求解区间B样条间断有限元的离散纵标中子输运方程;
S5,根据区间B样条间断有限元的离散纵标中子输运方程的求解结果,做为构建核反应堆的参数。
2.如权利要求1所述基于间断有限元的中子输运方程处理方法,其特征在于,所述步骤S3具体为:
将形状函数的行形式写为列形式,并带入经步骤S2得到的常规间断有限元的离散纵标中子输运方程,得到区间B样条小波间断有限元的离散纵标中子输运方程。
3.如权利要求2所述基于间断有限元的中子输运方程处理方法,其特征在于,所述步骤S4具体为:
S4.1,设定边界条件,并根据实际待解决问题,选择源项;
S4.2,根据相应的收敛准则,求得区间B样条间断有限元的离散纵标中子输运方程的解。
4.如权利要求3所述基于间断有限元的中子输运方程处理方法,其特征在于:
步骤S4.1中,所述边界条件包括真空边界条件和全反射边界条件;
步骤S4.2中,所述源项包括固定源和临界源。
5.如权利要求4所述基于间断有限元的中子输运方程处理方法,其特征在于:
步骤S1.1中,所述边值问题维数为一维,确定其一维方程为:
L(u(x))=f(x),Ω={x|x∈[c,d]}
其中,L表示微分算子,Ω为求解域,u为关于x的场函数,x为变量,c为x的下限值,d为x的上限值;
步骤S1.2具体为,通过网格剖分,将求解域Ω划分为多个子域Ωi,i=1,2,...,将任一子域映射为标准求解域Ωs={ξ|ξ∈[0,1]},其中,ξ为标准单元中的变量;
步骤S1.3具体为,以区间B样条尺度函数作为插值函数构造单元,并根据B样条尺度函数的阶数m和尺度j将单元分为n个部分,
n=2j+m-2;
确定单元节点总数为n+1;
步骤S1.4具体为,确定各单元节点实际坐标值xh为:
xh∈[x1,xn+1],1≤h≤n+1;
定义转换公式:ξ=(x-x1)/le,0≤ξ≤1,其中,le为实际单元的长度
将各单元节点的实际坐标值xh代入转换公式,得到各单元节点实际坐标值的映射值ξh:
ξh=(xh-x1)/le,0≤ξh≤1,1≤h≤n+1;
步骤S1.5具体为,确定未知场函数为:
其中,u(ξ)为场函数;i为尺度函数;为小波差值系数列向量,/>为m阶j尺度的尺度函数行向量,/>为小波系数,/>为小波尺度函数;
确定剖分单元节点处的场函数值ue为;
ue=[u(ξ1) u(ξ2)…u(ξn+1)]T
其中,u(ξ1)为单元ξ1节点的未知场函数值,u(ξ2)为单元ξ2节点的未知场函数值,u(ξn+1)为单元ξn+1节点的未知场函数值;
步骤S1.6具体为,将各单元节点实际坐标值的映射值ξh代入未知场函数中的ξ,得到:
ue=Reae;
其中,Re=[ΦT(ξ1) ΦT(ξ2)…ΦT(ξn+1)]T,ΦT(ξ1)为Φ的转置矩阵;
再将ue代入步骤S1.5中的未知场函数,得到:
u(ξ)=Φ(Re)-1ue=Neue;
其中,Ne=Φ(Re)-1=ΦTe,Ne为区间B样条函数单元的形状函数,Te为区间B样条尺度函数单元的转换矩阵。
6.如权利要求5所述基于间断有限元的中子输运方程处理方法,其特征在于:
步骤S2.1具体为:
确认稳态多群离散纵标中子输运方程为:
其中,r表示空间位置,k=1,…,M,g=1,…,G,M表示中子飞行方向总数,G表示能群总数;为第k个离散方向的第g群中子飞行方向,Σt g为第g群宏观中子总截面,/>为第k个离散方向的第g群中子角通量密度,Sg(r)为第g群的源项;
得到单群中子输运方程为:
其中,Ωk为第k个离散方向的中子飞行方向,Σt为宏观中子总截面,为第k个离散方向的中子角通量密度,S(r)为源项;
步骤S2.2具体为:根据Galerkin间断有限元理论,在单群中子输运方程的两侧同时乘以一个试函数v;
在单元空间域内进行积分,得到:
根据散度定理有:
其中,
Ve为单元空间域,Γe为单元空间域Ve的边界,n为边界单位外法向量,入射边界表示为Γe -,出射边界表示为Γe +;
得到:
对中子角通量密度和试函数v应用常规拉格朗日的形状函数进行展开,得到常规间断有限元的离散纵标中子输运方程:
7.如权利要求6所述基于间断有限元的中子输运方程处理方法,其特征在于:
步骤S4中,所述边界条件为真空边界条件或反射边界条件;所述实际待解决问题包括固定源问题和临界源问题。
8.一种计算机程序产品,包括计算机程序,其特征在于:该程序被处理器执行时实现权利要求1至7任一所述方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210062137.4A CN114510677B (zh) | 2022-01-19 | 2022-01-19 | 基于间断有限元的中子输运方程处理方法、计算机程序产品 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210062137.4A CN114510677B (zh) | 2022-01-19 | 2022-01-19 | 基于间断有限元的中子输运方程处理方法、计算机程序产品 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114510677A CN114510677A (zh) | 2022-05-17 |
CN114510677B true CN114510677B (zh) | 2024-06-11 |
Family
ID=81550026
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210062137.4A Active CN114510677B (zh) | 2022-01-19 | 2022-01-19 | 基于间断有限元的中子输运方程处理方法、计算机程序产品 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114510677B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115935770B (zh) * | 2022-12-16 | 2023-07-25 | 西安交通大学 | 针对核反应堆屏蔽设计的基于神经网络的中子输运方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH07140289A (ja) * | 1993-11-17 | 1995-06-02 | Toshiba Corp | 中性子輸送シミュレータ |
CN107273582A (zh) * | 2017-05-23 | 2017-10-20 | 西安交通大学 | 一种用于快中子反应堆中子输运燃耗耦合分析的计算方法 |
CN108733903A (zh) * | 2018-05-08 | 2018-11-02 | 西安交通大学 | 一种按中子能量进行个性化处理的中子输运数值模拟方法 |
CN109657408A (zh) * | 2019-01-10 | 2019-04-19 | 上海索辰信息科技有限公司 | 一种再生核粒子算法实现结构线性静力学仿真方法 |
KR20210059194A (ko) * | 2019-11-15 | 2021-05-25 | 울산과학기술원 | 확률론적 방법을 이용한 다군 군정수 생성방법 및 생성장치 |
CN113673116A (zh) * | 2021-09-01 | 2021-11-19 | 上海交通大学 | 针对均匀几何变分节块法的三维准输运加速方法 |
-
2022
- 2022-01-19 CN CN202210062137.4A patent/CN114510677B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH07140289A (ja) * | 1993-11-17 | 1995-06-02 | Toshiba Corp | 中性子輸送シミュレータ |
CN107273582A (zh) * | 2017-05-23 | 2017-10-20 | 西安交通大学 | 一种用于快中子反应堆中子输运燃耗耦合分析的计算方法 |
CN108733903A (zh) * | 2018-05-08 | 2018-11-02 | 西安交通大学 | 一种按中子能量进行个性化处理的中子输运数值模拟方法 |
CN109657408A (zh) * | 2019-01-10 | 2019-04-19 | 上海索辰信息科技有限公司 | 一种再生核粒子算法实现结构线性静力学仿真方法 |
KR20210059194A (ko) * | 2019-11-15 | 2021-05-25 | 울산과학기술원 | 확률론적 방법을 이용한 다군 군정수 생성방법 및 생성장치 |
CN113673116A (zh) * | 2021-09-01 | 2021-11-19 | 上海交通大学 | 针对均匀几何变分节块法的三维准输运加速方法 |
Non-Patent Citations (1)
Title |
---|
基于 MOOSE 的移动小堆中子学求解程序开发验证与应用;胡田亮 等;中国核科学技术进展报告;20211019;第1-9页 * |
Also Published As
Publication number | Publication date |
---|---|
CN114510677A (zh) | 2022-05-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Coirier et al. | An accuracy assessment of Cartesian-mesh approaches for the Euler equations | |
Sherbrooke et al. | Computation of the solutions of nonlinear polynomial systems | |
CN105868853B (zh) | 一种短期风电功率组合概率预测方法 | |
Yang et al. | MOCUM: A two-dimensional method of characteristics code based on constructive solid geometry and unstructured meshing for general geometries | |
CN114510677B (zh) | 基于间断有限元的中子输运方程处理方法、计算机程序产品 | |
Welch et al. | Isogeometric analysis for the multigroup neutron diffusion equation with applications in reactor physics | |
Ferrer et al. | The linear source approximation and particle conservation in the Method of Characteristics for isotropic and anisotropic sources | |
CN111177965B (zh) | 基于定常问题求解的多重分辨weno格式定点快速扫描方法 | |
CN114510775B (zh) | 一种复杂模型三维空间曲网格划分方法 | |
Hanson | Techniques in computational stochastic dynamic programming | |
Abedian | A symmetrical WENO-Z scheme for solving Hamilton–Jacobi equations | |
CN111695281A (zh) | 一种四面体网格划分有限元粒子模拟的粒子快速定位方法 | |
Eidnes et al. | Adaptive energy preserving methods for partial differential equations | |
CN112434451B (zh) | 一种基于分块并行计算的有限元分析方法 | |
Miao et al. | Local modified mesh deformation based on radial basis functions for fluid solid interaction in reactor core | |
Coirier et al. | A Cartesian, cell-based approach for adaptively-refined solutions of the Euler and Navier-Stokes equations | |
Merchel et al. | Fast formation of matrices for least-squares fitting by tensor-product spline surfaces | |
CN115272594A (zh) | 一种基于geotools的等值面生成方法 | |
CN113532308A (zh) | 数字图像相关中带初值的岭回归应变测量方法 | |
Bennett et al. | A surrogate model based on sparse grid interpolation for boiling water reactor subchannel void distribution | |
Moore et al. | Efficient energy evaluations for active B-Spline/NURBS surfaces | |
Zheng et al. | Track truncation method for improving the computational efficiency of 2D matrix MOC | |
Gong | Immersed-interface finite-element methods for elliptic and elasticity interface problems | |
Dominesey | Reduced-Order Modeling of Neutron Transport by Proper Generalized Decomposition | |
Faber et al. | Computer plotting of contour maps: An improved method |
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 |