# 2. 基于直接刚度法的杆系有限元方法

## 2.1 弹簧的力学分析原理

首先是大家都知道的胡克定律： $$F=k\delta$$。 但是这里引入了一个新的概念：刚度方程/矩阵。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-cf754424d27e394e9767c548d0901c3b14bffce9%2Fvlcsnap-2021-05-03-15h09m52s741.png?alt=media)

乍一看是画蛇添足之举，实则另有用处，这么处理之后的单元规范了，就可以和其他杆件连接了。

下一步分析两个弹簧的连接+存在一个外力$$F\_2$$作用在节点2上的实例：

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-80ceb6ce4b54f0e93be99f55326e7967d07f13ea%2Fvlcsnap-2021-05-03-15h14m51s011.png?alt=media)

建模过程视频很详细，略了，大致思路是：$$F\_2$$ 给右边的弹簧，弹簧内力设置为$$F\_{I\_2}$$，左右相互抵消，这个$$F\_{I\_2}$$将来要在方程组中约去。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-04a55ad0f36eab629f8768ed30e4d4ffaa458c3d%2Fvlcsnap-2021-05-03-19h42m04s254.png?alt=media)

这里扩充矩阵，使其标准化为3✖3的矩阵，方便相加，这里取了一个巧，正常我们有两个方程组，例如：$$A\_1x=b\_1$$和$$A\_2x=b\_2$$，要揉成一个，一般会取两个系数，比如$$\alpha\_1$$和$$\alpha\_2$$，加上系数以后相加：

$$(\alpha\_1A\_1 + \alpha\_2A\_2)x=(\alpha\_1b\_1 + \alpha\_2b\_2)$$

然后调整这些系数保证$$F\_{I\_2}$$约去，视频里直接让两个方程组b的第二项的内力互为相反数，这样直接加起来就约掉了。

最后考虑$$u\_1$$和$$u\_3$$都是0（约束条件），所以得到u2的计算公式：

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-dfb6f4984d6f14eac91bd2fe79308ae07d371495%2Fvlcsnap-2021-05-03-19h50m13s064.png?alt=media)

## 2.2 弹簧单元与杆单元的比较

杆单元和弹簧一样，都是线性元件，所以计算方式也是一样的。考虑弹性模量为E，长度为l，截面面积A的杆，等效于 $$k = \frac{EA}{l}$$ 的弹簧。

Q.E.D.

## 2.3 杆单元的坐标变换

之前的杆单元的位移u都是沿着杆的方向定义的，但是这个时候如果出现斜拉的杆件，计算就会显得有点复杂。

2.3节的主要思想，就是统一坐标系，全部分解到x-y或者x-y-z坐标系，这样刚度矩阵加起来就更加的流畅。

这一切，都会导致计算量大幅上升，但是这是值得的，因为带来的好处就是计算更加规范了，这样就利于编程实现，反正最后受苦的是CPU/GPU又不是我。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-1079a2cb1c0de5f367e7304e0736b699de0238b8%2Fvlcsnap-2021-05-03-21h33m39s221.png?alt=media)

这里就很清晰的展示了坐标变换的方式，虽然2个坐标变4个了，但是考虑小变形问题的话，$$\alpha$$ 是基本不会变化的，计算中省略其变化，所以秩其实还是2。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-42251bdc9d82c2656af5187824fa579ca4cb4d72%2Fvlcsnap-2021-05-03-21h35m54s041.png?alt=media)

这里我们进行矩阵化，我们用矩阵T联系了变化前和变化后的坐标。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-6973cdca608b5b383ecccb7d2ac7a328961aa31e%2Fvlcsnap-2021-05-03-21h36m26s415.png?alt=media)

我们也需要对刚度矩阵进行同等的变换才可以让刚度阵（其实还有F向量）适应新的坐标，直接思考刚度阵的变化是很费脑细胞的，所以这里使用（线性）代数进行等效变换，我们得到了新的刚度矩阵，过程如下：

* 首先把坐标$$\bold{q}^e$$替换成$$\bold{T}^e\bold{\overline{q}}^e$$
* 两边左侧乘以矩阵$$\bold{T}^{eT}$$。
  * 其实根据线性代数的理论这里其实乘什么都行，但是由于要对力F也进行坐标变换才乘以的这个矩阵。
  * 只有用$$\bold{T}^{eT}$$，右侧的式子才有物理意义了，代表的是力。
* 最后根据线性代数的结合律，我们把$$\bold{T}^{eT}\bold{K}^e\bold{T}^e$$捆绑一下，就得到了新的刚度矩阵。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-dcda4ecab9c4e17ac599784ac226bee0618d7b3a%2Fvlcsnap-2021-05-03-21h50m37s065.png?alt=media)

老师还列举了三维的形式，三维的刚度矩阵写出来可能一张PPT就满了（6x6），所以就略了。

## 2.4 一个四杆结构的实例分析

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-a04cf20e0d43c227992ca1da85ff5511b5520c2f%2Fvlcsnap-2021-05-03-23h29m03s501.png?alt=media)

以上是今天要分析的问题，四个线性元件组成的桁架结构的受力分析。第一件事情其实是体力活儿，列出所有杆和节点的信息：

| 杆子 | 节点    | l   | $$n\_x$$ | $$n\_y$$ |
| -- | ----- | --- | -------- | -------- |
| 1  | {1,2} | 400 | 1        | 0        |
| 2  | {3,2} | 300 | 0        | -1       |
| 3  | {1,3} | 500 | 0.8      | 0.6      |
| 4  | {4,3} | 400 | 1        | 0        |

| 节点 | x   | y   |
| -- | --- | --- |
| 1  | 0   | 0   |
| 2  | 400 | 0   |
| 3  | 400 | 300 |
| 4  | 0   | 300 |

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-da892bbc162bbfba2f2eaf369bbc1a8769463bb2%2Fvlcsnap-2021-05-03-23h38m14s659.png?alt=media)

随后就是看着图片和表格填刚度矩阵，这个活太细碎了……

曾老师是先填写了一个个小的矩阵，然后扩充填0以后相加，不过我觉得直接填写应该也OK。

这里的支反力是因为有了约束条件，当你约束一个坐标的时候，坐标是已知了，但相应的你就产生一个未知的力。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-aa87b75dc0a67cda725356b97aa1a58921edabc5%2Fvlcsnap-2021-05-03-23h38m40s776.png?alt=media)

这就是刚度矩阵，好特么大……（对于人类而言），其中的$$F\_x2$$，$$F\_x3$$，$$F\_y3$$都是刚才已知的力。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-57b02a9e2dab6f94a83e1d9026847271732f63c4%2Fvlcsnap-2021-05-03-23h40m00s293.png?alt=media)

考虑到$$u\_1$$，$$v\_1$$，$$v\_2$$，$$u\_4$$，$$v\_4$$全都是0，所以未知数只有3个，而且$$u\_2$$直接就可以肉眼算出，最后只剩下一个二元线性方程组。

![](https://3247607006-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-M6TU8XK0hd3CZdz_gQ_%2Fuploads%2Fgit-blob-bd160b4c8ba1fe50385959d89fadc0e701adcb6f%2Fvlcsnap-2021-05-03-23h40m13s974.png?alt=media)

节点位移求出以后，支反力也就知道了，这个是很自然的，因为位移知道以后，拉伸量也就知道了。

## 2.5 四杆结构的ANSYS实例分析

略