CN106126823B - 一种基于提高迭代法稳定性和收敛性的位移求解方法 - Google Patents

Abstract

本发明公开了一种基于提高迭代法稳定性和收敛性的位移求解方法，该方法基于有限元网格划分中包含奇异单元模型，从而单元刚度矩阵中会出现非对角占优的刚度矩阵为基础，对于非奇异单元及奇异单元中的非奇异节点采用传统的迭代技术进行处理，对于奇异单元中的畸形节点，可以在单元上进行改进，比如以预条件共轭梯度法为基础，采用体积作为惩罚函数进行预处理；本发明对非对角占优的有限元方程组进行修正，能够达到减少迭代次数，收敛速度更快，整个方法更加稳定的效果，其次采用先进行单刚和位移计算相乘后再叠加的方法，从而避免了组合总刚的计算，该计算量也是很庞大的；采用单刚的数据结构，方便数据的更新，并有利于并行计算。

Description

一种基于提高迭代法稳定性和收敛性的位移求解方法

技术领域

本发明涉及结构力学技术领域，具体涉及一种基于提高迭代法稳定性和收敛性的位移求解方法。

背景技术

通常在采用有限单元法对碰撞冲击问题进行数值分析时，特别是模拟高速冲击载荷下的结构响应时，例如弹打靶、爆炸过程中，一个实体变形前和拉伸变形后的情形如图6所示，节点A在进行弹打靶或爆炸前的位置如图1，一旦进行弹打靶或爆炸后，这个实体就发生了变形，继而点A的位置和方向就发生了变化，变成节点A’而线段

就是节点A在变形过程中产生的位移，那么该节点的位移怎么求？

如果要想描述整个物体位移，就要用有限单元的方法来描述，也即是用多个节点，及连接多个节点的网格(即是有限元网格)来代替原来的物体，那么求多个节点的位移就近似成求该物体的位移，也即是有限元的位移，那么该位移如何求？可以通过求解下面的有限元方程组来进行求解

Ku＝F (1)

其中K是单元刚度矩阵，F是物体在发生变形时所受到的力，u就是物体变形后的位置，即是位移。

有限元网格往往会随着材料的剧烈变形而发生畸变，导致计算中出现非常畸形的网格如图3所示。或者对于比较复杂的结构在划分网格的过程中，就会形成非常畸形的网格，那么这个非常畸形的网格就是奇异的单元，从而方程组(1)的刚度矩阵就会出现非对角占优，甚至接近奇异的现象，从而方程组(1) 就会出现条件数较差，求解困难的问题。那么对于包含这种非常畸形的网格的位移如何求解？现有的求解器方法如下：

1)采用总体刚度矩阵的奇异值分解方法求出矩阵的奇异值，再进行求矩阵的广义逆，从而求出位移，当矩阵很大时，这种方法计算量非常庞大，不适合在工程计算中应用。

2)采用总体刚度矩阵的预条件共轭梯度法，在采用预条件共轭梯度法求解位移时，针对系数矩阵选取不同的预优矩阵M，使矩阵的特征值分布更为集中，降低矩阵条件数，改善矩阵病态特性；以达到提高收敛速度的目的；

3)采用单元刚度矩阵的预条件共轭梯度法，对于有限元方程组的求解，通常在不形成总体刚度矩阵情况下，可以将大规模总体结构的计算分解成相对独立的中小规模的结构计算，且使这些计算可由各台相对独立的处理机并行完成，其优点是有效提高计算时间，该方法在单元刚度矩阵发生变化时可以直接更新数据，无需重新组合成总刚。

随着有限元技术的发展，越来越多的有限元问题涌现出来，再用以上三种方法求解有限元问题时，对于存在非常畸形的网格的求解位移时，就会出现迭代次数偏大，收敛速度慢，不收敛，有的甚至得到的解与真解误差较大的问题，那么原来方法的求解器就会失效。

发明内容

针对现有技术的不足，本发明的目的是提供一种基于提高迭代法稳定性和收敛性的位移求解方法，该方法基于有限元网格划分中包含奇异单元模型，从而单元刚度矩阵中会出现非对角占优的刚度矩阵为基础，对于非奇异单元及奇异单元中的非奇异节点采用传统的迭代技术进行处理，对于奇异单元中 的畸形节点，可以在单元上进行改进，比如以预条件共轭梯度法为基础，采用体积作为惩罚函数进行预处理。

