2D Truss单元计算实例

问题描述

如果所示为 $xy$ 平面内的桁架结构，结构几何信息如图所示。节点 1 处为固定铰支座，节点 4 处为滑动铰支座，节点 5、6、7 处分别受到 $-y$ 方向 $P=100000N$ 的集中力作用；结构中各杆件采用相同的规格，其中弹性模量 $E=200000MPa$，截面面积 $A=4532mm^2$。

编程思路

单元坐标（7个节点）：

单元（11个单元）：

结构自由度（9*2=14）

由于节点 1 处为固定铰支座，节点 4 处为滑动铰支座，所以约束自由度为 $[1,2,8]$。

初始化整体位移向量（数组大小为自由度14）：

$$ displacement = \left \{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \right \} $$

初始化整体载荷向量（数组大小为自由度14）：

$$ force = \left \{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \right \} $$

节点 5、6、7 处分别受到 $-y$ 方向 $P=100000N$ 的集中力作用，所以有：

$$ force = \left \{ 0, 0, 0, 0, 0, 0, 0, 0, 0, -100000, 0, -100000, 0, -100000 \right \} $$

定义结构总体刚度矩阵（$14\times14$）：

$$ stiffiness\_sum = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\end{bmatrix} $$

遍历单元，求解单元刚度矩阵：

对于单元1 = 【节点1，节点2】，节点1的坐标为(-4500,0)，节点2的坐标为(-1500,0)：

$$ \Delta x =(-1500) - (-4500) =3000 \quad \Delta y = 0 - 0 = 0 $$

$$ L = \sqrt{\Delta x * \Delta x + \Delta y * \Delta y} = 3000 $$

$$ \cos \theta = \frac{\Delta x}{L} = 1 \qquad \sin \theta = \frac{\Delta y}{L} = 0 $$

单元刚度矩阵：

$$ \left [ K \right ]=\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix} = \frac{EA}{L}\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix} $$

节点对应自由度： $eleDof = \{0, 1, 2, 3\}$

将单元刚度矩阵组装到整体刚度矩阵：

for(int i = 0; i < 4; ++i){
    for(int j = 0; j < 4; ++j){
        stiffness_sum[eleDof[i], eleDof[j]] += [K][i, j]
    }
}

$$ stiffiness\_sum =\begin{bmatrix} \frac{EA}{3000} & 0 & -\frac{EA}{3000} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{EA}{3000} & 0 & \frac{EA}{3000} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\end{bmatrix} $$

其它单元类似有（注意 $L$ 不全是 3000）：

$$ [K]_2 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix} $$

$$ [K]_3 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix} $$

$$ [K]_4 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix} $$

$$ [K]_5 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix} $$

$$ [K]_6 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & 0.4 & -0.2 & -0.4 \\ 0.4 & 0.8 & -0.4 & -0.8 \\ -0.2 & -0.4 & 0.2 & 0.4 \\ -0.4 & -0.8 & 0.4 & 0.8\end{bmatrix} $$

$$ [K]_7 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & 0.4 & -0.2 & -0.4 \\ 0.4 & 0.8 & -0.4 & -0.8 \\ -0.2 & -0.4 & 0.2 & 0.4 \\ -0.4 & -0.8 & 0.4 & 0.8\end{bmatrix} $$

$$ [K]_8 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & 0.4 & -0.2 & -0.4 \\ 0.4 & 0.8 & -0.4 & -0.8 \\ -0.2 & -0.4 & 0.2 & 0.4 \\ -0.4 & -0.8 & 0.4 & 0.8\end{bmatrix} $$

$$ [K]_9 =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & -0.4 & -0.2 & 0.4 \\ -0.4 & 0.8 & 0.4 & -0.8 \\ -0.2 & 0.4 & 0.2 & -0.4 \\ 0.4 & -0.8 & -0.4 & 0.8\end{bmatrix} $$

$$ [K]_{10} =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & -0.4 & -0.2 & 0.4 \\ -0.4 & 0.8 & 0.4 & -0.8 \\ -0.2 & 0.4 & 0.2 & -0.4 \\ 0.4 & -0.8 & -0.4 & 0.8\end{bmatrix} $$

$$ [K]_{11} =\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix}=\frac{EA}{L}\begin{bmatrix}0.2 & -0.4 & -0.2 & 0.4 \\ -0.4 & 0.8 & 0.4 & -0.8 \\ -0.2 & 0.4 & 0.2 & -0.4 \\ 0.4 & -0.8 & -0.4 & 0.8\end{bmatrix} $$

单元对应的自由度：

$eleDof_2 = \{2，3，4，5\}$、$eleDof_3 = \{4，5，6，7\}$、$eleDof_4 = \{8，9，10，11\}$、$eleDof_5 = \{10，11，12，13\}$、$eleDof_6 = \{0，1，8，9\}$、$eleDof_7 = \{2，3，10，11\}$、$eleDof_8 = \{4，5，12，13\}$、$eleDof_9 = \{2，3，8，9\}$、$eleDof_{10} = \{4，5，10，11\}$、$eleDof_{11} = \{6，7，12，13\}$

