一种基于分块并行计算的有限元分析方法与流程

[0001]
本发明属于工程设计与计算领域，具体涉及一种基于分块并行计算的有限元分析方法。


背景技术：

[0002]
有限元方法作为一种应用于解决工程学和数学物理学的数值方法，被广泛应用于结构分析等工程领域典型问题。通过将结构的特定研究区域离散为数个特定形状的单元，能够将最初描述物理现象的偏微分方程转变为某种形式的矩阵方程，利用特定离散点的计算结果来表征结构内任一点的物理量计算结果，从而简化问题。利用有限元方法对工程问题进行分析，能够为工程实际中问题的研究提供指导，缩短结构设计周期、节约试验成本。
[0003]
然而随着某些工程问题越来越复杂，有限元方法的求解规模也越来越大，传统的开源串行计算方法无论在性能还是求解速度上都难以满足要求，同时也对计算机内存提出了更高的要求。除此之外，随着eigen、mkl等线性代数运算库的飞速发展，如何将有限元过程与现有的功能强大的库函数相结合从而显著减少总运算时间，也成为了一个严峻的问题。因此亟需一种基于高效、便捷、与其他库函数结合方便的有限元并行分析方法，从而大幅缩短有限元方法中大型矩阵的形成时间，减少试验成本。


技术实现要素：

