CN113609720B - 一种有限元分析的主从自由度处理方法、设备及存储介质 - Google Patents

Abstract

本发明公开了一种有限元分析的主从自由度处理方法、设备及存储介质，包括步骤S1：预先获取有限元模型中各节点的约束关系，根据约束关系生成对应的系数矩阵；步骤S2：将系数矩阵的非零元沿对角线进行排列，并对行列重排后的系数矩阵进行主从自由度选择以输出主自由度和从自由度。本发明采取以高斯消元法为基础的主从自由度选择算法作为核心步骤，使得用户只需要输入约束方程就可以自动选择约束方程的主自由度和从自由度，自动解决由于用户输入主从自由度导致的约束方程出现矛盾的问题，提高运算稳定性；同时本发明将主从自由度选择与矩阵重排相结合，减少运算量与填充元，大幅度提高运算速度。

Description

一种有限元分析的主从自由度处理方法、设备及存储介质

技术领域

本发明涉及有限元分析技术领域，尤其涉及一种有限元分析中的主从自由度处理方法、设备及存储介质。

背景技术

目前，在有限元分析中，遇到接触或连接等问题时，通常使用约束方程来描述某些点的位移与其他点的位移产生某种相关关系。而对于简单的单元连接关系，例如只有两个自由度构成的耦合关系，通过约束方程描述耦合关系后，需要选择主从自由度来分析出自由度之间的依赖关系，若其中一自由度选择为主自由度，则意味着另一个自由度无需参与有限计算，其数值可以直接通过约束方程来获得。但是，将有限元分析应用在汽车制造或土木行业中时，节点单元之间的关系会变得相对复杂，会涉及成千上万个单元彼此耦合的情况，选择主从自由度就变成一项困难的任务。

而现有的主流有限元求解器如Ansys，Abaqus的自动选择功能主要有两种算法：第一种适用于系数矩阵比较简单的情况，此算法将循环系数矩阵所有的行，自动将每一行的第一个非零元选为从自由度，将其他元素选择为主自由度。但是，当系数矩阵较为复杂时，会出现约束方程彼此矛盾的情况；而此算法无法处理这种相互矛盾的约束方程，最终会导致主从自由度选择失败。针对第一种算法的缺点，提出了第二种算法，即通过系数分解的方法选择从自由度，以保证不会出现彼此矛盾的情况；同时对系数矩阵C采用full-pivoting，以保证选择合适的主元，但是该算法的运算速度较慢，使得有限元分析效率无法提高。

发明内容

为了克服现有技术的不足，本发明的目的之一在于提供一种有限元分析的主从自由度处理方法，可自动选择约束方程的主自由度和从自由度，减少出现矛盾的情况，同时提高运算速度，确保分析过程的稳定性。

本发明的目的之二在于提供一种电子设备。

本发明的目的之三在于提供一种存储介质。

本发明的目的之一采用如下技术方案实现：

一种有限元分析的主从自由度处理方法，包括：

步骤S1：获取有限元模型中各节点的约束关系，根据约束关系生成对应的系数矩阵，并对系数矩阵的非零元沿对角线进行排列；

步骤S2：调用行列重排后的系数矩阵进行主从自由度选择以获得主自由度和从自由度，并根据自动选择所获得的主从自由度对有限元模型进行结构仿真以输出有限元仿真结果。

进一步地，预先获取有限元模型中各节点的约束关系的方法为：

根据用户输入或者由连接关系得出的有限元模型中各节点的约束方程；其中，用户输入包括用户直接输入的约束方程或用户输入的指定文件格式的约束数据。

进一步地，对系数矩阵进行重排的方法为：

获取系数矩阵的行数和列数，对系数矩阵的行列数进行比较，将数值最大者所对应的行或列作为第一维度，将较小者所对应的行或列作为第二维度；

循环第一维度，并在循环过程中找到每行或每列非零元素的个数，并通过CuthillMckee稀疏反向排序方式找到第一维度的重排序列；

删除第一维度的重排序列中大于第二维度的行号或列号，得到第二维度的重排序列；

根据第一维度的重排序列和第二维度的重排序列对原系数矩阵进行运算获得重排后的系数矩阵。