为实现以上目的，本发明采取了以下的技术方案：

一种基于提高迭代法稳定性和收敛性的位移求解方法，以预条件共轭梯度法为基础，采用体积作为惩罚函数进行预处理其包括：

步骤1、取初始值u₀,矩阵乘向量部分采用单元刚度矩阵乘以单元向量的做法，这样避免了合成总体刚度矩阵，则初始残余力向量

r₀＝F-∑_eK_eu_0e

这里F是载荷，K_e是单元刚度矩阵，u_0e初始单元位移，r₀是残余力向量；

步骤2、解以下方程组

Mz₀＝r₀

求出z₀

及后面迭代步中，同样解以下方程组

Mz_k＝r_k

求出z_k。

这里M是总体刚度矩阵(简称总刚)的预处理矩阵，r_k是迭代步中的残余力向量。

所述步骤2包括以下步骤：

步骤21、对于非奇异单元及奇异单元中的非奇异节点，选取预条件M，这里的M可以是Jacobi预条件，也可以是松弛迭代法中的分裂因子的预优矩阵等等，也即是这里的M可以是任意的预优矩阵，下面以Jacobi预条件为例：

M＝diag(a₁₁,a₂₂,…a_nn)

并记

d＝(a₁₁,a₂₂,…,a_nn)

其中，a_ii是总体刚度矩阵的对角元素，则，每迭代步中的残余力向量的求解如下：

z_k＝r_k/d_k

这里k是第k次迭代；

步骤22、对于单元中的奇异节点，采用对单元中的奇异节点进行单独处理的方法，以改善刚度阵中非对角占优的情况，下面具体以采用体积惩罚函数方法进行改善的方法为例，体积惩罚函数求法如下：

首先找到奇异节点，及其被哪些单元包含，求出该单元的质心，并找到奇异夹角所在的面，然后求出各个面的质心。

从而求出这些面与对应单元的质心组成的立方体的体积，这就得到单元奇异节点的体积惩罚函数，我们记为v_1e,v_2e,…,v_ke，从而

z_ki＝∑_jv_kjez_kje

这里

z_kje＝r_kje/d_je

i是奇异节点号，j是奇异节点所在单元的序列号；k是迭代次数，r_kje是单元残向量，d_je是单元刚度矩阵对角元素，e是单元号；

步骤3下面求出位移向量的共轭梯度方向：p_k,k＝1,2,…n是迭代步；

if k＝1

p₀＝z₀

else

p_k＝z_k-1+β_kp_k-1

endif

步骤4、对迭代步k＝1,2,…计算

y＝∑_eK_ep_ke

u_k＝u_k-1+α_kp_k

r_k＝r_k-1-α_ky

再求出z_k；

步骤5、直到||r_k||<ε，也即是该方法收敛，则输出最终位移u_k，结束。

本发明与现有技术相比，具有如下优点：本发明首先对迭代方法中的采用单元刚度矩阵的形式进行计算；其次对于单元中奇异节点进行特殊处理，这样使对非对角占优的矩阵进行修正，减少迭代次数，收敛速度更快，整个算法更加稳定；计算实施步骤简洁，不需要修改迭代法中其它计算模块。同时无需组合总体刚度矩阵，这样避免了组合总刚后再乘以位移的计算，而是采用先进行单刚和位移计算相乘后再叠加的方法，从而避免了组合总刚的计算，该计算量也是很庞大的；再次，采用单刚的数据结构，方便数据的更新，而采用总刚的数据结构，每次都需先计算出总刚才能更新.采用这种方法，可以有效解决在有限元网格划分中存在奇异单元的位移求解。

附图说明

图1为本发明一种基于提高迭代法稳定性和收敛性的位移求解方法的流程示意图；

图2为奇异节点与所在单元的示意图；

图3为三维悬臂梁模型的示意图；

图4为迭代次数对比曲线；

图5为不同方法残向量r的对比图；

图6为实体变形前和拉伸变形后的示意图。

具体实施方式

下面结合具体实施方式对本发明作进一步的说明。

