20节点单元有限元的三维应力平衡问题

基于20节点单元有限元的三维应力平衡问题计算

作者：Planetarian@pku.edu.cn

介绍

写了一个三维有限元的程序，可以自动划分单元和节点、单元分析、组集、处理边界条件、计算可视化。但是本文不提供该程序，只是详细介绍了有限元的原理和推导。

本来要做视频的，但写完程序和文档太累了，只想丢个结果不想花心思做视频，于是只想发个文档。

第一部分

一、格点生成与单元划分

实现原理：

1.记录格点位置

对三个方向嵌套三层循环，生成27节点的网格，使用计数器判断处于单元的第几个点，使用oklst和计数器比较，判断该节点是否属于20节点单元的节点。如果是，则将该点的坐标存储于coordinates矩阵中，如果不是，则不记录该点位置，并将该点坐标编号记录在fakeindexlst中，用于后续修正格点编号。

2.将格点编号与单元关联

因为1-20格点在每个单元的位置是固定的，因此其编号在全局格点中有周期性。实现方法是生成27节点的单元网格，然后利用编号的周期性取需要的20节点。值得注意的是，这里的单元内格点编号和书上的20节点编号顺序并不相同，但我的程序与这里的格点编号是对应的。

3.修正格点编号

这样虽然形成了对应关系，但格点编号却是27节点的编号。因此最后一步是把格点编号修正为20节点，具体而言，是减去小于该节点编号的假节点个数，例如20号节点在修正前是第27号节点，减去27之前的7个假节点，修正为第20号。

测试：

1单元20节点的格点

 

Y方向有两个单元，XZ方向有一个单元生成的格点。

二、单元分析

处理数据的方式是先单元分析逐个计算刚度矩阵，再整体处理约束方程。

2.1 问题描述：

三维弹性介质的本构方程

三维弹性介质的几何方程

无重力情况三维弹性介质的平衡方程

所以

2.2 基于弱形式的有限元单元分析解

形函数选取20节点的形函数。

2.3 变换到自然坐标系计算刚度矩阵

导数的坐标变换：

是单元内20个格点在的全局坐标向量。

因此某单元的刚度矩阵为：

2.4 数值积分

使用三个高斯点，高斯点及对应权重。

三、组集

3.1 组集刚度矩阵

组集时建立节点局部自由度索引向量和节点总体自由度索引向量的关系，提取单元的20节点全局格点编号，计算每个单元的节点总体自由度索引向量。对m号节点，其总体自由度索引是。对每个单元的节点总体自由度索引向量，使用双重循环将单元的刚度矩阵的值赋给总体刚度矩阵的对应位置。

四、边界约束的考虑

对刚度矩阵和广义力的修正

记录有边界约束的格点全局自由度索引和位移值，然后将总体刚度矩阵对应位置的索引处改为1，该索引所在行其他元素设为0，将力的对应位置改为位移值。即

五、求解方程

一些技术细节：

避免直接求解逆矩阵，而是使用高斯消去法求解。MATLAB提供了这一功能。

使用稀疏矩阵存储总体刚度矩阵。MATLAB提供了这一功能。

结果展示我怕被人拿去交自己的作业，就不贴上来了。随便show几张图

结束。