进一步地，对系数矩阵进行重排的方法为：

循环系数矩阵中的所有行或列，将每一行或每一列中第一个非零元所在的行号或列号进行标记；

将所有行或所有列按照已标记的行号或列号从大到小的方式进行排列，以将系数矩阵重新排列为反对角排列的矩阵。

进一步地，对重新排列后的系数矩阵进行主从自由度选择的方法为：

循环重排后系数矩阵的每行元素，找到循环过程中第一个非零元素所在列，并在标志向量中将此列标记为从自由度；

实施高斯消元法对从自由度所在列的其他非零元素进行消除；

判断当前行当前列是否有且仅有一个非零元素，若是，则在标志向量中将此行对应的其他非零列标记为主自由度。

进一步地，当判断得出当前行当前列是否有且仅有一个非零元素时，在标志向量中将当前列标记为单点约束进行输出。

进一步地，在进行主从自由度选择时，还包括：预先根据重新排列的系数矩阵输入其对应的右端项系数；在实施高斯消元法对非零元素进行消除时，对右端项系数执行高斯消元操作；判断当前行的非零元个数是否为零，若是，再判断对应的右端项系数是否为零，若为零，则在标志向量中将此行标记为单点约束进行输出；若不为零，则执行用户输入有误的报警提示。

进一步地，在输出主自由度和从自由度之前，循环标记向量，根据标记向量的标记将全部自由度划分为主自由度、从自由度、与退化的单点约束进行输出。

本发明的目的之二采用如下技术方案实现：

一种电子设备，其包括处理器、存储器及存储于所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的有限元分析的主从自由度处理方法。

本发明的目的之三采用如下技术方案实现：

一种存储介质，其上存储有计算机程序，所述计算机程序被执行时实现上述的有限元分析的主从自由度处理方法。

相比现有技术，本发明的有益效果在于：

本发明采取以高斯消元法为基础的主从自由度选择算法作为核心步骤，使得用户只需要输入约束方程就可以自动选择约束方程的主自由度和从自由度，避免出现用户指定主从自由度导致约束方程出现矛盾的情况，提高运算稳定性；同时本发明将主从自由度选择与矩阵重排相结合，减少运算量与填充元，大幅度提高运算速度。

附图说明

图1为本发明有限元分析的主从自由度处理方法的流程示意图；

图2为本发明系数矩阵重排的流程示意图；

图3为本发明主从自由度运算的流程示意图；

图4为本发明主从自由度运算的逻辑示意图；

图5为建筑横梁结构的结构模型图；

图6为建筑横梁结构通过有限元法建模所得的节点关系图；

图7为本发明主从自由度自动选择算法性能测试结果图；

图8为本发明改进后Cuthill-Mckee算法性能测试的结果图；

图9为本发明经过改进后的Cuthill-Mckee重排的算法耗时测试结果图；

图10为未经过改进后的Cuthill-Mckee重排的算法耗时测试结果图；

图11为采用快速重排方法和采用改进后Cuthill-Mckee的测试比对结果图；

图12为采用快速重排方法后进行主从自由度选择运算与采用改进的Cuthill-Mckee重排方法后进行主从自由度选择算法的耗时测试结果图。

具体实施方式

下面，结合附图以及具体实施方式，对本发明做进一步描述，需要说明的是，在不相冲突的前提下，以下描述的各实施例之间或各技术特征之间可以任意组合形成新的实施例。

实施例一

本实施例提供一种有限元分析的主从自由度处理方法，可实现自动选择主从自由度的效果，提高运算稳定性和运算速度。

如图1所示，本实施例的主从自由度处理方法具体包括如下步骤：

步骤S1：预先获取有限元模型中各节点的约束关系，根据约束关系生成对应的系数矩阵，将系数矩阵的非零元沿对角线进行排列；

步骤S2：调用行列重排后的系数矩阵进行主从自由度选择以获得主自由度和从自由度，并根据自动选择所获得的主从自由度对有限元模型进行结构仿真以输出有限元模型的应力应变仿真结果。

