利用Mathematica进行有限元编程(三)：三角形单元分析

2023-12-11 09:29:27

三角形单元适用于具有复杂和不规则边界形状的问题，是一种常见的网格离散方式。

局部坐标系

之前的杆单元和桁架元的局部坐标系容易建立，即建在其自身上即可。而三角形单元的局部坐标系显然不能这样建立，其常采用一套无量纲的自然坐标系——面积坐标，如下图所示：

三角形123内部有一任意点P，P与顶点1、2、3组成3个子三角形，每个子三角形的面积与总面积之比记为$L_i$，即P点的面积坐标为$(L_1,L_2,L_3)$。
因为三个坐标相加为1，所以仅有两个独立的自然坐标，所以可以简记为：

注意到：面积坐标还有如下特点，在顶点1时，$xi=1$，其他坐标为0。其他顶点亦同，见上图三个顶点处的坐标值。因此，面积坐标还有形函数的功能。即：

写成矩阵形式为:

在局部坐标系下形函数对$xi,eta$的偏导数为：

代码为：

1 2	ShapeN = {xi, eta, 1 - xi - eta}; DeriveN = {D[ShapeN, xi], D[ShapeN, eta]}

结果为：

1	{{1, 0, -1}, {0, 1, -1}}

偏导数在不同坐标系下的变换：

其中，J为Jacobian矩阵。在上式中，当局部坐标系中明确给定函数$N_1$时，等式的左侧可以求出。同时，当基于局部坐标系给出x和y的显式表达时，基于局部坐标系也可以给出Jacobian矩阵的显式表达。具体过程如下：
根据等参元的性质，基于局部坐标给出的标准形函数与整体坐标的形函数完全形同，则：

所以，Jacobian矩阵为：

具体到这里使用的面积坐标，则有：

所以：

Jacobian矩阵的逆变换为：

而Jacobian行列式则为：

在求解刚度矩阵时，所对应的积分中要用到此行列式。积分过程中的变量和区域需要进行变换：

Jacobian矩阵的逆矩阵和行列式则可以通过Inverse和Det命令求得。

应力应变矩阵

对平面应力问题，应力应变矩阵D为：

对平面应变问题，应力应变矩阵D为：

应变位移矩阵

应变、位移的关系为：

代入位移函数表达式，可得：

其中应变位移矩阵的一部分为：

单元刚度矩阵

根据最小势能原理，求得单元刚度矩阵表达式：

模块分析

建立单元刚度矩阵

GenerateLinearTriangKm[cord_] := Module[{Bmatrix, DeriveN, J, km},
      n1 = xi; n2 = eta; n3 = 1 - xi - eta;
  bdv = {{1, 0, -1}, {0, 1, -1}};(*这里是局部坐标系下的偏导数*)
  J = bdv.cord;(*求解Jacobian矩阵*)
  bdv = Simplify[Inverse[J].bdv];(*转换成整体坐标系下的偏导数*)
  Bmatrix = {{bdv[[1, 1]], 0, bdv[[1, 2]], 0, bdv[[1, 3]], 0}, {0, 
              bdv[[2, 1]], 0, bdv[[2, 2]], 0, bdv[[2, 3]]}, {bdv[[2, 1]], 
              bdv[[1, 1]], bdv[[2, 2]], bdv[[1, 2]], bdv[[2, 3]], bdv[[1, 3]]}};
  Dmatrix = 
     ee/(1 - nu*nu) {{1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu)/2}};
  km = h Integrate[
       Det[J] (Transpose[Bmatrix].Dmatrix.Bmatrix), {xi, 0, 1}, {eta, 0,
                  1 - xi}]
                    ]

组装刚度矩阵

AssembleLinearTriangKm[p1_, p2_, p3_, m_] := 
 Module[{f}, f = {p1, p2, p3};
  For[j = 1, j <= 3, j++, For[k = 1, k <= 3, 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]

这里注意矩阵带不带Matrixform在计算时有很大区别。
注意这里的循环次数为3,是因为每个单元有3个节点。

二次三角形单元

二次三角形单元就是在每条边上还各有一个节点，如图：

具体分析过程跟双线性三角形单元相同，只不过形函数更加复杂。且在组装总刚时循环系数为6。

码农公寓