本发明可以提高迭代法稳定性和收敛性，该方法基于有限元网格划分中包含奇异单元模型，从而单元刚度矩阵中会出现非对角占优的刚度矩阵为基础，对于非奇异单元及奇异单元中的非奇异节点采用传统的迭代技术进行处理，对于奇异单元中的畸形节点，可以在单元上进行改进，比如以预条件共轭梯度法为基础，采用体积作为惩罚函数进行预处理，具体操作如下：

步骤1、取初始值u₀,矩阵乘向量部分采用单元刚度矩阵乘以单元向量的做法，这样避免了合成总体刚度矩阵，则初始残余力向量

r₀＝F-∑_eK_eu_0e (1)

这里F是载荷，K_e是单元刚度矩阵，u_0e初始单元位移，r₀是残余力向量；

步骤2、解以下方程组

Mz₀＝r₀ (2)

求出z₀

及后面迭代步中，同样解以下方程组

Mz_k＝r_k (3)

求出z_k。

这里M是总体刚度矩阵(简称总刚)的预处理矩阵，r_k是迭代步中的残余力向量。

步骤2采用以下方法：

对于非奇异单元及奇异单元中的非奇异节点，采用传统的预条件；对于奇异节点，先求出奇异节点所在单元的体积惩罚函数，再根据奇异节点所在的单元进行相乘，奇异节点所在单元的体积惩罚函数值越大，那么所占的比例也越 大，那么该单元对这个奇异节点的贡献也就越大，反而，奇异节点所在单元的体积惩罚函数值越小，那么所占的比例也越小，那么该单元对这个奇异节点的贡献也就越小；

所述步骤2包括以下步骤：

步骤21、对于非奇异单元及奇异单元中的非奇异节点，选取Jacobi预条件

M＝diag(a₁₁,a₂₂,…a_nn) (4)

并记

d＝(a₁₁,a₂₂,…,a_nn) (5)

其中，a_ii是总体刚度矩阵的对角元素，则，每迭代步中的残余力向量的求解如下：

z_k＝r_k/d_k (6)

这里k是第k次迭代；

另一方面此处的预条件M可以选取SSOR-CG方法，及其它预条件做预处理。

SSOR-CG方法是指对A＝D-L-U的自然分裂，取M＝(D-ωL)D^-1(D-ωL)^T.它对应于超松弛迭代法通用格式中的分裂阵因子M，ω是松弛因子(0<ω<2)。

步骤22、对于奇异节点，采用体积惩罚函数方法，体积惩罚函数求法如下：

首先找到奇异节点，及其被哪些单元包含，求出该单元的质心，并找到奇异夹角所在的面，然后求出各个面的质心。

如图2：首先找到奇异节点B，及其被哪些单元包含，求出该单元的质心，并找到奇异夹角所在的面:四边形ABHD和三角形BHC，

然后求出各个面的质心，质心坐标求法运用以下公式：

这里x_i、y_i、z_i,i＝1,2,…n是四边形ABHD和三角形BHC各顶点的坐标

从而求出这些面与对应单元的质心组成的立方体的体积，这就得到单元奇异节点的体积惩罚函数，我们记为v_1e,v_2e,…,v_ke，从而

z_ki＝∑_jv_kjez_kje (8)

这里

z_kje＝r_kje/d_je (9)

i是奇异节点号，j是奇异节点所在单元的序列号；k是迭代次数，r_kje是单元残向量，d_je是单元刚度矩阵对角元素，e是单元号，

步骤3下面求出位移向量的共轭梯度方向：p_k,k＝1,2,…n是迭代步；

if k＝1

p₀＝z₀

p_k＝z_k-1+β_kp_k-1

endif

步骤4、对迭代步k＝1,2,…计算

y＝∑_eK_ep_ke

u_k＝u_k-1+α_kp_k

r_k＝r_k-1-α_ky (11)

再根据公式(6)和(9)求出z_k。

步骤5、直到||r_k||<ε，也即是该方法收敛，则输出最终位移u_k，结束。ε是收敛准则，一般取一个接近于零的数。

以下实施例是一个三维悬臂梁问题，结构的整体几何尺寸为1*1*1m，如图3所示，将结构划分为两个实体单元，其中一个为六面体实体单元，一个为五面体实体单元，单元尺寸如图，对应ANSYS单元类型为SOLID65。模型采用一种材料，材料参数为：弹性模量E＝34.5GPa，泊松比μ＝0.2。模型如图2所示。

