利用Mathematica进行有限元编程(二):桁架元分析

本文是对Mathematica有限元分析与工程应用一书的学习笔记。

桁架元的特点

平面桁架元是既有局部坐标又有整体坐标的二维有限元,因此比起之前的杆单元,需要多一步坐标变换。
桁架元示意图如下:

指定整体坐标系为X-Y,局部坐标系为x-y。则两者之间的转换关系为:
$$
\begin{equation}
\begin{split}
U_{Xi}&=u_{xi}cos\theta,U_{Yi}=u_{xi}sin\theta \\
U_{Xj}&=u_{xj}cos\theta,U_{Yj}=u_{xj}sin\theta
\end{split}
\end{equation}
$$
即:
$$
U=Tu
$$
其中:
$$
U=
\begin{bmatrix}
U_{Xi} \\
U_{Yi} \\
U_{Xj} \\
U_{Yj}
\end{bmatrix}
,T=
\begin{bmatrix}
cos\theta & 0 \\
sin\theta & 0 \\
0 & cos\theta \\
0 & sin\theta
\end{bmatrix}
,u=
\begin{bmatrix}
u_{xi} \\
u_{xj}
\end{bmatrix}
$$
局部坐标系下的有限元方程为:
$$
f=Ku
$$
为了把有限元方程从局部坐标系变换到整体坐标系,可通过转换矩阵:
$$
u=T^{-1}U, f=T^{-1}F
$$
所以:
$$
F=TKT^{-1}U
$$
又因为转换矩阵T满足如下关系(可实际计算验证一下):
$$
T^{T}T=I
$$
所以:
$$
F=TKT^TU=\overline{K}U
$$
所以整体坐标系的刚度矩阵与局部坐标系的刚度矩阵关系为:
$$
\overline{K}=TKT^T
$$
由于局部坐标系下的单元刚度矩阵为:
$$
K=
\begin{bmatrix}
k & -k \\
-k & k
\end{bmatrix}
$$
其中$k=\frac{EA}{L}$。那么整体坐标系下的单刚为:
$$
\overline{K}=\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}
$$
其中$C=cos\theta,S=sin\theta$。

模块分析:

建立单元刚度矩阵(经过了坐标变换)

1
2
3
4
5
6
7
8
9
TrussElementKm[EE_, AA_, LL_, theta_] := Module[{},
x = theta*Pi/180;
w1 = Cos[x]^2;
w2 = Sin[x]^2;
w3 = Sin[x]*Cos[x];
y = EE*AA/
LL*{{w1, w3, -w1, -w3}, {w3, w2, -w3, -w2}, {-w1, -w3, w1,
w3}, {-w3, -w2, w3, w2}};
y]

组装整体刚度矩阵

1
2
3
4
5
6
7
8
AssembleSpringKm[p1_, p2_, m_] := Module[{j, k}, f = {p1, p2};
For[j = 1, j <= 2, j++, For[k = 1, k <= 2, k++,
GlobalK[[2 f[[j]], 2 f[[k]]]] += m[[2 j, 2 k]];
GlobalK[[2 f[[j]] - 1, 2 f[[k]]]] += m[[2 j - 1, 2 k]];
GlobalK[[2 f[[j]], 2 f[[k]] - 1]] += m[[2 j, 2 k - 1]];
GlobalK[[2 f[[j]] - 1, 2 f[[k]] - 1]] += m[[2 j - 1, 2 k - 1]];
]];
GlobalK]

这里的组装与之前的杆单元不同,注意此处每个节点上有两个自由度,但总体原则还是将单刚的每个元素叠加到总刚的对应位置上,只是自由度的多少决定了每个单刚的矩阵块的大小,所以得乘以适当的系数。
比如平面刚架元,其既考虑轴向变形,也考虑弯曲变形,每个节点上有三个自由度,其总刚组装时的系数同时变化,如图: