4.3  解题方略

4.3.1  结构有限元分析所用单元的类型选择

为了计算结果与实验结果相对照以验证计算模型的正确与否，所选择的单元其计算应力应能以总体坐标的形式来表示。叠片在轴向变形时，其位移量通常在1---2mm（角向位移时，许多节点的情况亦如此）这大大超过叠片的厚度（0.2—0.5mm )，因而其有限元应力分析应为几何非线性分析。基于叠片的几何特性，我们选用三角形板单元来模拟它的大部分区域，而在螺栓孔附近、叠片被垫片由预紧力压住，形成接触面的区域我们采用三维体元来模拟，这是因为在这种情况下，亦即存在着三维接触元的情况下，APOLANS程序指定必须使用三维体元，虽然在实际的几何尺寸中，单一三维体单元上厚度与长、宽尺寸差别较大，但在实际计算中，程序能够接受这种尺寸比例的体元，因此我们可以认为这种形式是适宜的。当然，最重要的还是体元和板元的交界部分如何对接。

如上图网格状阴影部分为接触区，即用体元来表示的区域，叠片的其余部分均用板元来表示，对于上图所指外圈处，亦即三维体元与板元的交界处，我们采用了如下图所示对接触方法。

如图所示，我们以12节点体元来模拟螺孔中心带磨擦部分的叠片实体，右边为三角形板元，其中第9、第10两个节点为三维体元和三角形板元共享的节点，这样就可以将两个部分结合在一起了。

当然，在三维体元和三解形板元的边界对接上是有很多种方法的，为了确保我们这种对接方式是最为有效、可靠的，我们做了大量的计算工作，尝试了诸多对接计算模型的验算，最后得出结论，只有上图所表示的接法才能够保证在载荷相同的情况下（即原载荷加在板元上，现在加在体元上），板元内部各点上所受力、力矩的方向和数量级才能保持一致，其余的各种方法，要么是力的数量级别产生较大差异，要么是力、力矩的方向发生变化，由此使我们最后认定以图4.3.2所示的对接方法是可反映实际情况的，后面将计算结果与实验数据相印证也证明了这种方法的可行性。

4.3.2  静力平衡方程的解法

对于静力问题最后建立的控制方程为代数方程组，对线性问题为线性代数方程组，对非线性问题为非线性代数方程组。

4.3.2.1  线性代数方程组的解法

一个线性方程组可写成

[K][X]=[R]                (4.3.1)

其中[K]为刚度阵，[X]为位移向量，[R]为载荷向量。

大多数情况下[K]阵是对称、正定的，所以有三角化分解：

[K]=[T]^T[T]                (4.3.2)

[T]为上三角阵，再进一步分解可得

[K]=[T]^T[D][T]                (4.3.3)

此时[T]为单位对角元的上三角阵，[D]为对角矩阵，由此得；

[T]^T[D][T][X]=[R]                (4.3.4)

即[T]^T[D][Z]=[R]                (4.3.5)

[Z]=[T][X]                (4.3.6)

式（4.3.5）为前代，(4.3.6）为回代，由此而计算出位移向量X ，有限元法形成的K矩阵是对称、正定、稀疏的，非零元主要分布在矩阵对角线两侧，一维变带宽存贮方法利用这一点，各列顺序存贮，但仅存贮第一个非零元以下的元素，由于对称仅存上三角矩阵，当实际工程模型较大时，内存放不下总刚度阵[K] , 则矩阵将被分块存贮于磁盘上，即程序规定给用户一部分内存空间A 作为工作数组用，不同计算阶段分配不同存贮内容，在组集总刚度阵时，存各元素信息及分块的总刚度阵，内、外存交互求解时放两个总刚度阵的块。总刚度阵即按A/2，数组尺寸自动分成若干块，存入磁盘文件。分解时按序调各块进入内存。

4.3.2.2  非线性代数方程组的解法

在叠片的带摩擦元的强度计算中，我们采用了带线性搜索的完全牛顿迭代方法。

一、牛顿拉斐逊法

结构的总位能为：

Π= U(q)-{q}^T{F}          (4.3.7)

其中，U(q)为结构的应变能，{q}^T{F}为外载所做的功。

由虚功原理使上式总位能变分为零可得有限元平衡方程：

Π=[W（q）][P(q)][F]=[K(q)][q]-[F]          （4.3.8）

其中，[P(q)]是与应力等效的节点力，[F]为外载。

我们有第i 次迭代得到 的［K_i（q_i）］可求到第i+1 次的未知变量｛q_i+1｝ 。

｛q_i+1｝={q_i′}-ω_i[K_i]^-1([K_i]{q_i}-{F})

为超松驰因子，用以加速收敛，即当我们有一个近似解{q_i}时，可用截断的阶泰勒级数得到一个改进的解：

上式左边称为雅可比矩阵，右边第一项为载荷位移空间的切线刚度，它是对称的，所以有：

[Δq_i]=-[K_i]^-1_T（[P（q_i）]-[F]）          ( 4.3.12)

其中（）内为第i 次迭代不平衡力。

二、线性搜索加速

在求解方程（4.3.8)时，，寻求一个n 阶矢量{q } ，实际上是在n 维空间找到一个方向和步长来满足该方程。对于非线性方程组，要经过多次迭代，每次迭代均要找到一个方向和步长，逐步达到平衡解，为了加速迭代的收敛，在由式（4.3.12）决定的力向上寻求一个优化的步长，即将｛△q｝乘以某一标量η而使迭代加快。决定η 的原则是使能量增量φ最小化，即：

实际上直接按上式寻找η很费事，每寻找一个｛△q｝，必须解一次方程（4.3.10）所以又有近似方法更简便，当不等式

|△{△q_i+1}^T{△F(q_i+△q_i+1)}|≤α|{△q_i+1}^T{△F(q_i)}|          (4.3.15)

成立时，即可得到满意的收敛率，其中α是收敛容差，一般在0.5-0.9之间，如上式不成立，则使用线性搜索以改善收敛速度。

|△{△q_i+1}^T{△F(q_i+η△q_i+1)}|≤α|{△q_i+1}^T{△F(q_i)}|          (4.3.16)

可用线性插值寻找η。

三、收敛准则

运用迭代法存在两次迭代间满足什么样的准则可以终止本次迭代问题，目前有三种收敛准则。

（一）位移准则

||△q_i+1||＜=α_D||q_i+q_i+1||          (4.3.17)

其中α_D为位移收敛容差，一般取0.1%--0.5%

取得是欧几里德范数

（二）不平衡准则

||△F_i||＜=α_D||F||

（三）能量准则

{△q_i}^T({F}{P_i})＜=α_E(U_i+T_i)

即把每次迭代时内能的增量（也就是不平衡力在位移增量上所做的功）同当前的总能量相比较。

其中，U_i为势能，T_i为动能

α_E是预定的能量容差。

单纯使用位移准则或不平衡准则不很理想，有些情况不会出现较大偏差，能量准则同时控制位移和力，由此达到平衡，是本文所采用的方法。

上一页

下一页