[0004]
本发明的目的在于提供一种基于分块并行计算的有限元分析方法，以解决上述技术问题，该算法主要用来快速并行完成用有限元方法分析的问题中的大规模矩阵的组装过程。该算法利用解释性语言开发，框架基于即时编译器搭建实现即时编译加速，能够高效便捷地实现有限元方法中网格节点单元信息的读取处理、大型矩阵的组装等内容，在保证算法效率的同时兼顾了与其他库函数的交互性、进一步开发的便捷性等。利用并行计算以及一些算法处理技巧，能够大幅缩短整个大规模矩阵的计算时间；同时能够与其他丰富的开源求解器、功能多样的包或是高效并行计算库相结合，从而高效地实现整个有限元分析流程。
[0005]
为了实现上述目的，本发明采用如下技术方案：
[0006]
一种基于分块并行计算的有限元分析方法，包括以下步骤：
[0007]
1)网格信息预处理，导入矩阵包括节点坐标矩阵、单元节点信息矩阵，以及节点位置向量；
[0008]
2)基于步骤1)形成单元节点相关性信息矩阵，得到每个单元包含的节点信息；基于步骤1)形成节点相关性信息矩阵，得到节点与节点是否处于同一单元的信息矩阵；
[0009]
3)基于步骤2)形成节点关联向量n
c
、i
c
和s
c
，得到某节点的关联节点数目、编号、总数信息；
[0010]
4)单元矩阵预处理，得到形成单元矩阵所需的利用有限元法形成矩阵所对应的插值函数，对于结构问题而言，还需要得到总体偏导矩阵d’e
；
[0011]
5)分块组装单元矩阵，对每个单元都并行形成与被研究问题有关的相应矩阵，于结构问题，分块并行形成单元刚度矩阵和单元质量矩阵；
[0012]
6)基于步骤1)中得到的节点位置向量n
ind
，遍历每个单元的节点形成总体矩阵；
[0013]
7)基于步骤6)的组装过程，形成行指针r
p
和列索引c
ind
；
[0014]
8)添加边界条件，对于结构问题，添加位移约束、恒定作用力载荷，然后进行问题求解。
[0015]
本发明进一步的改进在于，步骤1)的具体实现方法如下：
[0016]
导入矩阵包括节点坐标矩阵以及单元节点信息矩阵，其中节点坐标矩阵保存各节点坐标信息，维度为3
×
nnum，纵轴存储各节点x，y，z轴三个方向的坐标，横轴按顺序排列各个结点；单元节点信息矩阵保存各单元节点编号信息，对8节点六面体单元来说，该矩阵为8
×
enum的矩阵；其中nnum为网格节点个数，enum为网格单元个数；另外额外形成一个节点位置向量n
ind
，该索引向量将节点编号与节点序号联系起来，第i个元素存储了编号为i的节点的索引。
[0017]
本发明进一步的改进在于，步骤2)的具体实现方法如下：
[0018]
形成维度为enum
×
nnum的单元节点信息相关性矩阵e
re
＝[e
1τ
,e
2τ
,
…
,e
enumτ
]
τ
，该矩阵基于导出的单元节点信息矩阵利用稀疏矩阵的coo格式得到；其中nnum为离散后的节点个数，enum为离散后得到的单元个数，e
i
存储有第i个单元中的节点相关性信息；
[0019]
形成维度为nnum
×
nnum的节点相关性信息矩阵n，该矩阵表达式为n＝e
t
e；当第i个节点与第j个节点位于同一个单元之中时，元素n
ij
非零；当第i个节点与第j个节点位于同一单元之中时，元素n
ij
为零。
[0020]
本发明进一步的改进在于，步骤3)的具体实现方法如下：
[0021]
n矩阵的第i列非零元个数即为与某列索引对应的第i个节点处于同一单元中的节点个数第i列的第j个非零元列索引即为与该节点处于同一单元中的节点索引号为向量n
c
前i-1个元素之和；由此形成节点关联向量n
c
、i
c
和s
c
，最终n
c
为1
×
nnum的向量，第i列表示与第i个节点处于同一单元中的节点个数为i
c
为一个维度随具体问题变化的向量，它的大小为1
×
n
sum
，n
sum
为与某节点位于同一单元内的节点个数之和；s
c
为1
×
(nnum+1)的向量，第i列表示与前i-1个节点位于同一单元内的节点个数的总和。
[0022]
本发明进一步的改进在于，步骤4)的具体实现方法如下：
[0023]
对于结构问题，首先需要形成总体偏导矩阵d’e
，对某局部坐标(ξηζ)来说，局部偏导具有如下形式：
[0024][0025][0026]
从而雅各比矩阵j能够通过j＝ud
e
得到，其中u为一维度为3
×
8的节点位移矩阵，提前利用并行循环求得；形函数对总体坐标偏导为，d’e
表示成如下形式：
[0027]
d’e
＝[d’e1
,d’e2
,
……
,d’e7
,d’e8
]
τ
[0028][0029]
本发明进一步的改进在于，步骤5)的具体实现方法如下：
[0030]
组装单元矩阵，对于结构问题，采用8节点六面体单元时需要对各单元形成维度为24
×
24的刚度矩阵以及质量矩阵；基于步骤(4)中的d’e
，将原维度为6
×
24的应力应变矩阵a以及维度为6
×
6的应变位移矩阵d重新分块为：
[0031][0032][0033][0034][0035]
由此得到结构分析问题中某单元内第i个节点对第j个节点的刚度耦合矩阵计算结果为其中nip为高斯积分节点总数；对于单元质量矩阵采用类似的方法，得到以下某单元内第i个节点对第j个节点的质量矩阵计算结果
[0036]
本发明进一步的改进在于，步骤6)的具体实现方法如下：
[0037]
基于步骤1)中得到的位置向量n
ind
,遍历每个单元的节点形成总体刚度矩阵；采用块压缩行格式，对于对应节点索引相同的将它们在总体刚度矩阵相同的位置相加；在此基础上，按行优先的原则，遍历各个3
×
3的矩阵；对于稀疏总体刚度矩阵，存储的只是每个非零的3
×
3矩阵；因此最终得到的是一个维度为的三维数组；
[0038]
与刚度矩阵类似，对于质量矩阵采用压缩行的存储格式；具体为：遍历每个单元非零元素，按位置向量n
ind
得到它们在总体质量矩阵中的位置；对每个非零元只为其对角元素，并且对角元值相同，最终的取值向量v
m
为一维度为3s
c
(nnum+1)
×
1的向量，存有
各个m
ij
的值。
[0039]
本发明进一步的改进在于，步骤7)的具体实现方法如下：
[0040]
具体操作为：位置向量i
c
第个元素对应于列索引c
ind
的第个元素向量s
c
的第i个元素对应于行指针r
p
的第i个元素；对于质量矩阵，列索引c’ind
存有各个在取值向量v
m
中的位置，行指针r’p
保存有每行首个非零元在列索引c’ind
中的位置。
[0041]
本发明进一步的改进在于，步骤9)的具体实现方法如下：
[0042]
对于不同的研究问题，给定不同的约束条件以及载荷条件，并基于并行循环完成约束的添加；对于结构分析，采用位移约束，选出位移约束节点，利用刚度矩阵组装过程中对应的索引方式，将对应节点中对应方向的刚度值通过置大数法并行完成位移约束的添加；对于添加恒定作用力载荷，首先选出恒定力作用节点，这些节点编号已经在网格信息预处理过程中形成一个向量；基于刚度矩阵组装过程中的索引方式，在恒定作用力节点的对应方向上添加作用力，形成恒定作用力向量。
[0043]
与传统的计算方法相比，采用本发明的分析方法具有以下有益效果：
[0044]
本发明算法与框架的高效性与并行性使得大型矩阵组装速度以及荷载向量形成速度显著快于现有开源计算方法，对于大规模的有限元问题能够进行有效求解，矩阵计算速度基本与商业软件ansys持平。在步骤1)中采用了新的算法得到节点相关性信息矩阵，避免了大量循环的使用，较传统算法显著提升了计算效率；步骤2)中通过将矩阵重新分块而不再如传统算法一般直接计算，方便并行计算的同时便于后面组装时采用bsr存储格式存储稀疏矩阵；组装单元矩阵时，对每个单元都是并行独立的计算，通过原子操作防止线程相互干涉，提高了计算速度，缩短了整个流程所耗费的时间；将单元矩阵分块后采用bsr格式存储稀疏总体刚度矩阵。本算法避免了传统算法中大量不必要循环的使用，同时对于必要的矩阵组装循环采用并行、分块的方式进行，从而显著减少了计算时间与内存的占用。
[0045]
综上所述，本发明旨在提供一种基于分块并行计算的有限元分析方法。一方面比现有开源算法更加高效便捷，求解速度与商业软件ansys持平，利用稀疏矩阵算法得到节点相关矩阵，避免了传统算法中大量不必要循环的使用，可以显著缩短数值计算周期并节约试验成本；同时结合新的稀疏矩阵存储特点，重新设计了矩阵的分块计算过程，在减少内存的同时削减了形成总体矩阵所需时间，有利于高效准确的有限元分析
附图说明
[0046]
图1为本发明一种基于分块并行计算的有限元分析方法的流程图；
[0047]
图2为某12扇区透平叶盘的网格剖分示意图；
[0048]
图3为用于得到节点间相关信息的矩阵e结构示意图；
[0049]
图4为用于得到节点间相关信息的矩阵n结构示意图。
具体实施方式
[0050]
下面结合附图对本发明做进一步详细说明。
[0051]
请参阅图1所示，本发明提供的一种基于分块并行计算的有限元分析方法，包括以
下步骤：
[0052]
1)、网格信息预处理。给定某特定研究问题的网格节点单元信息，实现文件信息高效导入。导入矩阵包括节点坐标矩阵以及单元节点信息矩阵，其中节点坐标矩阵保存各节点坐标信息，维度为3
×
nnum，纵轴存储各节点x，y，z轴三个方向的坐标，横轴按顺序排列各个结点；单元节点信息矩阵保存各单元节点编号信息，对8节点六面体单元来说，该矩阵为8
×
enum的矩阵。其中nnum为网格节点个数，enum为网格单元个数。另外额外形成一个节点位置向量n
ind
，该索引向量将节点编号与节点序号联系起来，第i个元素存储了编号为i的节点的索引。
[0053]
2)、形成维度为enum
×
nnum的单元节点信息相关性矩阵e
re
＝[e
1τ
,e
2τ
,
…
,e
enumτ
]
τ
，该矩阵基于导出的单元节点信息矩阵利用稀疏矩阵的coo格式得到。其中nnum为离散后的节点个数，enum为离散后得到的单元个数，e
i
存储有第i个单元中的节点相关性信息。以8节点六面体单元为例，e
i
为一个1
×
nnum的向量，元素e
ij
为1表示第j个节点处在第i个单元中，元素e
ij
为0表示第j个节点与第i个单元无关。
[0054]
3)、形成维度为nnum
×
nnum的n矩阵，该矩阵表达式为n＝e
t
e。当采用压缩行的方式存储该稀疏矩阵时，每一行的首个非零元在列索引中的位置都通过向量存储，某非零元的列代表了与该行节点相关联的节点索引，从而最终得到存储节点相关性信息的向量。当第i个节点与第j个节点位于同一个单元之中时，元素n
ij
非零；当第i个节点与第j个节点位于同一单元之中时，元素n
ij
为零。
[0055]
4)、n矩阵的第i列非零元个数即为与某列索引对应的第i个节点处于同一单元中的节点个数第i列的第j个非零元列索引即为与该节点处于同一单元中的节点索引号为向量n
c
前i-1个元素之和。由此形成节点关联向量n
c
、i
c
和s
c
，最终n
c
为1
×
nnum的向量，第i列表示与第i个节点处于同一单元中的节点个数为i
c
为一个维度随具体问题变化的向量，它的大小为1
×
n
sum
，n
sum
为与某节点位于同一单元内的节点个数之和。s
c
为1
×
(nnum+1)的向量，第i列表示与前i-1个节点位于同一单元内的节点个数的总和。由此得到的节点相关性信息不需要使用大量循环，避免了不必要的时间浪费，因此大幅减少了有限元矩阵计算的时间。
[0056]
5)、单元矩阵预处理。对于结构问题，首先需要形成总体偏导矩阵d’e
。而对某局部坐标(ξηζ)来说，局部偏导很容易得到，具有如下形式：
[0057][0058][0059]
从而雅各比矩阵j能够通过j＝ud
e
得到，其中u为一维度为3
×
8的节点位移矩阵，可以提前利用并行循环求得。形函数对总体坐标偏导为d’e
＝d
e
j-1
。d’e
也可以表示成如下形式：
[0060]
d’e
＝[d’e1
,d’e2
,
……
,d’e7
,d’e8
]
τ
[0061][0062]
6)、组装单元矩阵。对于8节点六面体单元，需要对各单元形成维度为24
×
24的刚度矩阵以及质量矩阵。通过最小二乘法等变分方法求得的单元刚度矩阵具有以下表达式：[k
m
]＝∫∫∫[b][d][b]dxdydz。考虑到工程实际中单元形状较为复杂，因此采用数值积分的形式完成积分的求解。基于步骤(5)中的d’e
，将原维度为6
×
24的应力应变矩阵a以及维度为6
×
6的应变位移矩阵d重新分块为：
[0063][0064][0065][0066][0067]
由此得到结构分析问题中某单元内第i个节点对第j个节点的刚度耦合矩阵计算结果为其中nip为高斯积分节点总数。对每一个单元而言，该过程都利用并行化显著减少计算时间。对于单元质量矩阵采用类似的方法，得到以下某单元内第i个节点对第j个节点的质量矩阵计算结果由于单元中每一个节点对另一个节点的质量矩阵都是对角阵，并且对角元素相等，因此在保存时可以采用只保存一个元素的方式进行。通过分块并行，使得各单元矩阵的组装能够独立进行，并且小矩阵也有利于后面稀疏矩阵的存储，使得保存每个小矩阵元素的位置索引不再必要，因此进一步缩减了矩阵计算时间。
[0068]
7)、基于步骤1)中得到的位置向量n
ind
,遍历每个单元的节点形成总体矩阵。对于结构问题，采用块压缩行格式，对于对应节点索引相同的将它们在总体刚度矩阵相同的位置相加。在此基础上，按行优先的原则，遍历各个3
×
3的矩阵。对于稀疏总体刚度矩
阵，存储的只是每个非零的3
×
3矩阵。因此最终得到的是一个维度为的三维数组。与刚度矩阵类似，对于质量矩阵采用压缩行的存储格式。具体为：遍历每个单元非零元素，按位置向量n
ind
得到它们在总体质量矩阵中的位置。对每个非零元只为其对角元素，并且对角元值相同，最终的取值向量v
m
为一维度为3s
c
(nnum+1)
×
1的向量，存有各个m
ij
的值。
[0069]
8)、形成行指针r
p
和列索引c
ind
。具体操作为：位置向量i
c
第个元素对应于列索引c
ind
的第个元素向量s
c
的第i个元素对应于行指针r
p
的第i个元素对于质量矩阵，列索引c’ind
存有各个在取值向量v
m
中的位置，行指针r’p
保存有每行首个非零元在列索引c’ind
中的位置。通过采用这种存储方法，使得传统方法中保存每一个稀疏矩阵元素值不再必要，能够明显减少矩阵存储的时间消耗与内存占用。
[0070]
9)、添加边界条件。对于不同的研究问题，给定不同的约束条件以及载荷条件，并基于并行循环完成约束的添加。对于结构分析，一般采用位移约束。选出位移约束节点，利用刚度矩阵组装过程中对应的索引方式，将对应节点中对应方向的刚度值取为一个极大的数值比如1e20，即通过置大数法并行完成位移约束的添加。对于添加恒定作用力载荷。首先选出恒定力作用节点，这些节点已经在网格信息预处理过程中形成一个向量。基于刚度矩阵组装过程中的索引方式，在恒定作用力节点的对应方向上添加作用力，形成恒定作用力向量b。
[0071]
10)、问题求解。针对不同的研究问题，比如流场计算、结构分析等采用不同的方法进行求解，得到利用有限元方法求解问题的结果。这些问题包括但不限于静力分析、模态分析等。
[0072]
下面结合具体的实例对本发明一种基于分块并行计算的有限元分析方法进行说明。
[0073]
本实例以某带有12扇区的整圈透平叶片轮盘为例，直板叶片采用整锻式方法成形，顶部有围带连接。透平叶盘在稳定工况下，设置转速为100rad/s，约束轮缘进出汽侧上节点的全部位移。材料参数为弹性模量取2
×
105mpa，泊松比取0.3，密度取7.8
×
10-9
kg/mm3。网格如图2所示，在实际测试中采用了节点数约为5万到200万不等的6个算例。网格节点单元信息矩阵通过读入cdb文件得到。
[0074]
根据给出的各单元编号以及各单元中的8个节点编号的单元信息矩阵，能够建立如图3所示矩阵e，该矩阵维度为enum
×
nnum，当第j个节点位于第i个单元中时，元素e
ij
等于1；当第j个节点不在第i个单元中时，元素e
ij
等于0。注意到单元信息矩阵中已经包含了节点编号，因此只需在每一行将各节点编号所对应的列取值为1即可，避免大量循环的使用。然后得到矩阵n＝e
τ
e，该矩阵结构如图4所示，当第i个节点与第j个节点在同一单元中时，元素n
ij
为一非零元，在图4中用nz表示，其取值与两对应节点出现在同一单元中的次数有关；当第i个节点与第j个节点在同一单元中时，元素n
ij
为0。第i行非零元的个数就是与第i个节点位于同一单元节点的个数，第i行非零元的列索引就是与第i个节点在同一单元中节点的编号。通过将该稀疏矩阵按csr格式存储，能够非常容易地得到这些矩阵。与传统方法相比，该过程避免了大量的循环运算，仅仅通过稀疏矩阵索引关系来得到相应向量，大幅减少了
运算所需时间。
[0075]
在形成相关性矩阵以及向量之后，利用并行执行的for循环遍历每个单元，分块得到各个单元的刚度矩阵与质量矩阵，具体计算方法已经在前文给出。在单元矩阵计算完毕之后，完成总体刚度矩阵和总体质量矩阵的组装。由于在整体矩阵组装时不再需要遍历小矩阵所有元素的位置索引，只需要存储对应的列索引和行指针，相当于将原来需要存储的3
×
3矩阵转变为一个1
×
1的元素，因此减少了矩阵计算的时间以及矩阵的内存占用，有利于高效准确的有限元分析。