一种基于FETI的高精细流致振动模拟方法
一种基于feti的高精细流致振动模拟方法
技术领域
1.本发明属于反应堆流致振动数值模拟实现领域，具体涉及一种基于feti的大规模流致振动高精度模拟方法。
背景技术：
2.反应堆堆芯组件具有多样性和复杂性，全堆高精细结构力学模拟所需的计算资源和存储在小型集计算系统中难以实现。依托高性能计算技术来实现高精细的数值计算就成为研究反应堆安全性的必要途径，依托超级计算机，采用数学物理计算方法，对反应堆堆中子物理、结构力学、燃料性能等物理过程进行精细化数值模拟，已成为国际上一个前沿研究热点。
3.燃料棒组件在高温、高压以及强辐照下运行，并且长期承受冷却剂的冲刷。冷却剂的循环流动会引起燃料棒组件的长期微幅振动，且总是伴随着反应堆的运行而存在，这些振动使燃料棒在与格架接触的界面上产生相对位移，并在支撑处发生包壳磨损(grid-to-rod fretting，gtrf)。燃料棒的这种渐进磨损是影响燃料组件在正常运行和事故工况下结构完整性的关键因素，根据全世界燃料棒失效率的统计数据，gtrf磨损失效事故占比高达78％。为确保结构在使用寿命期内的完整性，只能通过优化结构设计和调整流速，将流致振动引起的结构振动响应控制在可接受的范围之内。因此，燃料组件流致振动数值分析与计算是安全分析不可缺少的重要内容。
4.高精细的分析和模拟燃料棒的流致振动，如对于第四代快堆带绕丝的燃料棒的高精度数值模拟，往往需要生成千万单元以上规模的网格数据。考虑湍流模型时所需网格规模更是达到了数亿乃至上十亿网格节点，这对计算机的硬件配置、计算算法及其前处理过程的时空效率都提出了严峻的挑战。使用串行算法求解问题，在内存空间和时间上都变的不可接受。以开源软件code_aster为例，仅对一螺距单根带绕丝燃料棒(约为10万网格单元)的流致振动模拟采用单核求解，其用时需达2小时。以此估算，针对亿量级单元的串行求解，在内存空间和时间上都变的不可接受。当前的研究大多针对小规模算例，并未对千万单元以上规模算例进行分析。利用高性能计算机，开发高效并行的多领域数值模拟求解算法及与之对应的并行前后处理算法是迎接这一挑战的必由之路。开展流致振动并行求解方法研究，是提高流致振动求解效率的关键所在。
5.有限元法是求解复杂工程问题的一种近似数值解法，现已广泛应用到力学、热学、电磁学等各个学科。传统的有限元法一般以串行的方式进行求解，称为串行有限元法。串行有限元法在面对大型复杂工程问题时，往往会出现求解时间过长、内存占用过高等现象。随着支持并行计算的并行计算机问世，出现适用于将串行有限元法并行化的三种技术，分别是方程组并行求解、ebe(element-by-element)有限元法以及有限元区域分解法。基于舒尔补(schur complement)的有限元区域分解法，最早是由farhat等人提出，称为撕裂有限元法(finite element tearing and interconnecting method，feti)，常用于求解线性工程问题中出现的大型线性系统的数值解。在feti方法中，将一个主体分解为几个不重叠的子
域，并且这些子域之间的连续性由拉格朗日乘数强制执行。feti方法使用对偶理论，可以导出较小且条件相对较好的对偶问题，并通过共轭梯度算法的合适变体有效解决。
技术实现要素：
6.为了解决上述问题，本发明提出了一种基于feti的高精细流致振动模拟方法。本发明涉及基于feti方法的newmark积分法模拟高精细流致振动和适用于feti方法的域边界平衡的图二分算法，feti方法的newmark积分法求解高精细流致振动问题过程分成两个阶段：第一阶段大规模并行的feti方法被用于将上亿网格单元的求解域划分为多个非重叠求解子域进行并行求解；第二阶段动力学newmark方法用于实现流致振动过程的时间步更新。
7.本发明方法的具体步骤是：
8.步骤1.读入相应的算例数据，主要以网格数据为主，以及材料参数等。
9.步骤2.采用newmark方法对流致振动动力学过程进行数值离散。
10.步骤3.采用feti方法对离散后的方程进行并行分解,划分子域，子域在划分边界处由拉格朗日乘子λ进行粘合,并满足相应边界条件。根据feti划分子域边界特点,提出了域边界平衡的图二分算法,均衡各个子域中的单元量和计算量,保证进程间负载均衡。
11.步骤4.采用预处理共轭梯度法进行迭代求解,迭代求解得到拉格朗日乘子λ。根据拉格朗日乘子λ得到相应的位移，记录位移值。
12.步骤5.根据位移值，newmark方法实现流致振动过程的时间步更新。
13.本发明有益效果：本发明采用feti方法与newmark积分法相结合，完成了数亿规模的网格数据求解，提高了求解大规模流致振动问题的效率，实现了流致振动的快速高效模拟。采用基于域边界平衡的图二分算法进行区域分解，保证了进程间的负载均衡，进一步提高了求解大规模流致振动问题的效率。
附图说明
14.图1是基于feti方法的newmark积分法求解高精细流致振动问题具体流程；
15.图2是halo权重说明；
16.图3是单次迭代求解时间比例；
17.图4是多线程并行操作求解矩阵向量乘；
18.图5是弱可扩展性测试。
具体实施方式
19.以下结合附图对本发明作进一步说明，请参阅图1；图1给出了本发明求解高精细流致振动问题具体流程。主要分成两个阶段：第一阶段为feti求解阶段，该阶段主要可以分成两个步骤，第一步为采用域边界平衡的图二分算法进行子域划分；第二步为采用预处理共轭梯度法进行迭代求解。第二阶段为newmark方法进行时间步的更新。
20.下面分别对上述实施步骤进行详细说明：
21.步骤1：采用cefr单盒燃料棒组件为例进行数值模拟。燃料棒束底部固支，模拟以流体压力的形式作用在燃料棒和绕丝的外表面上。分布式读取相应网格参数，网格文件格式为vtk格式，以四面体网格和六面体网格为主。设置需要读入的相应材料参数，主要为材
料密度、弹性模量、泊松比、时间步长、划分子域数等。此次模拟材料密度为7890kg/m3、泊松比为0.3、弹性模量为2.1*105mpa，模拟12个时间步长，每个时间步长为0.01秒，划分子域数为512个子域。
22.步骤2：流致振动问题通过对平衡方程、质量矩阵、阻尼矩阵、几何方程、物理方程和位移边界条件采用虚功原理推导得到整体动力方程，用数值算法(newmark算法)根据当前时刻及前几个时刻体系的动力反应值来建立下一时刻的动力反应方程，并进行求解，针对燃料棒流致振动问题建模得到下面的整体动力方程：
[0023][0024]
其中，m为质量矩阵，c为阻尼矩阵，k为刚度矩阵，f(t)是节点载荷，u为求解的燃料棒的位移向量，u的一阶导数为速度，u的二阶导数为加速度。采用newmark方法对该动力学过程进行数值离散，newmark方法针对线性加速度假定做修正，在δ
t
+t时刻的速度和位移表达式中引入两个参数η和γ得到newmark算法的离散格式：
[0025][0026][0027]
将(2)式与(3)式代入原始二阶微分方程(1)式中，可得：
[0028][0029]
将(4)式简化为下面关于求解未知量u
n+1
的线性方程形式，n为时间：
[0030]
au
n+1
＝hnꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0031]
其中矩阵a和矩阵hn分别为
[0032][0033][0034]
步骤3：采用feti方法对离散后的方程进行并行分解，划分子域。在feti方法进行子域划分时，采用本发明提出的基于域边界平衡的图二分算法对求解域进行子域划分，将求解域分解为i个独立的非重叠子域,这些子域在划分边界处由拉格朗日乘子λ进行粘合,并满足边界条件。这样公式5可转化成如下形式：
[0035][0036]
其中,
[0037][0038]
表示分解后的各子域位移向量，bi为一个有符号布尔矩阵，定义子域与相邻子域的连通性，由1、-1和0组成。如果子域x、y连通，设置子域x到子域y的系数为1，子域y到子域x的系数为-1；如果子域x、y不连通，设置系数为0。o为零矩阵。
[0039]
feti方法特殊的边界处理方式，会使得传统的点分割方法产生严重的边划分不均衡现象以及边的缺失。子域边界的节点数量影响着迭代求解效率，构建的求解矩阵计算成本随着边的计算成本的增加而增加。因此，保证每个子域边界的节点数量均衡至关重要。在feti的区域分解时，采用基于域边界平衡的图二分算法不仅可以保证各个子域内部节点均衡,还能确保这个子域边界节点均衡，从而提高了对应的求解效率。
[0040]
在传统的贪心图增长算法基础上，对顶点增加两个权重：对属于边界节点的顶点设置halo权重为1，非边界节点的顶点设置halo权重为0，对子图内部非边界节点的顶点设置halo权重为1，属于边界节点的顶点设置halo权重为0。将halo权重为1的顶点添加至halo点集vh。请参阅图2，左图黑点为分隔符，将初始图划分两张子图，halo权重都为4。右图为两张子图划分为4张子图，子图的halo权重为别为4,4,4,4。
[0041]
如下伪代码给出了域边界平衡的图二分算法具体流程：
[0042][0043]
首先在halo点集vh中选择两个种子顶点ω1和ω2，使得ω1与ω2尽可能保持距离，
以保证ω1与ω2,节点增长较晚的相遇。将初始图划分为两张子图，两张子图的点集分别用p1和p2表示。p1和p2最初为空，由种子顶点ω1和ω2生长出来。从ω1和ω2出发分别对剩下顶点做划分，根据p1和p2两个点集中的halo平衡情况进行调节。如果p1(或p2)中的halo权重为1的顶点比p2(或p1)中少，则尽可能选择halo权重为1的目标顶点加入p1(或p2)。如果具有更多的halo权重为1的顶点，则优先选择非halo权重为1的顶点加入p1(或p2)。考虑完halo平衡后，则p1选取离ω1最近和ω2最远的顶点(dist(v,ω1)-dist(v,ω2)最小值的顶点)。这个过程需要满足以下公式：
[0044][0045][0046]
其中s代表边分隔符，vh代表求解域中的halo点集。保持最小的边分隔符和halo节点均衡可以保证子域边界的节点均衡以及各子域之间边的负载均衡。选择不同的种子顶点进行多次增长，根据最小边分隔符和halo节点均衡进行分割。
[0047]
如表1、表2、表3所示可以发现采用本发明提出的基于域边界平衡的图二分算法与开源软件metis的对比。表1对单元总数为29026600的六面体网格进行划分。表中avg值为理想条件下的划分数，max、min分别代表子域中网格最大数与网格最小数。从表1中可以看出来，本发明提出的方法划分的最大最小子域数更接近理想值。表2对单元总数为102853760的四面体网格进行划分同样证明了本发明提出的方法进行子域划分效果明显较好。本方法对四面体网格算例划分成512、256、128个子域用时与开源图划分软件metis用时对比，两者用时较为接近，子域划分效率接近。
[0048]
表1六面体网格算例子域划分对比
[0049][0050]
表2四面体网格算例子域划分对比
[0051][0052]
表3子域划分效率对比
[0053][0054]
步骤4：矩阵a和矩阵b均为非奇异矩阵，使得未知量λ具有唯一解。根据矩阵广义逆得：
[0055]un+1
＝a
+
(h
n-b
t
λ)+rα s.t.ar＝o or r
t
a＝o
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0056]
其中,矩阵r为矩阵a零空间的基，a
+
是a的伪逆矩阵。将边界条件bu＝0代入至(12)式中可得：
[0057][0058]
这样,原问题转化为对关于未知量λ，α的方程组(13)的数值求解。这里采用预处理共轭梯度法对(13)式进行迭代求解。首先引入预处理矩阵p，对(13)式的第一个式子两边同时左乘矩阵p，用于消除α。从而得到(14)式：
[0059][0060]
其中，
[0061]
p＝i-g(g
t
g)-1gt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0062]
对(14)式采用迭代求解得到λ,再代回(11)式解出α。将λ和α代入(12)式，得到每个求解域的当前时间步的位移向量，保存当前时间步计算结果u
n+1
。
[0063]
请参阅图3，图3给出了单次迭代求解时间比例。程序的计算时间主要由程序的迭代次数和单次迭代时间决定。在迭代求解过程中矩阵求解计算时间占到单次迭代时间的95％。每轮迭代都会执行以下步骤：投影、计算beta、矩阵向量乘、计算alpha、更新结果向量、更新残差向量以及迭代数与迭代精度的判断。
[0064]
基于feti方法求解流致振动算法流程的伪代码如下：
[0065][0066]
伪代码中8-10行为预处理共轭梯度法的迭代求解过程的投影步骤，投影存在大量的矩阵向量乘操作。在各个子域中计算ak、λk以及rk，需要先计算矩阵向量操作f
pk
。矩阵向量乘相互独立，可以进行并行处理。可以将伪代码中的8-10步与17步移植至国产加速器中来提高求解效率。针对曙光先进计算平台架构(具体实验环境如表4)，对矩阵向量乘部分进行了异构移植，采用hip并行编程模型将具有较高计算密度的稠密矩阵向量乘计算交由国产加速器和cpu共同完成。采用openmp并行编程模型进行多线程并行操作。图4给出了多线程并行操作求解矩阵向量乘的过程，一个线程用于控制国产加速器，调用rocblas线性代数库进行稠密矩阵计算。利用线程调用hipmalloc函数，申请显存空间，将存储在cpu上的数据通过hipmemcpy传递给国产加速器显存，通过调用hiplaunchkernel函数启动核函数，再调用hipblas线性代数库进行稠密矩阵计算。在国产加速器上的计算完成后再用hipmemcpy将结果传递给cpu，等待所有设备完成计算，同步数据，再进入下一次迭代。其他线程控制cpu核心，通过调用intel mkl数学库计算剩余稠密矩阵向量。在迭代过程中，为了避免cpu与国产加速器之间频繁数据交互带来的时间开销，在求解开始前就将所有稠密矩阵都传递至国产加速器中。在求解过程中，cpu与国产加速器根据各自任务的编排执行对应任务。
[0067]
表4实验环境
[0068][0069]
在多个计算节点并行计算时，国产加速器需要从内存中读取数据到显存中，在计算结束后，需要将计算结果写入内存。随着问题规模变大，这样的访存操作时间开销会逐渐增大。为充分利用国产加速器计算性能，本发明利用内存和显存之间的双缓冲技术，来减少访存的时间。此外，为合理利用国产加速器显存，采用数据对齐方法以减少显存的不必要占用。在对齐边界中创建数据，确保缓存线大小与向量数据对齐，避免编译器额外计算，从而导致性能的损失。并利用内存合并访问技术，获取更高的吞吐量，以提高计算性能。
[0070]
步骤5：根据迭代求解得出的下一时间步的位移值代入(7)式更新hn，从而进入下一个时间步求解。newmark方法的计算精度取决于时间步长δt的大小。
[0071]
重复上述步骤2至步骤5直至达到所设定的时间步长数(注：域边界平衡的图二分算法在第一个时间步长划分，后续不进行划分)。如图5所示，本发明在曙光先进计算平台架构测试，计算量从0.48亿扩展到7.6亿,计算节点从128个扩展到2048个,每个节点处理64个子域,并行效率达到54％，可以明显发现本发明具有良好的可扩展性，求解时间都在分钟级别，具体实验数据见表5所示。
[0072]
表5弱可扩展性实验
[0073]