本实施例中用户可预先在系统中直接输入约束方程，系统根据用于输入的约束方程自动选择主自由度和从自由度；此外，系统还可提供其他输入形式让用户输入对应的约束数据，根据约束数据自动选择主从自由度。例如：用户可输入BDF文件格式的数据来表达各单元节点之间的约束关系。若用户输入BDF文件格式的约束数据，系统可根据用户输入的数据转化为对应的约束方程，再根据约束方程转化为对应的系数矩阵；在此，不对用户的输入方式进行限定。

本实施例为了减少主从自由度选择运算过程中的运算量，需对系统矩阵进行重排操作。而本实施例的重新排列原则为将矩阵的非零元沿着矩阵正对角线或反对角线进行排列；如图2所示，本实施例将矩阵的非零元尽量沿着反对角线进行排列的方法为：

步骤a：将系数矩阵输入系统后，系统获取系数矩阵的行数和列数，对系数矩阵的行列数进行比较，将较大者所对应的行或列作为第一维度，将较小者所对应的行或列作为第二维度；即若该系数矩阵的行数大于列数，则将系数矩阵的行数作为第一维度，将列数作为第二维度；若该系数矩阵的列数大于行数，则将列数作为第一维度，将行数作为第二维度。

步骤b：循环第一维度，并在循环过程中找到每行或每列的degree，degree等于该行或该列的非零元素的个数，并对非零元素进行Cuthill-Mckee稀疏反向排序操作以找到第一维度的重排序列T₁；而所述稀疏反向排序方法为现有技术，在此不对其进行详细描述。

步骤c：在第一维度的重排序列中，删除第一维度的重排序列中大于第二维度的行号或列号，得到第二维度的重排序列T₂；

步骤d：通过公式C′＝T₁*C*T₂将第一维度的重排序列、第二维度的重排序列与原系数矩阵进行运算获得重排后的系数矩阵C′，并将重排后的系数矩阵C′输出以进行主从自由度选择运算。

本实施例在矩阵重排过程中，对现有的Cuthill-Mckee矩阵重排算法进行改进，使得原Cuthill-Mckee矩阵重排算法不再局限于处理对称的方阵，也可处理原Cuthill-Mckee矩阵重排算法无法处理的非对称的矩形阵。本实施例将矩阵的非零元尽量沿着反对角线进行排列，可以在重新排列非零元的同时，减少矩阵的带宽，大幅度减少运算量，提高运算速度。

本实施例对系统数据进行重新排列后，通过主从自由度节点选择算法对主从自由度进行自动选择，如图3、图4所示，具体方法为：

步骤A：输入重排后的系数矩阵，循环系数矩阵的每行元素，查找循环过程中第一个非零元素所在列的位置，选择第一个非零元素作为从自由度，并在标志向量中将此列标记为从自由度；

步骤B：判断从自由度所在列是否有其他非零元素，若有，则实施高斯消元法对从自由度所在列的其他非零元素进行消除；

步骤C：判断当前行当前列是否有且仅有一个非零元素，如果当前行当前列有且只有一个非零元素，则在标志向量中将此行对应的其他非零列标记为主自由度，从而完成主从自由度的自动选择功能。

此外，本实施例自动选择主自由度和从自由度的同时，还可自动输出单点约束。即在步骤C中当判断得出当前行当前列是否有且仅有一个非零元素时，在标志向量中将当前列标记为单点约束。

步骤D：循环标记向量，根据标记向量的标记输出主自由度、从自由度和单自由度，从而将全部自由度自动划分为主自由度、从自由度和单点约束，实现自动选择主从自由度的目的。

本实施例的主从自由度选择运算过程中，在输入系数矩阵的同时，还需确定该系数矩阵对应的右端项，其右端项是根据系数矩阵所对应的约束方程来确定的，右端项可体现约束方程中两个自由度之间的位移差；本实施例通过其右端项的值来确定实施用户输入的约束方程是否出现矛盾。图4为本实施例主从自由度选择运算的逻辑图，如图4所示，在输入系统矩阵的同时输入系数矩阵的右端项，在循环系数矩阵的每行元素过程中，对从自由度所在列的其他非零元素通过高斯消元法进行消除后，将右端项执行同样的操作；此时再判断当前行的非零元个数是否为零，若是，再判断对应的右端项系数是否为零，若为零，则在标志向量中将此行标记为单自由度进行输出；若不为零，代表着用户输入的约束方程出现矛盾，此时执行用户输入有误的报警提示操作。

