注:此处的矩阵位移法使用后处理法
杆件单元可以简称「单元」,本文中有时也简称「构件」
在结构力学中——当材料变形很小时,其「应力」与「应变」的关系接近线性。这个特性被总结为线弹性小变形。
当「节点」发生位移或是转角,都会在「杆件」的连接下,在另一个「节点」上产生应力。(注意,位移和力的关系是线性的)
『矩阵位移法』中的「位移法」的意义是,当每个节点的位移和转角都达到合适的位置时——「节点」在「杆件」约束下产生的应力刚好抵消荷载。而正是因为「线性」,这种力学关系可以用「矩阵」表达出来。
$$ {\begin{bmatrix} \bar{F_{Ni}} \ \bar{F_{Qi}} \ \bar{M_i} \ \bar{F_{Nj}} \ \bar{F_{Qj}} \ \bar{M_j} \ \end{bmatrix}
\begin{bmatrix} \frac{EA}{l} & 0 & 0 & {\color{Red}-}\frac{EA}{l} & 0 & 0\ 0 & \frac{12EI}{l^3} & \frac{6EI}{l^2} & 0 & \frac{{\color{Red}-}12EI}{l^3} & \frac{6EI}{l^2}\ 0 & \frac{6EI}{l^2} & \frac{4EI}{l} & 0 & \frac{{\color{Red}-}6EI}{l^2} & \frac{2EI}{l}\ {\color{Red}-}\frac{EA}{l} & 0 & 0 & \frac{EA}{l} & 0 & 0\ 0 & \frac{{\color{Red}-}12EI}{l^3} & \frac{{\color{Red}-}6EI}{l^2} & 0 & \frac{12EI}{l^3} & \frac{{\color{Red}-}6EI}{l^2}\ 0 & \frac{6EI}{l^2} & \frac{2EI}{l} & 0 & \frac{{\color{Red}-}6EI}{l^2} & \frac{4EI}{l} \end{bmatrix}
\begin{bmatrix} \bar{u_i} \ \bar{v_i} \ \bar{\theta_i} \ \bar{u_j} \ \bar{v_j} \ \bar{\theta_j} \ \end{bmatrix}}\ \mathbf{\bar{F}^e} = \mathbf{\bar{K}^e} \mathbf{\bar{\delta}^e} $$
-
$$u$$ 、$$v$$、$$\theta$$——节点的轴向位移、法向位移、转角(轴向与法向是相对于构件) -
$$F_N$$ 、$$F_Q$$、$M$——节点受构件约束产生的力,分别是轴向力、法向力、弯矩 -
$$_i$$ 、$$_j$$——代表杆件的两端 -
$$\mathbf{\bar{F}^e}$$ 、$$\mathbf{\bar{\delta}^e}$$、$$\mathbf{\bar{K}^e}$$——「单元坐标系」中的杆端力、杆端位移、单元刚度矩阵
「单元刚度矩阵」是描述单个「杆件」的矩阵(因为「杆件」代表了两个节点之间的约束关系、所以单元刚度矩阵也是在表达两个节点之间的关系)
以矩阵的第二行为例,可以得到以下关系 $$ \left. \begin{bmatrix} 0 & \frac{12EI}{l^3} & \frac{6EI}{l^2} & 0 & \frac{{\color{Red}-}12EI}{l^3} & \frac{6EI}{l^2}\end{bmatrix} \right.
\begin{bmatrix} \bar{u_i} \ \bar{v_i} \ \bar{\theta_i} \ \bar{u_j} \ \bar{v_j} \ \bar{\theta_j} \ \end{bmatrix} \ = 0 \cdot \bar{u_i} + \frac{12EI}{l^3} \cdot \bar{v_i}
- \frac{6EI}{l^2} \cdot \bar{\theta_i}+ 0 \cdot \bar{u_j}
- \frac{{\color{Red}-}12EI}{l^3} \cdot \bar{v_j} + \frac{6EI}{l^2} \cdot \bar{\theta_j} = \bar{F_{Qi}} $$ 可以看到$$\overset{i端的剪力}{\bar{F_{Qi}}}$$与$$\overset{i端的轴向位移}{\bar{u_i}}$$、$$\overset{j端的轴向位移}{\bar{u_j}}$$无关,而与$$\overset{i端的法向位移}{\bar{u_i}}$$、$$\overset{i端的转角}{\bar{\theta_i}}$$、$$\overset{j端的法向位移}{v\bar{_i}}$$、$$\overset{j端的转角}{\bar{\theta_j}}$$相关。
$$ \mathbf{T} = \left.\begin{bmatrix} \cos \alpha & \sin \alpha & 0 & 0 & 0 & 0\ -\sin \alpha & \cos \alpha & 0 & 0 & 0 & 0\ 0 & 0 & 1 & 0 & 0 & 0\ 0 & 0 & 0 & \cos \alpha & \sin \alpha & 0\ 0 & 0 & 0 & -\sin \alpha & \cos \alpha & 0\ 0 & 0 & 0 & 0 & 0 & 1\ \end{bmatrix}\right.\
令 \quad \mathbf{\bar{F}^e} = \mathbf{T}\mathbf{{F}^e} \qquad \mathbf{\bar{\delta}^e} = \mathbf{T}\mathbf{{\delta}^e} \qquad \mathbf{\bar{K}^e}=\mathbf{T} \mathbf{K^e} \mathbf{T}^T\
\Rightarrow \mathbf{{F}^e} = \mathbf{T}^T\mathbf{\bar{F}^e} \qquad \mathbf{{\delta}^e} = \mathbf{T}^T\mathbf{\bar{\delta}^e} \qquad \mathbf{{K}^e}=\mathbf{T}^T \mathbf{\bar{K}^e} \mathbf{T}\
\Rightarrow \mathbf{{F}^e} = \mathbf{K^e} \mathbf{{\delta}^e} $$
-
$$\alpha$$ ——杆件与$$x$$轴的夹角 -
$$\mathbf{T}$$ ——坐标转换矩阵是一个正交矩阵($$\mathbf{T}^T = \mathbf{T}^{-1}$$) -
$$\mathbf{{F}^e}$$ 、$$\mathbf{{\delta}^e}$$、$$\mathbf{{K}^e}$$——「结构坐标系」中的杆端力、杆端位移、单元刚度矩阵
因为杆件不一定平行于$x$轴,所以需要利用坐标转换矩阵。
单元定位向量$$\mathbf{\lambda^e}$$
==todo==
==todo==
本项目基于矩阵位移法(后处理法)进行内力计算。
矩阵位移法是求弹性构件内力精确解的方法。
本项目使用后处理法,与先处理法相比
- 优点
- 编码方便
- 稍微更好理解一些
- 缺点
- 资料少,同济的《结构力学》没教
- 通用性差(但我用了一些特殊的方法修复了通用性)
矩阵位移法的核心是总刚度矩阵,在后处理法中,总钢矩阵的长和宽都是节点数的三倍。每一行每一列各代表节点的一个分量($\nu, \mu, \theta$)。
单元是附加在节点上的联系,这种联系用单元刚度矩阵描述(两个节点之间相对位移产生的力),代表着两个节点之间位移与力之间的关系。在总钢矩阵中,单元刚度矩阵会叠加在单元对应的两个节点上(因为不同单元的对节点的作用是叠加的)。
组装完成的总钢矩阵代表着所有节点位移与力的关系,求解该矩阵与荷载,其实就是求解一个状态——所有节点都在合适的位置,在这个位置上拉扯单元产生的力恰好与荷载相同。
-
『相比弯矩分配法』弯矩分配法不考虑水平和垂直位移的分配,在非对称结构/荷载下与矩阵位移法算出的结构误差较大(可能达到10%)。对所有节点都增加水平和垂直位移的约束(考虑到对轴力图的影响以及杆件轴向变形比较小,垂直约束可以不用加),与弯矩分配法的数值差可以减少两个左右的数量级。
-
『相比D值法』
-
D值法的抗侧移刚度是通过本层梁柱进行计算的。而矩阵位移法中,柱子受到上层框架和地面的(转角)约束越大,抗侧移刚度更大。越接近顶层抗侧移刚度越小。
-
如果只在某两层之间施加剪力,在D值法中只有这两层之间存在层间水平位移。但在矩阵位移法中,因为转角的变形协调,会使得其他层也产生层间水平位移(可以达到存在剪力层的30%)。
在风荷载计算中,用矩阵位移法计算时,因为其他层的影响,层间水平位移会比仅本层有剪力的情况大100%左右。但是最终算出的水平位移和D值法大多偏差只有5%。(换而言之,D值法的抗侧移刚度会考虑其他层剪力对本层的影响。但这个考虑是静态的,不会随着实际其他楼层情况而改变(不会实际计算其他层对本层抗侧移刚度的影响)。)(但首层(D值法偏大)和顶层(D值法偏小)会差30%~40%)。
-
虽然抗侧移刚度没有考虑其他层对本层端部转角约束的影响,但是在内力计算时,反弯点的位置会随着层数/相邻层高比而变化。所以内力图比侧移量精确得多。
-
-
『坐标手性』将坐标系从左手坐标系转到数学常用的右手坐标系不需要修改坐标转换矩阵/单元刚度矩阵等...因为力和位移的坐标系同时被转换了。
-
『通用性增强』一般后处理法只能处理两端刚节点的单元。但是只需要根据两端约束情况写适配的单元刚度矩阵,即可适配所有情况。(虽然两端6个分量总共有64种情况,但是轴向与法向和转角无关,分出来为4*16种情况。这16种情况中,只有5种是超静定。(静定结构不因位移产生内力))
-
『边界条件』边界条件使用置零法,将无意义分量对应行列(总钢矩阵)设为零,将主对角线值设为1,对应的荷载分量设为0。
-
原本想写弯矩分配法,但没成功(因为弯矩分配法不适合写通用求解器)
-
Python与Ansys apdl有限元系列二:矩阵位移法计算桁架结构_王.伟的博客-CSDN博客
大概看了一眼,知道实现矩阵位移法的代码量不大,但因为程序的结构和功能差异过大,无法借鉴。
-
这篇很有用,特别是关于边界条件的处理
-
一个基于矩阵位移法的结构力学计算器,我编写的过程中,数据就是与这个计算器做对比。 (我当时不知道这个结构力学计算器可以精确编辑EA/EI和杆件长度(不过建模大框架时,还是比较麻烦),所以我自己写了一个)
(这个结构计算器在计算定向链杆连接时似乎会出错)
plot.py
是绘图模块,利用旋转矩阵和偏移向量将函数将求解后的应力曲线绘制到合适的位置。solver.py
是求解模块,定义了Node
、Element
、Struction
三个主要的类,分别定义节点
、杆件单元
和整体结构
。utils.py
是工具模块,放置一些通用的工具函数。