三维有限元编程教学 从入门到开门 （FORTRAN）

小编是水工结构研究生
最近前前后后花了一个暑假的时间完成一个六面体八节点单元的有限元计算程序，这里打算分几个部分来写下我的心得和前前后后的总结，希望给有需要的小伙伴一些提示和建议。

这里首先说明几点建议（与其说是建议不如说是心得）吧：

1.首先你需要对有限元方法有一定的认识，先从简单的基础知识学起《数值分析》，《有限元原理及应用》陈老写的这本书可谓是一本敲门砖。一定要先看书再开始编程，研究了里面的数学知识就好办了许多。
2.虽然学习编程的方法，但是面对长达近千行的代码还是需要有耐心，需要掌握调试debug的技巧，要会设置断点，耐心的看里面的每一步的计算成果。当然对于vs编译器，对于工科生真的很不友好，跳出的问题几乎看不懂，都得去网上搜。（吐槽一下）
3.自己编程后对有限元的计算的过程会理解很多，这也是我刚开始编程的目的，当然程序还有不足的地方我也还在修改，这里就先写出来希望有小伙伴看到给我指正。

一、整体结构

准备文件 控制变量 （文章整体思路）
数学知识补充 刚度矩阵
数学知识补充 求解
数学知识补充 三种荷载
结果应力输出

对于有限元编程的一个整体思路是，事先将结构的网格画好，编号节点顺序（按照顺序走不要跳，不然会影响），按照准备文件的格式准备；程序读入数据，程序生成刚度矩阵（从形函数，到雅克比矩阵，弹性矩阵D，应变矩阵B，高斯积分点）；各个单元刚度矩阵组装；形成荷载数组（包括自重，面力，体力三种荷载的生成与组装）；分解迭代求解；位移数据后处理。
下面介绍一下主要的控制变量（因为小编是按照陈老教材上的二维程序改成的三维所以变量和子程序都差不多这里就截图表示吧，不一样的地方后续程序里都会注释说明）

二、准备文件

这里展示准备文件的样式

文件只准备了一个单元类型，并且荷载类型也只有有一个自重荷载（这个是小编检查程序用的准备文件所以很简单，可惜最后计算结果有点问题，但是但是但是，我觉得大致思路是对的，希望有老板看到可以给我指正一下哪里出现了小毛病）
接下来我来我来详细的给我大家介绍一下准备文件里的几个部分：

1. 总结点数NP、总单元数NE、材料数目NM、约束节点数NR

2. 每个节点的xyz坐标

3. 节点序号、单元材料类型、单元的各个节点编号（一定要按照顺序写，顺时针或者逆时针，右手螺旋定则）

4. 每个节点的约束情况（0代表自由，1代表约束）

5. 序号、材料的弹性模量、泊松比、密度、厚度（三维空间单元用不到）

6. 集中力、自重、面力（重头戏来了） 0代表没有，1代表存在

   如果集中力为1则：需要输入：
   节点号、x方向的力PX、y方向的力PY、Z方向的力PZ

   如果面力是1则，输入：
   序号、单元数目（有那几个单元在这个面上收到力了）、1（静水压）/2（均布荷载）、液体密度/法相压强大小、Z轴坐标（只在静水压下才会读取）、单元面号（因为有六个面，所以不同的面用不同的序号代替，荷载那里会细说）
   单元序号（对应上一行的单元数目下的单元序号）

三、文件读取（fortran）

小伙伴们等着急了吧，来这里上高汤：

SUBROUTINE INPUT (JR,COOR,AE,MEL)
    DIMENSION COOR(3,*),JR(3,*),AE(4,*),MEL(9,*)
    COMMON /CMN1/ NP,NE,NM,NR
    COMMON /CMN2/ N,MX,NH
    

    DO 70 I= 1,NP