基于Householder变换的无约束结构静力分析方法

摘要

本发明公开了基于Householder变换的无约束结构静力分析方法，其包括：步骤1、构建结构的刚体位移模式矩阵X；步骤2、根据刚体位移模式矩阵X构造6个相应的n×n阶Householder矩阵P

说明书

技术领域

本发明涉及有限元的结构分析，具体涉及一种无约束状态下的结构力学问 题的有限元分析方法。

背景技术

通常做线性静力分析需要保证结构没有刚体位移，否则求解器没有办法 计算。但是在空间中很多问题没有足够的约束，如飞机在飞行时，轮船在航 行时，或者货品在起重机上吊装时，要想计算结构上的应力分布，需要采用 特殊的无约束结构静力分析方法。现有的分析方法有以下几类。

1）惯性释放法，其基本思路是用结构的惯性(质量)力来平衡外力。尽 管结构没有约束，分析时仍假设其处于一种“静态”的平衡状态。采用惯性释 放功能进行静力分析时，需要对一个节点进行6个自由度的约束(虚支座)。 针对该支座，程序首先计算在外力作用下每个节点在每个方向上的加速度， 然后将加速度转化为惯性力反向施加到每个节点上，由此构造一个平衡的力 系(支座反力等于零)。求解得到的位移描述所有节点相对于该支座的相对运 动。

2）动力松弛法，动力松弛法是用显式动力学原理将结构的静力学问题 转化为动力学问题近似求解线性或非线性系统平衡状态的数值方法，在计算 时添加虚拟的质量和阻尼，构成质量-刚度-阻尼的平衡体系，添加的阻尼应 足够大，以使得结构的动能在动力松弛分析完成后足够小。

随着有限元技术的推广，各种各样的无约束问题涌现出来，在用传统有 限元法对其进行求解时，会遇到刚度阵奇异的问题。而现有的无约束问题求 解法如惯性释放法需要引入虚支座的概念，由于虚支座不是真实的位移约 束，与真实情况相比会产生较大的误差，且需要求解方程，效率有时偏低； 动力松弛法由于采用显式分析方法，因此效率高但精度低。

对于有限元方程组的求解，通常结构有限元分析中形成的总体刚度阵都 具有稀疏质，对于这类刚度阵对应的线性方程组非常适合采用共轭梯度法进 行求解，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解 大型非线性最优化最有效的优化方法之一。在各种优化方法中，共轭梯度法 是非常重要的一种。其优点是所需存储量小，具有有限步收敛性，稳定性高， 而且不需要任何外来参数。

发明内容

针对现有技术存在的问题，本发明的首要目的是提供一种基于 Householder变换的无约束结构静力分析方法，该方法以空间无约束模型为 基础，根据Householder变换理论，用结构的6个的刚体位移模式构造6个 Householder矩阵，用此Householder矩阵去除结构的刚体模态，并采用修正后 的共轭梯度法求解去刚体模态后的刚度方程。

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

基于Householder变换的无约束结构静力分析方法，其包括：

步骤1、构建结构的刚体位移模式矩阵H：

H=[h₁ h₂ h₃ h₄ h₅ h₆]    （1）

其中，h₁～h₃为结构的平动刚体位移模式，h₄～h₆为结构的转动刚体位移模式， h₁～h₆均为n维列向量，n为结构整体自由度；

步骤2、根据h₁～h₆构造6个相应的n×n阶Householder矩阵P_i，i=1，2……6； 设结构离散化后计算得到的未指定位移边界条件的有限元整体刚度阵为K，则 此刚度阵K具有6重零特征值，且具有6个对应于零特征值的正交特征向量， 这6个特征向量张成的空间为刚度阵K的零空间，这个零空间与6个结构刚体 位移模式张成的空间一致，即满足

K*h_i=0    （2）

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

步骤21、首先将平动刚体位移模式h₁单位化，且与单位向量e₁（e₁为第一 个元素为1，其它元素都为0的单位向量）一起用来构造Householder矩阵P₁， 使得：

P₁*h₁=e₁    （3）

所述Householder矩阵P₁为正交对称矩阵；

使用Householder矩阵P₁对刚度阵K进行相似变换，将刚度阵K的第一行和 第一列均化为0，即：

K₁=P₁*K*P₁    （4）

其中，K₁是刚度阵K的第一行和第一列均为0的相似刚度矩阵；

步骤22、构造Householder矩阵P₂，使用Householder矩阵P₂对相似刚度矩 阵K₁进行相似变换，可使相似刚度矩阵K₁的前两行和前两列均为0，即：