以l＝0.025*a为例，将上述模型运用一种提高迭代法稳定性和收敛性的方法与总体刚度矩阵的预条件共轭梯度法(这里预条件仍旧选取Jacobi预条件)求解位移结果对比见表1所示。

表1位移求解结果与总体刚度矩阵法求解的结果对比

对于上述模型更改l占a的比例，即是改变奇异单元的畸形程度，也即是比例越小越畸形，从而形成的刚度矩阵就是非对角占优矩阵，以l₁＝1/10*a，l₂＝1/20*a，l₃＝1/30*a，l₄＝1/40*a，l₅＝1/60*a，l₆＝1/80*a，l₇＝1/100*a为例，查看运用将上述模型运用一种基于提高迭代法稳定性和收敛性的位移求解方法，与运用原来的总体刚度矩阵的共轭梯度算法和原来的单元刚度矩阵预条件共轭梯度法(都采用Jacobi预条件)，求解位移结果的迭代次数对比见图4所示。

从图4可以看出，在对于包含奇异楔形单元的单元刚度矩阵的有限元方程组求解，当楔形角度越来越尖锐时，也就是楔形体的畸形度越来越大时，用一种基于提高迭代法稳定性和收敛性的位移求解方法的迭代次数远远小于运用原来的总体刚度矩阵的共轭梯度算法和原来的单元刚度矩阵预条件共轭梯度法，从而对于包含奇异楔形单元的单元刚度矩阵的有限元方程组求解，当畸形度比较大的时候，运用一种基于提高迭代法稳定性和收敛性的位移求解方法具有较好的稳定性，及收敛速度比较快的的优点。

图5是不同方法残余力向量的比较，这里以l₄＝1/40*a为例，从图中可以看出，运用一种基于提高迭代法稳定性和收敛性的位移求解方法，在最后几步，残余力向量下降的非常快，从而表明该算法收敛的也比较快。

上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。

Claims (2)

1.一种基于提高迭代法稳定性和收敛性的位移求解方法，其特征在于，用于采用有限元方法模拟高速冲击载荷下的结构时，对于奇异单元中的畸形节点，可以在单元上进行改进，采用体积作为惩罚函数进行预处理，包括步骤：

步骤1、取初始值u₀，矩阵乘向量部分采用单元刚度矩阵乘以单元向量的做法，这样避免了合成总体刚度矩阵，则初始残余力向量

r₀＝F-∑_eK_eu_0e

上式中，F是载荷，K_e是单元刚度矩阵，u_0e初始单元位移，r₀是残余力向量；

步骤2、解以下方程组

Mz₀＝r₀

求出z₀，及后面迭代步中，同样解以下方程组求出z_k：

Mz_k＝r_k

所述步骤2包括以下步骤：

步骤21、对于非奇异单元及奇异单元中的非奇异节点，选取Jacobi预条件

M＝diag(a₁₁，a₂₂，…a_nn)

并记

d＝M＝(d₁，d₂，…，d_n)

则，残差的求解如下：

z_k＝r_k/d_k

这里k是第k次迭代，此处的预条件选取Jacobi-CG方法；

步骤22、对于奇异节点，采用体积惩罚函数方法，体积惩罚函数求法如下：

首先找到奇异节点，及其被哪些单元包含，求出该单元的质心，并找到奇异夹角所在的面，然后求出各个面的质心，质心坐标求法运用以下公式：

从而求出这些面与对应单元质心组成该单元的体积，这就得到体积惩罚函数，记为v_1e，v_2e，…，v_ke，从而

z_ki＝∑_jv_kjez_kje

这里

z_kje＝r_kje/d_je

i是奇异节点号，j是奇异节点所在单元的序列号；r_kje是单元残差，d_je是单元刚度矩阵对角元素。

2.根据权利要求1所述的基于提高迭代法稳定性和收敛性的位移求解方法，还包括步骤：

步骤3、对k＝1，2，…计算

y＝∑_eK_ep_ke

u_k＝u_k-1+α_kp_k

r_k＝r_k-1-α_ky

z_ki＝∑_jv_kjez_kje

步骤4、直到||r_k||＜ε，输出u_k，结束。

Images