数字旗手

电气化、自动化、数字化、智能化、智慧化

0%

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

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

桁架元的特点

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

指定整体坐标系为X-Y,局部坐标系为x-y。则两者之间的转换关系为:

即:

其中:

局部坐标系下的有限元方程为:

为了把有限元方程从局部坐标系变换到整体坐标系,可通过转换矩阵:

所以:

又因为转换矩阵T满足如下关系(可实际计算验证一下):

所以:

所以整体坐标系的刚度矩阵与局部坐标系的刚度矩阵关系为:

由于局部坐标系下的单元刚度矩阵为:

其中$k=\frac{EA}{L}$。那么整体坐标系下的单刚为:

其中$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]

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