计算力学(八)总体刚度矩阵和边界条件
总体刚度矩阵
在前面给出的代码中,可以看到在计算插值函数值的时候给单元中的各个节点编了号,同时这些节点在整体中也有一个整体编号.为了计算总应变能,需要将单元应变能相加,但由于\(W^e=\dfrac{1}{2}\mathbf{V}^T_e\mathbf{K}_e\mathbf{V}_e\)是在单元节点编号下得到的,不能直接相加,需要转换到整体编号下再相加.
理论上数学推导应用单元指示矩阵(或布尔矩阵)对单元刚度矩阵进行操作,但在实际程序中,直接对自由度进行操作.
单元刚度矩阵的性质
- 单元刚度矩阵的对角线元素表示要使单元的第\(i\)个节点产生单位位移,而其它节点位移为0时,需在节点\(i\)所施加的力.
- 单元刚度矩阵的非对角线元素\(K_{ij}(i\neq j)\)表示要使单元的第\(j\)个节点产生单位位移,而其它节点位移为0时,需要在第\(i\)个节点所施加的力.
- 单元刚度矩阵是对称(symmetry)的
- 单元刚度矩阵是半正定(positive Semi-definite)的.
- 单元刚度矩阵是奇异(singularity)的.所以在未加位移约束条件以前有限元刚度方程的解不是唯一的.
- 刚度矩阵的任一行(或列)代表一个平衡力系;当节点位移全部为线位移时(即为\(C_0\)型问题),任一行(或列)的代数和应为零.
刚度矩阵的性质
- 对称性
- 奇异性
- 半正定性
- 稀疏矩阵(spars matrix)
- 非零元素显现带状性(banded).
总体刚度矩阵的存储
节点自由度只影响所在单元的其它节点,所以总体刚度矩阵中会有很多0元素,这里采用一维变带宽存储来节省储存空间,其中带宽指从每一行的第一个非零元素到该行最后一个非零元素,其间所有的元素个数,由于矩阵的对称性,实际上只用存储从第一个非零元素到主对角元素即可.
一维变带宽存储虽然更能节省存储空间,但必须定义用于主对角元素定位的辅助数组,下面给出存储思路
- 遍历所有单元,找出各个节点所在单元中与之编号相差最大的节点
- 遍历所有节点的所有自由度,得到其半带宽
- 将所有自由度的半带宽相加,计算主对角元在一维存储中的位置
最后要取到总体刚度矩阵中的值,根据主对角元的位置和自由度插值就可得到.
总体编号的选择也会影响到存储空间,下面是同样的网格划分,但采取了不同的总体节点编号.红色的点代表可能为非零的元素,右图中元素更好地聚集在矩阵的主对角线附近,两张图片中的非零条目数量是相同的.
有限元平衡方程
结构总位能 \[ \Pi=\frac{1}{2}\mathbf{V}^T\mathbf{KV}-\mathbf{V}^T\mathbf{P}\\ \] 其中\(\mathbf{V}\)是节点位移向量,\(\mathbf{K}\)是总体刚度矩阵,\(\mathbf{P}\)是节点等效载荷.根据最小势能原理,得到有限元平衡方程 \[ \mathbf{KV}=\mathbf{P}\\ \]
边界条件的处理
- 处理边界条件的直接法
- 处理边界条件的置“1”法
- 处理边界条件的乘大数法
- 支反力的计算
- 处理耦合边界条件的拉格朗日(Lagrange)乘子法
- 处理耦合边界条件的罚函数法
下面给出我用到的方法.
对于第\(k\)个自由度的位移为\(a\)的情况,将总体刚度矩阵的\((k,k)\)号元素置为1第\(k\)行其他元素置为0,将对应第\(k\)行的等效载荷置为a.通过对应的第\(k\)哥方程可以解出第\(k\)个自由度的位移为\(a\). \[ \begin{bmatrix} K_{11}&\cdots&K_{1k}&\cdots&K_{1n}\\ \vdots&\ddots&\vdots&\ddots&\vdots\\ 0&\cdots&1&\cdots&0\\ \vdots&\ddots&\vdots&\ddots&\vdots\\ K_{n1}&\cdots&K_{nk}&\cdots&K_{nn}\\ \end{bmatrix} \begin{bmatrix} V_1\\ \vdots\\ V_k\\ \vdots\\ V_n\\ \end{bmatrix}= \begin{bmatrix} P_1\\ \vdots\\ a\\ \vdots\\ P_n\\ \end{bmatrix} \] 为了保持刚度矩阵的对称性,将上述方程改写成如下形式 \[ \begin{bmatrix} K_{11}&\cdots&0&\cdots&K_{1n}\\ \vdots&\ddots&\vdots&\ddots&\vdots\\ 0&\cdots&1&\cdots&0\\ \vdots&\ddots&\vdots&\ddots&\vdots\\ K_{n1}&\cdots&0&\cdots&K_{nn}\\ \end{bmatrix} \begin{bmatrix} V_1\\ \vdots\\ V_k\\ \vdots\\ V_n\\ \end{bmatrix}= \begin{bmatrix} P_1-K_{1k}a\\ \vdots\\ a\\ \vdots\\ P_n-K_{nk}a\\ \end{bmatrix} \]
下面给出代码(明天补)
下一篇写施加载荷/求解方程组/应力的计算
#! https://zhuanlan.zhihu.com/p/696564490