编制三角形单元的有限元程序并用该程序分析一平面问题算例

好的，我会尽力回答你的问题。首先，编写三角形单元的有限元程序需要以下步骤：

1. 首先需要确定问题的几何形状和边界条件，例如平面问题可以是一个矩形区域，或者任意形状的区域。

2. 划分网格，将区域分成若干个三角形单元，可以采用Delaunay三角剖分算法。

3. 确定每个节点的坐标和单元节点编号，建立节点和单元之间的关系。

4. 计算每个单元的刚度矩阵和载荷向量，可以采用数值积分方法进行计算。

5. 组装刚度矩阵和载荷向量，得到整个系统的刚度矩阵和载荷向量。

6. 施加边界条件，例如固定边界条件和力边界条件等。

7. 解线性方程组，得到节点的位移和应力分布。

8. 计算应力或变形的最大值和平均值等。

下面，我将用一个平面问题算例来说明如何使用三角形单元的有限元程序。

假设有一个长方形区域，左边界为固定边界，右边界施加一个水平向右的力，其他边界为自由边界。区域的尺寸为 $L=20$，$H=10$，力的大小为 $F=100$。材料的弹性模量为 $E=200\times 10^9$，泊松比为 $\nu=0.3$。

首先，将区域划分为若干个三角形单元，如下图所示：


节点的坐标和单元节点编号如下表所示：

| Node | x | y |
|------|---|---|
| 1 | 0 | 0 |
| 2 | 20| 0 |
| 3 | 20| 10|
| 4 | 0 | 10|
| 5 | 10| 0 |
| 6 | 20| 5 |
| 7 | 10| 10|
| 8 | 0 | 5 |

| Element | Node 1 | Node 2 | Node 3 |
|---------|--------|--------|--------|
| 1 | 1 | 5 | 8 |
| 2 | 2 | 6 | 3 |
| 3 | 6 | 7 | 3 |
| 4 | 5 | 2 | 6 |
| 5 | 8 | 5 | 7 |
| 6 | 1 | 4 | 8 |
| 7 | 7 | 4 | 3 |
| 8 | 5 | 6 | 7 |

然后，可以使用数值积分方法计算每个单元的刚度矩阵和载荷向量。以第一个单元为例，其刚度矩阵和载荷向量为：

$$
\begin{aligned}
& K_1=\frac{E}{1-\nu^2} \begin{pmatrix}
1 & \nu & 0 \\
\nu & 1 & 0 \\
0 & 0 & \frac{1-\nu}{2}
\end{pmatrix}\frac{1}{4A_1}\begin{pmatrix}
2 & 1 & -1 \\
1 & 2 & -1 \\
-1 & -1 & 2
\end{pmatrix} \\
& f_1=F\frac{H}{2}\begin{pmatrix}
1 \\
1 \\
1
\end{pmatrix}
\end{aligned}
$$

其中，$A_1=25$是第一个单元的面积。

接下来，将所有单元的刚度矩阵和载荷向量组装成整个系统的刚度矩阵和载荷向量。施加边界条件，得到线性方程组：

$$
\left[\begin{pmatrix}
K_{11} & K_{12} \\
K_{21} & K_{22}
\end{pmatrix} \right]
\left[\begin{pmatrix}
u_1 \\
u_2
\end{pmatrix} \right]
=
\left[\begin{pmatrix}
f_1 \\
f_2
\end{pmatrix} \right]
$$

其中，$K_{11}$和$K_{22}$是固定边界上的节点的位移已知条件所对应的子矩阵，$K_{12}$和$K_{21}$是未知节点和已知节点之间的子矩阵。

最后，解线性方程组，得到节点的位移和应力分布。可以计算出右边界上的节点的位移为：

$$
u_2=\frac{FL}{2EH}(3\nu+1)
$$

应力分布如下图所示：


其中，最大应力出现在右下角的单元内，为$5.3\times 10^7$。

以上就是使用三角形单元的有限元程序分析平面问题的方法。希望能对你有所帮助。