一文搞懂静力超单元方法的公式详细推导，结果一致

在前面很多文章中讲到了超单元的应用技巧和方法，其中超单元最重要的价值就是通过缩聚模型自由度来缩短求解时间，接下回顾下超单元的主要理论。

一、超单元方法的概念

图1 超单元模型的应用流程

二、超单元方法的意义

三、超单元的理论方程

（2）

(3)

（4）

(5)

(6)

（7）

即将{Ua}前部分采用 表示，等号右边第一项采用 表示，即方程（7）可以表示如下：

（8）

      (9)

     (10)

将式（10）第一行展开，并乘以 ，可得式（11）。

 (11)

（1）定义 为转换边界

（2）定义 为固定边界位移

（3）定义 为自由边界位移

（4）整体的内部位移为固定位移和自由边界位移之和，即 展开式（10）第二行可得

 (12)

将 代入式（12）

（1）定义 为转换边界（界面）刚度矩阵。

（2）定义 为转换边界（界面）载荷矩阵，通过转换可求得界面自由度位移，如式（13）所示。

(13)

(10)

由此可知，要求解界面点Ua的位移可通过求解内部节点Uo的位移得到，即通过将矩阵的自由度进行缩减，即降维。当求解得到uo后，再通过式（16）可得到界面点的位移Ua结果。

1.4 案例说明

根据下式可知，要求解Uo需要对Kaa求逆，然后再求解得到相应变量。

如：某泊松方程系数矩阵有如形式

Fo=[0;-1.1111;-1.1111;0]，Fa=[-2.2222;-2.2222;-2.2222]

在Matlab中实现如下：

1）Matlab中的矩阵求逆可采用inv函数

2）Matlab中的矩阵转置可采用Koa.'或transpose（Koa）

3）Matlab中的矩阵求和和乘积与常规一样，但是当对两个矩阵进行加减运算时，需要确保这两个矩阵具有相同的维度，即它们的行数和列数必须相同。

4）在Matlab中每一列的元素用空格或逗号分隔，而行与行之间则用分号分隔。

4）首先输入Koo、Koa和Kaa的矩阵，同时求Koa的转置得到Kao，以及Kaa的逆矩阵Kaa_1，如下所示。

Koo=[1 0 0 0;-0.1 -1.4 -0.1 0;0 -0.1 -1.4 -0.1;0 0 0 1]

Koa=[0 0 0;0.8 0.8 0;0 0.8 0.8;0 0 0]

Kaa=[-1.6 0 0;0 -1.6 0;0 0 -1.6]

Kao=Koa.'

Kao=[0 0.8 0 0;0 0.8 0.8 0;0 0.8 0 0]

Kaa_1=inv(Kaa)=[-0.6250 0 0;0 -0.6250 0;0 0 -0.6250]

5）其次输入

Fo=[0;-1.1111;-1.1111;0]，Fa=[-2.2222;-2.2222;-2.2222]

6）根据输入的各参数求解出Kre和Fre，即

Kre=Koo-Koa*Kaa_1*Kao，得到

Kre = [1.0000  0 0 0;-0.1000 -0.6000 0.3000 0;0 0.3000 -0.6000 -0.1000;0 0 0 1.0000]

Fre=Fo-Koa*Kaa_1*Fa，得到

Fre =[0;-3.3333;-3.3333;0]

7）最后根据缩聚公式进行求解得到Uo，即

Kre_1=inv(Kre)，

Uo=Kre_1*Fre

8）根据Uo求解出Ua，即

Ua=Kaa_1*(Fa-Kao*Uo)

9）不采用超单元方法，即直接通过对式（14）求解得到Ua和Uo，即

（14）

Kcc=inv([Koo Koa;Kao Kaa])，得到

再求解{Uo;Ua}，即得到

Ud=Kcc*[Fo;Fa]

即Uo=[0;11.1110;11.1110]，Ua=[6.9444；12.4999；6.9444]

结论：该结果与超单元缩聚后结果一致，即通过对系统刚度矩阵缩聚可以快速得到内部节点位移Uo，进而得到界面节点Ua，即为静力超单元缩聚的由来。

需要更多详情及代码请点击阅读原文。