将单元刚度矩阵组装到整体刚度矩阵：

$$ stiffiness\_sum =\begin{bmatrix} 356181 & 108095 & -302133 & 0 & 0 & 0 & 0 & 0 & -54047.3& -108095 & 0 &0 &0 &0 \\ 108095 & 216189 & 0 & 0 & 0 & 0 & 0 & 0 & -108095 & -216189 & 0 &0 &0 &0 \\ -302133 & 0 & 712361 & 0 & -302133 & 0 & 0 & 0 & -54047.3& 108095 & -54047.3 &-108095 &0 &0 \\ 0 & 0 & 0 & 432378 & 0 & 0 & 0 & 0 & 108095 & -216189 & -108095 &-216189 &0 &0 \\ 0 & 0 & -302133 & 0 & 712361 & 0 & -302133 & 0 & 0 & 0 & -54047.3 &108095 &-54047.3 &-108095 \\ 0 & 0 & 0 & 0 & 0 & 432378 & 0 & 0 & 0 & 0 & 108095 &-216189 &-108095 &-216189 \\ 0 & 0 & 0 & 0 & -302133 & 0 & 356181 & -108095 & 0 & 0 & 0 &0 &-54047.3 &108095 \\ 0 & 0 & 0 & 0 & 0 & 0 & -108095 & 216189 & 0 & 0 & 0 &0 &108095 &-216189 \\ -54047.3& -108095 & -54047.3& 108095 & 0 & 0 & 0 & 0 & 410228 & 0 & -302133 &0 &0 &0 \\ -108095 & -216189 & 108095 & -216189 & 0 & 0 & 0 & 0 & 0 & 432378 & 0 &0 &0 &0 \\ 0 & 0 & -54047.3& -108095 & -54047.3& 108095 & 0 & 0 & -302133 & 0 & 712361 &0 &-302133 &0 \\ 0 & 0 & -108095 & -216189 & 108095 & -216189 & 0 & 0 & 0 & 0 & 0 &432378 &0 &0 \\ 0 & 0 & 0 & 0 & -54047.3& -108095 & -54047.3& 108095 & 0 & 0 & -302133 &0 &410228 &0 \\ 0 & 0 & 0 & 0 & -108095 & -216189 & 108095 & -216189 & 0 & 0 & 0 &0 &0 &432378 \end{bmatrix} $$

根据约束自由度为 $[1,2,8]$（节点 1 处为固定铰支座，节点 4 处为滑动铰支座），删除刚度矩阵的第 0，1，7行和列，得到。

$$ stiffiness\_sum =\begin{bmatrix} 712361 & 0 & -302133 & 0 & 0 & -54047.3& 108095 & -54047.3 &-108095 &0 &0 \\ 0 & 432378 & 0 & 0 & 0 & 108095 & -216189 & -108095 &-216189 &0 &0 \\ -302133 & 0 & 712361 & 0 & -302133 & 0 & 0 & -54047.3 &108095 &-54047.3 &-108095 \\ 0 & 0 & 0 & 432378 & 0 & 0 & 0 & 108095 &-216189 &-108095 &-216189 \\ 0 & 0 & -302133 & 0 & 356181 & 0 & 0 & 0 &0 &-54047.3 &108095 \\ -54047.3& 108095 & 0 & 0 & 0 & 410228 & 0 & -302133 &0 &0 &0 \\ 108095 & -216189 & 0 & 0 & 0 & 0 & 432378 & 0 &0 &0 &0 \\ -54047.3& -108095 & -54047.3& 108095 & 0 & -302133 & 0 & 712361 &0 &-302133 &0 \\ -108095 & -216189 & 108095 & -216189 & 0 & 0 & 0 & 0 &432378 &0 &0 \\ 0 & 0 & -54047.3& -108095 & -54047.3& 0 & 0 & -302133 &0 &410228 &0 \\ 0 & 0 & -108095 & -216189 & 108095 & 0 & 0 & 0 &0 &0 &432378 \end{bmatrix} $$

根据 $[F] = [K][S]$ 计算得 $[S] = [F][K]^{-1}$ 即 $displacement = force * stiffiness\_sum^{-1} $ (只计算激活得自由度，即位移矢量和集中力矢量要去除第0，1，7个元素)。

$$ displacement= \begin{bmatrix}0.0 & 0.0 & 0.248 & -1.587 & 0.662 & -1.587 & 0.910 & 0.0 & 0.786 & -1.087 & 0.455 & -1.922 & 0.124 & -1.087\end{bmatrix} $$

单元的内力：

$$ F = EA\varepsilon=EA\frac{u_2-u_1}{L}=EA\begin{bmatrix}-C & -S & C & S\end{bmatrix}\begin{Bmatrix}u_1 \\ v_1 \\ u_2 \\ v_2\end{Bmatrix} $$

遍历单元计算可以得到轴力：

$$ axialForce=[[75000], [125000], [75000], [-100000], [-100000], [-167705], [-55902], [55902], [55902], [-55902], [-167705]] $$