K₂=P₂*K₁*P₂    （5）

其中，K₂是相似刚度矩阵K₁和刚度阵K的前两行和前两列均为0的相似刚度矩 阵；

构造Householder矩阵P₂的方法是：

根据公式（2）以及正交对称矩阵特性可知：

K*h_i=0    （6）

P₁*P₁=I_n    （7）

其中，I_n为n×n阶单位矩阵；

由公式（6）、（7）可得：

K₁*P₁*h₂=0    （8）

令：

h₂′=P₁*h₂    （9）

h₂′为新构造的平动刚体位移模式，由于相似刚度矩阵K₁的第一列元素均为 零，且根据特征向量两两相互正交的特点，将h₂′的第一个元素赋值为零，即：

h₂′(1)=0    （10）

将新构造的平动刚体位移模式h₂′单位化与单位向量e₂一起构造 Householder矩阵P₂：

P₂*h₂′=e₂    （11）

步骤23、根据步骤22中构造Householder矩阵P₂的方法逐个对平动刚体位 移模式h₃以及转动刚体位移模式h₄～h₆进行操作，构造新的平动刚体位移模式 h₃′以及新的转动刚体位移模式h₄′～h₆′，进一步构造Householder矩阵P₃～P₆；

步骤3、利用步骤2中构造的Householder矩阵P_i对结构原始刚度矩阵做正 交相似变换，得到去除刚体模态后的结构的刚度矩阵K_p：

K_p=P₆P₅P₄P₃P₂P₁KP₁P₂P₃P₄P₅P₆    （12）

简化得：

K_p=P^TKP    （13）

其中，P=P₁P₂P₃P₄P₅P₆；

步骤4、根据通过有限元离散并用常规方法得到结构静力问题分析的线性方 程：

KU=F    （14）

其中K为刚度阵，U为待求解的位移向量，F为外载右端项；

通过Householder矩阵的正交相似变换去除结构总体刚体位移模态后上述方 程（14）为：

K_PU_P=F_P    （15）

其中，K_P=P^TKP，U_P=P^TU，F_P=P^TF，P=P₁P₂P₃P₄P₅P₆。

步骤1包括：假设将结构整体沿x方向平移一个单位位移，则在物体内不 会产生任何应变，而此时结构的刚体位移模式为只在第6j+1个自由度上有位移 值1其它都为0，具体可写为：

h₁=[1 0 0 0 0 0 1 0 0 0 0 0 … 1 0 0 0 0 0]^T    （16）

同理，将结构整体分别沿y方向和z方向平移一个单位位移，构造平动刚 体位移模式h₂和h₃。

假设将结构上的节点m绕z轴旋转一个微小角度θ，此时m点x方向的位 移为：

$\left(\begin{array}{c}u={x_m}^{\prime }-x_m\\ =-\delta \sin \alpha \\ =-r\sin \mathrm{\alpha \theta }\\ =-y_m\theta \end{array}\right)---\left(17\right)$

同理，y方向的位移为：

v=x_mθ    （18）

z方向的位移为零：

w=0    （19）

此时m点的位移模式可写为：

h₆=[-y_m θ x_mθ 0 0 0 1θ]^T    （20）

在上式中约去θ，并将其扩展到整个结构，即结构整体绕z轴旋转一个微小 角度θ，而在物体内不产生任何应变，此时结构的刚体位移模式可写为：

h₆=[-y₁ x₁ 0 0 0 1 -y₂ x₂ 0 0 0 1 … -y_n/₆ x_n/₆ 0 0 0 1]^T    （21）

其中，x₁～x_n/6是节点1到节点n/6的初始x轴坐标，y₁～y_n/6是节点1到节点n/6 的初始y轴坐标，1≤m≤(n/6)；

同理，将结构整体分别绕x轴和y轴旋转一个微小角度θ，构造转动刚体位 移模式h₄和h₅。

Householder变换理论可以变换表述为：对于任意给定的两个n阶单位向量 a和b，定义向量v满足：

c=a+b    （22）

v=c/|c|    （23）

其中|c|为向量c的长度。

Householder矩阵可以构造为

P_i=2v_iv_i^T-I_n    （24）

其中v_i为归一化后的n维Householder向量。矩阵P为对称正交矩阵，满足 如下三个条件：

P_i=P_i^T    （25）

P_i=P_i^-1    （26）

P_i*a=b    （27）

在本方案中，Householder矩阵P_i均为对称正交矩阵，而由这些Householder 矩阵P_i组成的Householder矩阵P则只满足正交性，不满足对称性。