为保证算法的稳定性，本实施例采取以高斯消元法为基础的主从自由度选择算法作为核心程序，为保证算法的速度，在实施高斯消元之前，对系数矩阵进行必要的重排，以减少运算量与填充元；同时，本实施例在矩阵重排时，对现有的Cuthill-Mckee矩阵重排算法进行一系列的修改，使之可以处理系数矩阵这类矩形阵；本实施例采用主从自由度选择与矩阵重排相结合，大幅度提高运算效率和稳定性。

考虑土木行业的案例，在土木行业中，有很多桁架结构，在钢结构顶点的位置，会采用各种方式，将左右横梁连接在一起。以具体一横梁结构作为例子对本实施例的运算处理方法进行说明：

如图5所示，图5中的横梁结构包括一号单元、二号单元和三号单元，中间三号单元在横向的变形梁较小，因而可以简化为刚体单元，使用有限元法将上面的模型建模如图6所示，将上述模型抽象为三个单元，通过约束方程组合在一起。其节点编号如图6所示，其中1节点到2节点为梁单元，2节点到3节点为刚体单元，10节点到4节点为梁单元。3节点与10节点采用铰接。

用户输入如下BDF文件格式的约束数据来表示节点间的约束关系：

表1用户输入的BDF文件格式的约束数据表

MPC编号	节点编号	自由度编号	系数	节点编号	自由度编号	系数
							1	3	1	1.0	10	1	1.0
1	3	2	1.0	10	2	1.0
							1	3	3	1.0	10	3	1.0
1	2	1	1.0	3	1	1.0
							1	2	2	1.0	3	2	1.0
1	2	3	1.0	3	3	1.0

由于梁单元每个节点有六个自由度，因此将上述单元节点自由度编号如下：

表2节点编号与自由度编号关系表(ID array)

Node Number	u<sub>1</sub>	u<sub>2</sub>	u<sub>3</sub>	u<sub>4</sub>	u<sub>5</sub>	u<sub>6</sub>
							1	1	2	3	4	5	6
2	7	8	9	10	11	12
							3	13	14	15	16	17	18
4	19	20	21	22	23	24
							10	25	26	27	28	29	30

根据上述自由度编号，可以将约束方程写为：

u₁₃-u₂₅＝0；

u₁₄-u₂₆＝0；

u₁₅-u₂₇＝0；

u₇-u₁₃＝0；

u₈-u₁₄＝0；

u₉-u₁₅＝0；

则系数矩阵可以写为：

经过cuthill-mckee重排后：

可以看到，系数矩阵的带宽被明显减少了。

再经过主从自由度选择，主自由度为：{u₇，u₈，u₉}，从自由度为：{u₁₃，u₁₄，u₁₅，u₂₅，u₂₆，u₂₇}，相应的约束方程被转换为：

u₁₃-u₇＝0；

u₁₄-u₈＝0；

u₁₅-u₉＝0；

u₂₅-u₇＝0；

u₂₆-u₈＝0；

u₂₇-u₉＝0；

本实施例利用具有矩阵重排的主从自由度处理方法自动选取合适的主自由度与从自由度，相当于重新确定目标物有限元模型中各单元节点之间的约束关系，并根据自动选择的主从自由度转换为约束方程的形式输入至ZwSim结构仿真软件中进行结构仿真以输出应力应变等仿真结果；经过实验可知，经本实施例主从自由度选择后仿真结果与现有NxNastran软件的计算结果相一致。

本实施例采用如下系数矩阵来测试以上算法的性能：

该系数矩阵的右端项g＝0。

通过选择不同的n，来测试该算法的性能。经过测试，主从自由度选择算法和本实施例改进后的Cuthill-Mckee重排算法的性能如下图7、图8所示，图7、图8中横轴代表不同大小的矩阵C，纵轴代表所消耗的时间，单位为s；经过测试可知，本实施例中主从自由度选择与矩阵重排算法的时间复杂度都为线性，本实施例采用具有矩阵重排的主从自由度自动选择算法在运算速度上具有较大的优势。此外，再对运算速度进行测试，经过测试可知，经过矩阵重排和未经矩阵重排算法的系统运行速度存在明显的差异，具体如下：

