专利名称：金属体积塑性成形中有限元分析刚度矩阵存储与生成方法
技术领域：
本发明涉及金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，适合于锻造、挤压、轧制、旋转锻造等体积成形过程的有限元数值模拟。
背景技术：
在工业生产中，锻造、挤压、轧制等金属体积塑性成形制造技术支撑国民经济发展与国防建设的主要技术之一，在汽车、航天航空、装备制造、兵器、能源、造船等行业具有广泛用途。体积塑性成形过程就是将室温或高温下(冷、温、热成形)的金属材料(坯料) 在锻压设备及其安装其上的模具作用下发生塑性变形，并成形为所需形状和性能的零件的过程，体积塑性成形材料一般被假设为刚(粘)塑性材料，因此，坯料在模具型腔内发生大而剧烈的塑性变形，存在材料非线性、几何非线性和复杂边界条件，成形过程的求解十分复杂。成形工艺和模具设计对于成形件的质量及其能否顺利成形具有重要影响，也是影响节能降耗和精确成形的关键。而掌握金属材料在模具型腔内的塑性变形行为和流动规律是设计成形工艺和模具的关键理论依据。有限元方法作为工程问题数值分析的重要方法，已在工程领域发挥了重要作用，金属体积塑性成形过程有限元数值模拟方法是研究金属材料塑性变形行为和流动规律的主要方法，可获得金属材料流动规律和各种场变量分布结果，验证工艺设计、参数选择和模具设计方案的可行性，进而实现成形工艺和模具的优化优化，提高产品开发的速度、质量和降低开发成本。对于金属体积塑性成形问题的有限元分析，首先通过塑性变形控制方程，建立金属体积塑性变形力学问题的数学模型，将几何模型(模具与变性材料)进行有限元法离散， 建立有限元刚度方程，实现成形过程的数值模拟，有限元等数值分析方法已经在工程领域， 特别是在金属成形制造和结构分析等领域发挥了重要作用。有限元数值求解金属体积塑性成形问题的步骤为建立问题的CAD几何模型(变形材料的几何形状)，将CAD几何模型进行有限单元离散，得到几何模型的有限元网格模型，该网格模型由有限数量的网格单元组成；首先建立网格模型中每个单元的刚度方程， 将所有单元的刚度方程进行组装，形成关于网格模型所有节点速度的非线性总体刚度方程组；通过对关于节点速度的非线性刚度方程的线性化处理，获得关于节点速度增量的线性刚度方程组；采用Newton-Raphson迭代法对速度场进行迭代求解。由于对节点速度增量进行求解，因此，对于一个金属体积塑性成形过程，需要将成形时间分解为上百甚至几百个时间步(长)，在每个时间步长内均需要计算单元刚度矩阵和节点力向量，组装成节点速度增量的总体刚度方程，施加速度边界条件，求解节点速度增量的总体刚度方程(大型线性方程组)，得到速度修正量，再用速度修正量对速度场进行修正，然后重复上述步骤，对速度场进行迭代求解，直到速度场收敛，一般情况下，至少需要7至10次迭代方可获得收敛的速度场；根据获得的收敛速度场，刷新变形工件形状和模具位置，才能完成一个时间步长内的数值计算。然后再进行下一个步长内的数值计算，直至达到所要求的变形程度。由此可见，对于金属体积塑性成形过程的数值模拟，需要频繁的求解总体刚度方程(大型线性方程组)。例如对于中等复杂的三维金属体积塑性成形问题，一般需要200次以上的模拟步数，每个模拟步内大致需要迭代9次速度场，而每次速度场迭代就需要求解一次大型线性方程组，如果一个金属体积塑性成形问题的有限元节点数目为5万个(中等复杂的三维问题)，就需要求解由15万个方程构成的大型线性方组1800次。如果分析更复杂问题，其计算工作量更加巨大。近年来，随着重大装备、船舶、汽车和风电等大型工程的需要，对大型、重型金属体积塑性成形零件的需求日益增多，需要分析的成形问题日益复杂，几何尺寸大，有限元模型的节点和单元数量多，数据量庞大，有限元计算效率低下，所需计算机内存剧增，甚至在分析大型和超大型问题时，经常出现存储空间不足的问题，这大大限制了有限元数值方法在金属成形过程分析中的应用。不仅金属体积塑性成形问题，而且很多其它工程问题的有限元分析最终都会归结为求解大型线性方程组，有限元的求解效率主要取决于线性方程组刚度矩阵的存储、生成及其刚度方程的求解效率。金属体积塑性成形问题有限元刚度矩阵为大型、带状稀疏矩阵，矩阵的非零元素一般集中在对角线区域，针对这种情况和为了提高计算效率，人们提出了许多节约存储空间的存储方法，例如，半带宽存储方法、变带宽存储方法、稀疏方程组预处理迭代求解方法等，尽管这些方法减少了刚度矩阵的存储空间，但也无可避免地存储了部分非零元素，浪费了部分储存空间，同时这些方法均依赖于节点编号，需要对节点编号顺序进行优化，从而也浪费了计算时间。现有的一维变带宽压缩存储方法采用一维数组按行的顺序保存刚度矩阵的非零元素，采用两个辅助数组记录列号和每列起始元素在刚度矩阵数组的位置，节约了存储空间，避免了节点编号优化问题，但所采用的辅助数组与刚度矩阵存储数组具有同样长度，仍存在消耗存储空间和浪费计算时间的问题。因此，采用有限元分析金属体积塑性成形等工程问题时，有限元刚度矩阵元素的高效压缩存储与生成方法成为迫切需要解决的问题。
发明内容
本发明的目的是为金属体积塑性成形问题有限元数值模拟过程中存在的刚度矩阵数据存储空间大和刚度方程计算效率不高的问题，提供一种新的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，它不仅具有显著节约刚度矩阵存储空间和提高计算效率的特点，而且具有随金属体积塑性成形问题有限元网格模型节点与单元数量的增加，其节约空间越显著和求解效率越高的特点。本发明是通过下面的技术方案来实现的一种金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，包括以下步骤(1. 1)首先采用现有CAD设计软件或者三坐标测量仪等反求设备获得金属体积塑性成形问题的三维几何模型，采用现有网格生成方法生成成形工件的有限元网格模型，其次根据塑性成形工件的有限元网格模型拓扑结构，确定其网格节点的“广义相邻节点对”， 建立“广义相邻节点对”与工件有限元网格模型总体刚度矩阵非零子矩阵之间的对应关系；(1. 2)基于“广义相邻节点对”，储存非零子矩阵，开辟一个非零子矩阵列号存储数组J即整型索引数组，记录与工件有限元网格模型中“广义相邻节点对”对应的总体刚度矩阵非零子矩阵的列号，开辟另一个整型索引数组K，记录工件有限元网格模型刚度矩阵中每行首个非零子矩阵在数组J中的位置，确定工件有限元网格模型刚度矩阵中非零子矩阵的位置与数量，根据非零子矩阵与非零元素的对应关系，确定工件有限元网格模型刚度矩阵非零元素的位置和数量；(1. 3)计算工件有限元网格模型中每个单元刚度矩阵元素在工件有限元网格模型总体刚度矩阵中所在的非零子矩阵，利用索引数组J和K中的数据，确定单元刚度矩阵元素在网格模型总体刚度矩阵存储数组A中的位置(具体确定步骤见4. 1 4. 9)；按照网格模型总体刚度矩阵的生成规则，将网格模型的每个单元刚度矩阵元素对号入座地组装到网格模型的总体刚度矩阵存储数组A中；遍历网格模型中的所有单元后，最终生成金属体积成形问题网格模型的总体刚度矩阵存储数组；(1. 4)采用现有的对称正定系数矩阵松弛预处理共额梯度法，对网格模型的总体刚度方程(大型线性方程组)进行求解，获得金属体积成形过程网格模型中所有节点的速度增量场的数值解。所述步骤(1. 1)中，确定网格节点的“广义相邻节点对”的方法为(2. 1)在体积成形问题的有限元网格模型中，每个单元内的任意两个节点构成的节点对(包括节点自身与自身构成的节点对)，称为“广义相邻节点对”；将相邻节点有顺序的节点对称为“有序广义相邻节点对”；反之，称为“无序广义相邻节点对”，若网格模型中某个单元两个节点的全局节点编号为(I，J)，则构成两个“有序广义相邻节点对”^^P^j， 若不考虑节点的顺序，则构成一个“无序广义相邻节点对” 即认为^7和云 是同一个节点对，记为^或^工。为了应用和表示的方便，可采取下标由小到大的形式来表示“无序广义相邻节点对”，即统一记为BuG ^ J)；(2. 2)在网格模型中，每个“有序广义相邻节点对”对应网格模型的总体刚度矩阵中的一个与所分析问题维数相同的非零子矩阵，非零子矩阵的位置由“有序广义相邻节点对”的全局节点序号决定，例如，在网格模型中的“有序广义相邻节点对” ^7对应网格模型的总体刚度矩阵中第I行、第J列的非零子矩阵；(2. 3)在网格模型中，每个“无序广义相邻节点对”对应网格模型的总体刚度矩阵的上三角矩阵中的一个与所分析问题维数相同的非零子矩阵，非零子矩阵的位置由 “无序广义相邻节点对”的全局节点序号决定；例如在网格模型中的“无序广义相邻节点对”^ J)，对应网格模型的总体刚度矩阵的上三角矩阵中第I行、第J列的非零子矩阵；(2. 4)由于金属体积塑性成形过程有限元分析的总体刚度矩阵为对称矩阵，为了节约空间可以只存储其上三角矩阵，因此本发明中提到的网格模型的总体刚度矩阵存储是指其上三角矩阵的存储，由于网格模型中“无序广义相邻节点对”与网格模型的总体刚度矩阵的上三角矩阵中的非零子矩阵一一对应，并且节点对的全局节点编号与非零子矩阵的行号、列号对应，因此通过这种对应关系，由网格模型中所有的“无序广义相邻节点对”的数量和全局节点号即可确定网格模型刚度矩阵的上三角矩阵中非零子矩阵的数量和位置(行号与列号)，并可进一步确定组成非零子矩阵的非零元素的数量和位置。所述步骤(1. 2)中，还包括以下步骤
(3. 1)采用一个浮点型数组A按行顺序存储网格模型总体刚度矩阵上三角矩阵中的各非零元素；(3. 2)采用整形数组J记录与“无序广义相邻节点对”对应总体刚度矩阵上三角矩阵中的非零子矩阵的列号；(3. 3)采用整形数组K记录总体刚度矩阵上三角矩阵中每行第一个非零子矩阵在数组J中的编号，并以此确定非零子矩阵在总体刚度矩阵上三角矩阵的行号，因为K[i]为总体刚度矩阵上三角矩阵中第i行第一个非零子矩阵在数组J中的编号，而K[i+1]为总体刚度矩阵上三角矩阵中第i+1行第一个非零子矩阵在数组J中的编号，由此可知在数组J 中编号从K[i]到K[i+1]-1对应的所有非零子矩都在第i行。所分析问题的几何模型划分后的有限元网格模型为三维的六面体单元、四面体单元或其二者的混合单元模型；或者所分析问题的几何模型划分后的有限元网格模型为二维的四边形单元、三角形单元或其二者的混合单元模型，都可以采用本发明的方法进行网格模型总体刚度矩阵的存储。所述步骤(1.3)中，对于网格模型中的一个单元，若单元的节点数为η、物理量的维数为m，记该单元η个节点的整体编号为I[i]，i = 1，2，···，η，该单元刚度矩阵为U(mXn， mXn)，则该单元刚度矩阵的任意元素为Uijkl = U[(i-1) Xm+k] [(j_l) Xm+1]其中，i，j = 1,2,…，η为该单元节点的局部编号，k，l = 1,2,…，m为物理量维
数编号；确定Uljkl在网格模型总体刚度矩阵存储数组A中序号的过程如下(4. 1)确定Uijkl的整体节点编号I = I[i]，J = I[j]，其中I、J为UijklK在的非零子矩阵Au中的在网格模型总体刚度矩阵上三角矩阵中的行号和列号；非零子矩阵的位置如图1所示。(4. 2)确定非零子矩阵A11之前非零子矩阵的数量=N1 = K[I]-1 ；(4. 3)确定非零子矩阵A11之前非零元素的数量NN1 = N1Xm2-(I-I) XmX (m-l)/2其中，m2为每个非零子矩阵非零元素的数量，m(m-l)/2为每行的对角非零子矩阵引起的非零元素的减少项；如图2中的非零子矩阵A11为m = 3的三维子矩阵，只存储了 6 个非零元素，与其它子矩阵相比少存3个。(4. 4)确定第I行非零子矩阵中第k行之前非零元素的数量NNik= (K[I+1]-K[I]) XmX (k-1)(4.5)确定第I行中，第J列之前的非零子矩阵的数量。搜索数组J从标号K[I] 到κ[1+1]中J[II] = J的数组标号II，第I行中，第J列之前的非零子矩阵数量为Nt = II-K[I]；(4. 6)确定Uijkl在总体刚度矩阵中所在行中的编号NNjl = SjXm-TT+1其中，TT为由于对角非零子矩阵引起的非零元素的减少项，TT = k(k-l)/2 ；(4. 7)确定Uijkl在网格模型总体刚度矩阵储存数组A中的编号NNukl = NNj+NN^+NNj!
(4. 8)根据网格模型总体刚度矩阵的生成规则和Uijkl在网格模型总体刚度矩阵储存数组A中的位置编号NNukl，将Uijkl组装或累加到网格模型的总体刚度矩阵储存数组A 中；(4. 9)重复步骤(4. 1) (4. 8)，遍历网格模型中的所有单元刚度矩阵中的每一个元素，即可获得网格模型的总体刚度矩阵一维存储数组A。求解网格模型所采用的数值分析方法最终归结为稀疏、对称总体刚度矩阵的存储、生成及其刚度方程(即线性方程组)求解，所分析的网格模型的总体刚度方程求解的参
量为全量或者增量。本发明的有益效果是根据本发明提出的”广义相邻节点对”及其与网格模型总体刚度矩阵非零子矩阵之间的对应关系，存储网格模型总体刚度矩阵中的非零子矩阵及其位置，改变了以往存储网格模型总体刚度矩阵中非零元素及其位置的方法；本发明的压缩存储方法所需的存储空间随网格模型节点数量的增加呈线性函数增加，不同于现有压缩存储方法所需的存储空间随网格模型节点数量的增加呈二次函数增加，因此，本发明的方法不仅显著节约了网格模型总体刚度矩阵的存储空间和提高了运算效率，而且随网格模型节点和单元数量的增加，其节约存储空间的能力越显著；基于存储非零子矩阵的网格模型总体刚度矩阵方程的松弛预处理共额梯度法，无需对对称逐步超松弛迭代的分裂矩阵进行求逆和另外开辟数组存放分裂矩阵，可直接对分裂矩阵的非零元素进行运算，求解和运算效率高；因此，特别适合于金属体积成形这种需要反复进行增量求解和迭代求解的大型复杂工程问题的数值分析。
图1为总体刚度矩阵中非零子矩阵的分布示意图；图2为总体刚度矩阵中非零元素的位置示意图；图3为金属刚(粘)塑性成形网格模型总体刚度矩阵方程压缩存储和生成算法流程图；图4为由一个八节点六面体单元构成的网格模型；图5为由两个八节点六面体单元构成的网格模型；图6a为锻造过程有限元数值模拟；图6b为锻造过程有限元数值模拟；图6c为锻造过程有限元数值模拟；图6d为锻造过程有限元数值模拟。
具体实施例方式下面结合附图与实施例，对本发明做进一步说明。图3为金属刚(粘)塑性成形网格模型总体刚度矩阵方程压缩存储和生成算法流程图。根据图3所示，金属刚(粘)塑性成形网格模型总体刚度矩阵方程压缩存储和生成算法流程如下根据已经生成的金属塑性成形问题的几何模型和有限元网格模型，确定整个网格模型的“广义相邻节点对”并计算其数量；开辟网格模型总体刚度矩阵存储数组A，开辟非零子矩阵列号存储数组J，开辟每行首个非零子矩阵位置存储数组K ；对网格模型中的所有节点进行循环，对于任一节点i，搜索网格模型中包含节点i的所有单元；搜索包含节点i的所有单元所包含的所有节点，确定节点i的“广义相邻节点对”;依据包含节点i “广义相邻节点对”的编号，生成索引数组J ；依据包含节点i “广义相邻节点对”的数量，生成索引数组K ；判断网格模型中所有节点的循环是否结束，如果还未结束，则对下一个节点重复上述步骤，直至对网格模型中的所有节点处理完毕。对网格模型中的所有单元进行循环，生成每个单元的刚度矩阵U，判断单元刚度矩阵每个元素在网格模型总体刚度矩阵上三角矩阵中所在的非零子矩阵的位置，计算所在非零子矩阵所在行之前的非零子矩阵及非零元素的数量；计算单元刚度矩阵元素对应的网格模型总体刚度矩阵上三角矩阵的元素所在行之前的所有非零元素的数量；计算单元刚度矩阵元素对应的网格模型总体刚度矩阵上三角矩阵元素所在行的顺序编号；计算单元刚度矩阵元素对应的网格模型总体刚度矩阵上三角矩阵元素在网格模型总体刚度矩阵存储数组A中的位置编号，依据网格模型总体刚度矩阵的生成规则，将该单元刚度矩阵元素装入网格模型总体刚度矩阵上三角矩阵存储数组。判断该单元矩阵中所有元素是否计算完毕，如果尚未完毕，则对该单元中的下一个元素进行上述处理，直至对该单元中所有元素处理完毕。然后判断是否仍有未生成刚度矩阵的单元，如果还有，则进行下一个单元的循环计算，直至所有单元的刚度矩阵按照上述步骤生成完毕， 从而获得整个网格模型总体刚度矩阵的上三角矩阵及其刚度方程。最后，采用松弛预处理共额梯度法，求解整个网格模型的总体刚度方程，获得速度场的增量解，然后利用此增量解刷新网格模型在该时刻的速度场。网格模型中“广义相邻节点对”的数量由网格模型的拓扑结构决定，属于同一个单元的任意两个节点都可以组成一个“广义相邻节点对”，因此“广义相邻节点对”数量的确定，就是计算整个网格模型中，不重复的属于同一个单元的节点对数量，对于不同单元类型的网格模型，确定“广义相邻节点对”的原理相同，但具体方法有所不同。由于“广义相邻节点对”与网格模型的总体刚度矩阵上三角矩阵中非零子矩阵一一对应，可以由“广义相邻节点对”的数量确定总体刚度矩阵上三角矩阵中非零子矩阵的数量，进而确定非零元素的数量。然而由于网格模型的总体刚度矩阵采用上三角矩阵存储， 所以对角线上的每个非零子矩阵也只需存储上三角子矩阵。若设m为所求解问题的维数， 对角线非零子矩阵的非零元素为m(m+l)/2，而非对角线非零子矩阵数据为m2，由此可以确定网格模型总体刚度矩阵上三角矩阵的非零元素的数量。以由有限数量的八节点六面体单元构成的网格模型为例对上述计算过程加以说明。对于由八节点六面体单元组成的网格模型，其“广义相邻节点对”分为三个类型， 一是边相邻节点对；二是面相邻节点对；三是体相邻节点对。边相邻节点对由单元内相邻节点连线的数量决定，每条连线对应一个“广义相邻节点对”;面相邻节点对的数量由单元面的数量决定，每个面有两条对角线，每条面对角线对应一个“广义相邻节点对”；对于体相邻节点对，每个单元有四条体对角线，每条体对角线对应一个“广义相邻节点对”。设N1为节点数量(即为节点与自身生成的节点对数量)，N2为网格模型中相邻节点连线的数量(即为边相邻节点对的数量)，N3为网格模型中单元的面数(每个面有两个面相邻节点对，2队表示面相邻节点对的数量)，N4为网格模型中单元的数量(每个单元有两四个体相邻节点对， 4N4表示体相邻节点对的数量)，则网格模型“广义相邻节点对”的数量为
Ssm = N!+N2+2 X N3+4 X N4其中N1为节点自身构成的节点对的数量，该类型节点对对应的非零子矩阵行号与列号一致，位于对角线上，所以N1为对角线非零子矩阵的数量。对于三维问题(m = 3)，则每个对角线非零子矩阵需要存储的数据为m(m+l)/2 = 6，非对角线非零子矩阵需要存储的数据为m2 = 9。因此网格模型的总体刚度矩阵上三角矩阵非零元素的数量为Snz = 6 X N1+9 X (N2+2 X N3+4 X N4)下面以两个具体网格模型为例对上述计算过程作进一步说明。对于附图4所示的由一个六面体单元构成的网格模型,N1 = 8,N2 = 12,N3 = 6,N4 =1，采用上三角矩阵存储时，其非零子矩阵数量为Ssm = N1+ (N2+2 X N3+4 X N4) = 8+(12+2 X 6+1 X 4) =36图4所示的网格模型需要保存的总体刚度矩阵上三角矩阵非零元素数量为
Snz = 6 X N1+9 X (N2+2 X N3+4 XN4) = 6 X 8+9 X 28 = 300同理，对于附图5所示的由两个六面体单元构成的网格模型，N1 = 12、N2 = 20、N3 =IUN4 = 2。采用上三角矩阵存储时，其非零子矩阵数量为Ssm = N1+ (N2+2 X N3+4 X N4) = 12+ (20+2 X 11+2 X 4) = 62图5所示的网格模型需要保存的总体刚度矩阵上三角矩阵非零元素数量为Snz = 6 X N1+9 X (N2+2 X N3+4 XN4) = 6X 12+9X50 = 522确定了非零子矩阵和非零元素的数量，即可准确地分配网格模型的总体刚度矩阵上三角矩阵的存储空间，进而实现网格模型的总体刚度矩阵上三角矩阵的存储。下面以图5所示的由两个六面体单元构成的网格模型的总体刚度矩阵为例，说明本发明总体刚度矩阵的压缩存储优势。本发明所涉及的总体刚度矩阵存储方法中的J数组的数据长度为总体刚度矩阵非零子矩阵的数量，即Ssm = 62，而K数组的数据长度为单元数量12，则该网格模型的总体刚度矩阵存储数据量为62+12 = 74。如果采用目前常用的CSR(Compressed Sparse Row format)格式存储方法，则其J数组的数据长度为所有非零元素的数量，即Snz = 522，K数组的数据长度为总体刚度方程的维数，即为12X3 = 36，因此，总数据量为558。如果只比较两个辅助数组J和K，则本发明所涉及的总体刚度矩阵存储方法与CSR方法相比，其节约的存储空间为(558-74)/558 = 87%。如果浮点型数组A采用单精度浮点数保存，则本发明所涉及的总体刚度矩阵存储方法与CSR方法相比，总体刚度矩阵存储所节约的空间为 1-(74+552X 2)/(558+552X 2)=四％。可见，对于图5所示的例子，本发明节约辅助数组 J和K的存储空间达到87%，节约总体刚度矩阵存储空间四％，空间节约效果明显。下面通过一个具体的锻造实例，对本发明的有益效果做进一步说明。附图6a_图6d为一锻件在锻造过程中某一时刻采用不同的单元划分数量得到的有限元数值模拟结果，其中所有网格模型的单元全部为八节点六面体单元，图6a中的锻件网格模型含有1565个单元和2045个节点，图6b中的锻件网格模型含有3284个单元和4119 个节点，图6c中的锻件网格模型含有5339个单元和6557个节点，图6d中的锻件网格模型含有8530个单元和10217个节点。下面从内存占用和计算效率两个方面，将本发明的压缩存储迭代求解算法与半带存储直接求解的算法进行比较。其中所采用的计算机为CPU:Intel(R) Core (TM) 2, Duo, CPU E8400，3. OOGHz，内存3. 25，操作系统Windows XP。在内存占用方面，锻件网格模型的总体刚度方程的元素采用单精度浮点数保存， 即每个数据占用4个字节；辅助数组采用整形数据保存，即每个数据占用2个字节。在不同的节点和单元数量下，对采用不同方式存储的锻件网格模型总体刚度矩阵所占存储空间进行比较，结果如表1所示。表1在不同存储方法下锻件网格模型总体刚度矩阵存储空间的比较
权利要求
1.一种金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，其特征在于包括以下步骤(1. 1)首先采用现有CAD软件或者三坐标测量仪获得金属体积塑性成形问题的三维几何模型；采用现有网格生成方法生成成形工件的有限元网格模型；其次根据塑性成形工件的有限元网格模型拓扑结构，确定其网格节点的“广义相邻节点对”，建立“广义相邻节点对”与工件有限元网格模型总体刚度矩阵非零子矩阵之间的对应关系；(1.2)基于“广义相邻节点对”，储存非零子矩阵，开辟一个非零子矩阵列号存储数组J 即整型索引数组，记录与工件有限元网格模型中“广义相邻节点对”对应的总体刚度矩阵非零子矩阵的列号，开辟另一个整型索引数组K，记录工件有限元网格模型刚度矩阵上三角矩阵中每行首个非零子矩阵在数组J中的位置，确定工件有限元网格模型刚度矩阵中非零子矩阵的位置与数量，根据非零子矩阵与非零元素的对应关系，确定工件有限元网格模型刚度矩阵非零元素的位置和数量；(1. 3)计算工件有限元网格模型中每个单元刚度矩阵元素在工件有限元网格模型总体刚度矩阵中所在的非零子矩阵，利用索引数组J和K中的数据，确定单元刚度矩阵元素在网格模型总体刚度矩阵存储数组A中的位置；按照网格模型总体刚度矩阵的生成规则，将网格模型的每个单元刚度矩阵元素对号入座地组装到网格模型的总体刚度矩阵存储数组A 中；遍历网格模型中的所有单元后，最终生成金属体积成形问题网格模型的总体刚度矩阵存储数组；(1. 4)采用现有的对称正定系数矩阵松弛预处理共额梯度法，对网格模型的总体刚度方程即大型线性方程组进行求解，获得金属体积成形过程网格模型中所有节点的速度增量场的数值解。
2.根据权利要求1所述的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法， 其特征在于，所述步骤(1. 1)中，确定网格节点的”广义相邻节点对”的方法为(2. 1)在体积成形问题的有限元网格模型中，每个单元内的任意两个节点构成的节点对，即包括节点自身与自身构成的节点对，称为“广义相邻节点对”；将相邻节点有顺序的节点对称为“有序广义相邻节点对”；反之，称为“无序广义相邻节点对”，若网格模型中某个单元两个节点的全局节点编号为(I，J)，则构成两个“有序广义相邻节点对” ^7和云 ，若不考虑节点的顺序，则构成一个“无序广义相邻节点对”即认为^7和云 是同一个节点对，记为&^或I1 ；为了应用和表示的方便，可采取下标由小到大的形式来表示“无序广义相邻节点对”，即统一记为知(1彡J)；(2. 2)在网格模型中，每个“有序广义相邻节点对”对应网格模型的总体刚度矩阵中的一个与所分析问题维数相同的非零子矩阵，非零子矩阵的位置由“有序广义相邻节点对”的全局节点序号决定；(2. 3)在网格模型中，每个“无序广义相邻节点对”对应网格模型的总体刚度矩阵的上三角矩阵中的一个与所分析问题维数相同的非零子矩阵，非零子矩阵的位置由“无序广义相邻节点对”的全局节点序号决定；(2. 4)金属体积塑性成形过程有限元分析的总体刚度矩阵为对称矩阵，只存储其上的三角矩阵；网格模型中“无序广义相邻节点对”与网格模型的总体刚度矩阵的上三角矩阵中的非零子矩阵一一对应，并且节点对的全局节点编号与非零子矩阵的行号、列号对应，因此通过这种对应关系，由网格模型中所有的“无序广义相邻节点对”的数量和全局节点号确定网格模型刚度矩阵的上三角矩阵中非零子矩阵的数量和位置，即行号与列号，从而进一步确定组成非零子矩阵的非零元素的数量和位置。
3.根据权利要求1所述的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法， 其特征在于，所述步骤(1. 中，还包括以下步骤(3. 1)采用一个浮点型数组A按行顺序存储网格模型总体刚度矩阵上三角矩阵中的各非零元素；(3. 2)采用整形数组J记录与“无序广义相邻节点对”对应总体刚度矩阵上三角矩阵中的非零子矩的列号；(3. 3)采用整形数组K记录总体刚度矩阵上三角矩阵中每行第一个非零子矩阵在数组 J中的编号，并以此确定非零子矩阵在总体刚度矩阵上三角矩阵的行号，因为K[i]为总体刚度矩阵上三角矩阵中第i行第一个非零子矩阵在数组J中的编号，而K[i+1]为总体刚度矩阵上三角矩阵中第i+Ι行第一个非零子矩阵在数组J中的编号，由此在数组J中编号从 K[i]到K[i+1]-1对应的所有非零子矩都在第i行。
4.根据权利要求1所述的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法， 其特征在于，所分析问题的几何模型划分后的有限元网格模型为三维的六面体单元、四面体单元或其二者的混合单元模型；或者所分析问题的几何模型划分后的有限元网格模型为二维的四边形单元、三角形单元或其二者的混合单元模型。
5.根据权利要求1或4所述的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，其特征在于，所述步骤(1.3)中，对于网格模型中的一个单元，若该单元的节点数为η、 维数为m，记该单元η个节点的整体编号为I [i]，i = 1,2,…，n，单元刚度矩阵为U(mXn， mXn)，则该单元刚度矩阵的任意元素为Uijkl = u[(i-l) Xm+k] [(j-1) Xm+1]其中，i，j = l，2，…，η为该单元节点的局部编号，k，l = 1，2，…，m为物理量维数编号；确定Uijkl在网格模型总体刚度矩阵存储数组A中序号的过程如下(4. 1)确定Uijkl的整体节点编号I = I[i]，J = I[j]，其中1、为Uijkl所在的非零子矩阵Au中的在网格模型总体刚度矩阵上三角矩阵中的行号和列号；(4. 2)确定非零子矩阵A11之前非零子矩阵的数量=N1 = K[I]-1 ；(4. 3)确定非零子矩阵A11之前非零元素的数量NN1 = N1Xm2-(I-I) XmX (m_l)/2其中，m2为每个非零子矩阵非零元素的数量，m(m-l)/2为每行的对角非零子矩阵引起的非零元素的减少项；(4. 4)确定第I行非零子矩阵中第k行之前非零元素的数量NNik = (K[I+1]-K[I]) XmX (k_l)(4. 5)确定第I行中，第J列之前的非零子矩阵的数量；搜索数组J从标号K[I]到 Κ[Ι+1]中J[II] = J的数组标号II，第I行中，第J列之前的非零子矩阵数量为Nj = II-K[I]；(4. 6)确定Uukl在网格模型总体刚度矩阵中所在行中的编号 NNjl = SjXm-TT+1其中，TT为由于对角非零子矩阵引起的非零元素的减少项，TT = k(k-l)/2 ； (4. 7)确定Uukl在网格模型总体刚度矩阵储存数组A中的编号 NNljkl = NNj+NNj.+NNj!(4. 8)根据网格模型总体刚度矩阵的生成规则和Uijkl在网格模型总体刚度矩阵储存数组A中的位置编号NNukl，将Uijkl组装或累加到网格模型的总体刚度矩阵储存数组A中；(4.9)重复步骤(4. 1) 8)，遍历网格模型中的所有单元刚度矩阵中的每一个元素，即可获得网格模型的总体刚度矩阵一维存储数组A。
6.根据权利要求1所述的金属体积塑性成形中有限元分析刚度矩阵存储与生成方法， 其特征在于分析网格模型所采用的数值分析方法最终归结为稀疏、对称总体刚度矩阵的存储、生成及其刚度方程即线性方程组求解，所分析的网格模型的总体刚度方程求解参量为全量或者增量。
全文摘要
本发明涉及一种金属体积塑性成形中有限元分析刚度矩阵存储与生成方法，可显著节约刚度矩阵存储空间和提高刚度方程求解效率，它包括以下步骤(1.1)根据体积塑性成形工件有限元网格模型，确定网格节点的“广义相邻节点对”及其与网格模型总体刚度矩阵中非零子矩阵的对应关系；(1.2)基于“广义相邻节点对”，确定网格模型总体刚度矩阵中非零子矩阵的位置与数量，存储非零子矩阵；(1.3)确定网格模型每个单元刚度矩阵元素在总体刚度矩阵中所在的非零子矩阵及其在总体刚度矩阵存储数组中的位置，将其组装到网格模型的总体刚度矩阵存储数组中；(1.4)采用对称正定系数矩阵松弛预处理共额梯度法，对网格模型的总体刚度方程进行求解。
文档编号G06F17/50GK102184298SQ20111013007
公开日2011年9月14日 申请日期2011年5月18日 优先权日2011年5月18日
发明者王忠雷, 赵国群, 马新武 申请人:山东大学