步骤5、采用修正共轭梯度法求解线性方程组（15）：

去除刚体模态后的结构的刚度矩阵K_p中前六行的元素均为零，因此相应的 要把变换后的外载右端项的前六个元素设置为零，即：

F_P(1)=F_P(2)=...=F_P(6)=0    （28）

为避免直接计算去除刚体模态后的结构的刚度矩阵K_p，而造成占用大量空 间，并避开共轭梯度法求解的过程中计算矩阵向量乘法的时候需要花费大量的 计算时间，对常规的共轭梯度法进行相应的修正，以便充分发挥共轭梯度法的 优点。

对于占用空间大的问题，可以仅需存储具有稀疏特性的刚度阵K和 Householder变换矩阵P_i以代替去除刚体模态后的结构的刚度矩阵K_p，由 Householder变换的性质可知

P_i=2v_iv_i^T-I_n    （29）

更进一步地也不需要存储6个Householder变换矩阵，仅需要存储Householder 变换矩阵P_i对应的向量vi，避免存储满阵K_p占用大量空间（其中，刚度阵K具 有稀疏特性，其为稀疏矩阵，而K_p为满阵）。这样的做法可以减少大量不必要 的计算机存储空间。

另一方面，在共轭梯度法中计算矩阵向量乘法，可以采用分步乘积的方式， 具体形式如下：

K_p*p_k=P^T*K*P*p_k    （30）

其中，p_k为共轭梯度法中的修正方向。从上式可以看出，上式的矩阵向量乘法 可以分解为三步进行：

s_k=P*p_k    （31）

s_k=K*s_k    （32）

s_k=K_p*p_k=P^T*s_k    （33）

从上面的公式（32）可以看到，做矩阵向量乘法运算的时候，总体刚度矩阵 K的稀疏特性可以大大提高计算效率。

更进一步，计算P^T*s_k时也可以均分解成6次乘法进行，每一次均对向量 前乘一个Householder变换矩阵，由Householder变换矩阵的特性可知，对于任 意Householder矩阵P和任一向量g，其乘积可以表示为

P*g=(2vv^T-I_n)g=2v(v^T_g)-g    （34）

由上式可以看出Householder矩阵P与向量做乘积运算时，实际上可以用一 次点积、一次标量向量乘法和一次向量加减法运算代替，这无疑会大大减少计 算量，特别是对于大规模、节点数庞大的结构分析问题，上述分析方法可以大 大提高计算效率。

本发明与现有技术相比，具有如下优点：本发明可以在无需构造任何“虚” 边界条件的情况下，通过去除结构的刚体位移模态，并控制共轭梯度法的指定 最小误差阈值，精确求解结构响应；计算实施步骤简洁，不需要修改目前通用 的有限元计算框架；采用修正共轭梯度法使得结构总体刚度矩阵的稀疏特性能 够得到很好的利用，分步求解使得整体求解过程占用空间少、计算效率高。

附图说明

附图1为求旋转刚体位移模式的说明图；

附图2为实施例的壳单元示意图；

附图3为实施例的位移误差曲线。

具体实施方式

下面结合附图和具体实施方式对本发明的内容做进一步详细说明。

实施例：

本实施例为一个如图2所示的四边形壳单元构成的平面问题，一种材 料。Abaqus单元类型：S4R。单元尺寸为1*1m，厚度为0.01m，材料参数 为：弹性模量E=1.0，泊松比μ=0.3。

在Abaqus中，对上述单元的左端1、3节点施加固定边界条件，右端2、 4节点施加向右的拉力。计算得出单元各节点的支反力和支反力矩，作为本 方法求解的与Abaqus中有约束的静力等效的无约束问题的外载右端项F。

根据壳单元的位移公式及壳的平面应变和弯曲应变为0，写出上述单元 的刚体位移，原理请参照图1所示。假设将结构整体沿x方向平移一个单位位 移，则在物体内不会产生任何应变，而此时结构的刚体位移模式为只在第6j+1个 自由度上有位移值1其它都为0，具体可写为：

h₁=[1 0 0 0 0 0 1 0 0 0 0 0 … 1 0 0 0 0 0]^T    （35）

同理，将结构整体分别沿y方向和z方向平移一个单位位移，构造平动刚 体位移模式h₂和h₃。

假设将结构上的节点m绕z轴旋转一个微小角度θ，此时m点x方向的位 移为：