如图9和图10所示，图9和图10均为经过本实施例改进后的Cuthill-Mckee重排与未经Cuthill-Mckee重排的算法耗时对比图，图9以折线图方式来呈现二者的耗时测试结果，图10以列表方式来呈现二者的耗时测试结果，从图9和图10可见，经过矩阵重排后的算法耗时迅速下降，总耗时比未经重排降低几百倍，时间复杂度有显著降低，由此可以说明，矩阵重排在该流程中的重要作用，大幅度提升运算速率。

实施例二

本实施例在实施例一的基础上对系数矩阵的重排方法进行改进，即本实施例对系数矩阵进行重新排列的方法还可以采用快速重排算法来实现，而快速重排算法包括行重排和列重列两种方式，其中快速行重排算法包括如下步骤：

第一步，循环所有的行，找出每一行第一个非零元所在的列号，记为MC(MinimumColumn)；

第二步，将所有行按照MC由大到小排列，即MC越大的行，越靠上，MC越小的行越靠下。同时，对应的右边项也要按同样的顺序进行重排。这样，一个任意的矩阵，便按照反对角的趋势完成了重新排列。

以上的快速重排算法针对的是行重排，还有一种快速重排算法针对列重排；列的快速重排实现方法与行的快速重排基本一致，效率和时间复杂度也基本一致，如下是快速列重排的步骤：

第一步，循环所有的列，找到每一列第一个非零元对应的行号，记为MR(MinimumRow)；

第二步，将所有列按照MR由大到小排列，即MR越大的列，越靠近矩阵的左边，MR越小的列越靠近右边；

第三步，与行的快速重排相比，列的快速重排不需要对右边项进行重排，但是需要将原来的自由度替换为新的重拍后的自由度。

与上述对现有的Cuthill-Mckee矩阵重排算法进行改进的重排方法相比，快速重排的优点是重排速度较快，也可以处理非方阵问题，缺点是无法减少矩阵的带宽，因此在系数矩阵带宽较小的情况比较有优势。

本实施例由于主从自由度算法的特质，对于反对角排列的矩阵，运算量最小，而对角排列的矩阵，运算量最大。因此，快速重排正是利用了这一特点，尽可能将矩阵按照反对角重新排列，以减少主从自由度选择的运算量。

此外，主从自由度算法将循环顺序将从左往右改为从右往左，可以达到同样的效果。这时，将系数矩阵的非零元沿主对角排列，也可以大幅减少计算量。而以上三种重排算法在不改变原理的情况下，可以稍加改变，即可将系数矩阵沿反对角重排变为沿对角线重排。举个例子：

1)对Cuthill-Mckee进行改进的重排方案中，可以将第一维或第二维的重排序列进行反序排列，再执行同样的操作，即可将重排方向由反对角变为对角。

2)快速行重排或者快速列重排也可以将之前形成的MR序列或者MC序列反向排列，这样重排之后的非零元也会沿着对角线排列。

综上，所有按照实施例一或实施例二所记载的重排方法延伸而成的沿主对角线排列的算法也属于本专利保护范围。

本实施例对快速重排方法的运算耗时进行测试，其测试结果如图11所示，与实施例一种改进的Cuthill-Mckee重排，本实施例快速重排的优点是重排速度较快，也可以处理非方阵问题；但其缺点是无法减少矩阵的带宽，因此在系数矩阵带宽较小的情况比较有优势。同时，本实施例还对快速重排方法后进行主从自由度选择进行耗时测试，如图12可知，本实施例采用快速重排方法后进行主从自由度选择运算的耗时与实施例一中采用改进的Cuthill-Mckee重排方法后进行主从自由度选择算法的耗时相当，同样可实现加快计算速度的目的。

实施例三

本实施例提供一种电子设备，其包括处理器、存储器及存储于所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现实施例一中的有限元分析的主从自由度处理方法；另外，本实施例还提供一种存储介质，其上存储有计算机程序，所述计算机程序被执行时实现上述的有限元分析的主从自由度处理方法。

本实施例中的设备及存储介质与前述实施例中的方法是基于同一发明构思下的两个方面，在前面已经对方法实施过程作了详细的描述，所以本领域技术人员可根据前述描述清楚地了解本实施中的系统的结构及实施过程，为了说明书的简洁，在此就不再赘述。

上述实施方式仅为本发明的优选实施方式，不能以此来限定本发明保护的范围，本领域的技术人员在本发明的基础上所做的任何非实质性的变化及替换均属于本发明所要求保护的范围。

Claims (8)

1.一种有限元分析的主从自由度处理方法，其特征在于，包括：

步骤S1：获取有限元模型，所述有限元模型包括横梁结构模型，横梁结构模型中具有若干个单元节点；并获取模型中各单元节点间耦合关系，根据耦合关系确定各单元节点的约束方程，并根据约束方程生成对应的系数矩阵，并对系数矩阵的非零元沿对角线进行排列；对系数矩阵进行重排的方法包括：

获取系数矩阵的行数和列数，对系数矩阵的行列数进行比较，将数值最大者所对应的行或列作为第一维度，将较小者所对应的行或列作为第二维度；

循环第一维度，并在循环过程中找到每行或每列非零元素的个数，并通过CuthillMckee稀疏反向排序方式找到第一维度的重排序列；

删除第一维度的重排序列中大于第二维度的行号或列号，得到第二维度的重排序列；

根据第一维度的重排序列和第二维度的重排序列对原系数矩阵进行运算获得重排后的系数矩阵；

步骤S2：调用行列重排后的系数矩阵进行主从自由度选择以获得主自由度和从自由度，其方法为：循环重排后系数矩阵的每行元素，找到循环过程中第一个非零元素所在列，并在标志向量中将此列标记为从自由度；

实施高斯消元法对从自由度所在列的其他非零元素进行消除；

判断当前行当前列是否有且仅有一个非零元素，若是，则在标志向量中将此行对应的其他非零列标记为主自由度；

步骤S3：将自动选择的主从自由度转换为约束方程，并将其输入至结构仿真软件中进行结构仿真以输出模型的应力应变仿真结果。

2.根据权利要求1所述的有限元分析的主从自由度处理方法，其特征在于，预先获取有限元模型中各节点的约束关系的方法为：

根据用户输入或者由连接关系得出的有限元模型中各节点的约束方程；其中，用户输入包括用户直接输入的约束方程或用户输入的指定文件格式的约束数据。

3.根据权利要求1所述的有限元分析的主从自由度处理方法，其特征在于，对系数矩阵进行重排的方法为：

循环系数矩阵中的所有行或列，将每一行或每一列中第一个非零元所在的行号或列号进行标记；

将所有行或所有列按照已标记的行号或列号从大到小的方式进行排列，以将系数矩阵重新排列为反对角排列的矩阵。

4.根据权利要求1所述的有限元分析的主从自由度处理方法，其特征在于，当判断得出当前行当前列是否有且仅有一个非零元素时，在标志向量中将当前列标记为单点约束进行输出。

5.根据权利要求1所述的有限元分析的主从自由度处理方法，其特征在于，在进行主从自由度选择时，还包括：预先根据重新排列的系数矩阵输入其对应的右端项系数；在实施高斯消元法对非零元素进行消除时，对右端项系数执行高斯消元操作；判断当前行的非零元个数是否为零，若是，再判断对应的右端项系数是否为零，若为零，则在标志向量中将此行标记为单点约束进行输出；若不为零，则执行用户输入有误的报警提示。

6.根据权利要求1所述的有限元分析的主从自由度处理方法，其特征在于，在输出主自由度和从自由度之前，循环标记向量，根据标记向量的标记将全部自由度划分为主自由度、从自由度、与退化的单点约束进行输出。

7.一种电子设备，其特征在于，其包括处理器、存储器及存储于所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现权利要求1～6任一所述的有限元分析的主从自由度处理方法。

8.一种存储介质，其特征在于，其上存储有计算机程序，所述计算机程序被执行时实现权利要求1～6任一所述的有限元分析的主从自由度处理方法。

Images