$\left(\begin{array}{c}u={x_m}^{\prime }-x_m\\ =-\delta \sin \alpha \\ =-r\sin \mathrm{\alpha \theta }\\ =-y_m\theta \end{array}\right)---\left(36\right)$

同理，y方向的位移为：

v=x_mθ    （37）

z方向的位移为零：

w=0    （38）

此时m点的位移模式可写为：

h₆=[-y_mθ x_mθ 0 0 0 1θ]^T    （39）

在上式中约去θ，并将其扩展到整个结构，即结构整体绕z轴旋转一个微小 角度θ，而在物体内不产生任何应变，此时结构的刚体位移模式可写为：

h₆=[-y₁ x₁ 0 0 0 1 -y₂ x₂ 0 0 0 1 … -y_n/6 x_n/6 0 0 0 1]^T    （40）

其中，x₁～x_n/6是节点1到节点n/6的初始x轴坐标，y₁～y_n/6是节点1到节点n/6 的初始y轴坐标，1≤m≤(n/6)；

同理，将结构整体分别绕x轴和y轴旋转一个微小角度θ，构造转动刚体位 移模式h₄和h₅，最后本实施例的刚体位移模式为：

$\left(\begin{array}{c}H=\left(\begin{array}{cccccc}{h}_1& h_2& h_3& h_4& h_5& h_6\end{array}\right)\\ =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& {-y}_1\\ 0& 1& 0& 0& 0& x_1\\ 0& 0& 1& y_1& -x_1& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·& ·\\ 1& 0& 0& 0& 0& -y_4\\ 0& 1& 0& 0& 0& x_4\\ 0& 0& 1& y_4& -x_4& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\end{array}\right)---\left(41\right)$

其中，x₁～x₄是节点1到节点4的初始x轴坐标，y₁～y₄是节点1到节点4的初 始y轴坐标；

首先将平动刚体位移模式h₁单位化，且与单位向量e₁（e₁为第一个元素为1， 其它元素都为0的单位向量）一起用来构造Householder矩阵P₁，使得：

P₁*h₁=e₁    (42)

然后，构造Householder矩阵P₂，值得注意的是，Householder矩阵P₂的构 造方式与Householder矩阵P₁略有不同，构成新的特征向量h₂′，并将h₂′的第一 个元素赋值为零，即

h₂′=P₁*h₂    （43）

h₂′(1)=0    （44）

将新构造的平动刚体位移模式h₂′单位化与单位向量e₂一起构造 Householder矩阵P₂：

P₂*h₂′=e₂    （45）

依次类推，根据构造Householder矩阵P₂的方法逐个对平动刚体位移模式h₃ 以及转动刚体位移模式h₄～h₆进行操作，构造新的平动刚体位移模式h₃′以及新 的转动刚体位移模式h₄′～h₆′，进一步构造Householder矩阵P₃～P₆。

利用以上所构造的6个Householder矩阵去除结构的刚体模态，即用6个 Householder矩阵对结构原始刚度矩阵做正交相似变换，得到去除刚体模态后的 结构的刚度矩阵K_p：

K_p=P^TKP    (46)

其中，K为原始的结构刚度矩阵；P=P₁P₂P₃P₄P₅P₆。

根据上式生成的去除刚体模态后的结构的刚度矩阵K_p以及从Abaqus中 生成的右端项F（为与Abaqus计算结果对比，采用Abaqus生成的外载F， 包含支反力、支反力矩与施加的外力），做结构的静力分析：

KU=F    (47)

其中，F为施加在结构上的力，即外载右端项；U为结构在力F下产生的位移。

将F与K_p代入上式：

K_PU_P=F_P    (48)

其中，U_P=P^TU；F_P=P^TF。

采用修正的共轭梯度法对上式进行求解，得到位移：

U=P*U_P    (49)

由于K_p的前6行6列都为0，并强制使F_P的前6个值为0，从而不求解式 （15）中的前6个方程，这样一来，U_P的前6个值可以为任意值，进而求解出 来的位移U就会包含刚体位移，最后只需减去其刚体位移，就可与真实值进行 比较，即：

U_n=U-k_ih_i(i=1、2…6)    (50)

其中，h_i为第i个刚体位移模式，k_i为系数。

将上式（50）求得的减去刚体位移后的位移与Abaqus求解的位移结果 对比见表1所示。

表1位移求解结果与Abaqus求解的结果对比

把Abaqus求解的结果作为位移真值，去刚体位移的位移误差见图3所 示。由此可见对于壳单元，本方法的结果是相当精确的，完全可以满足工程 计算需要。

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