3.8 平面三角形单元有限元的MATLAB程序

3.8.1 程序总体功能设定及主程序函数

1．程序总体功能设定

根据有限元法分析力学问题的主要步骤（见图1.3），有限元程序至少具有“数据读入”“形成结构总刚度矩阵”“形成不同工况总载荷矢量”“求解有限元方程”“计算节点位移与单元应力”等五项基本功能。为了更有效便捷地利用应用程序来解决实际问题，程序还应具备文件管理、模型显示、模型数据校验、变形状态及应力分布的图形显示及后处理等方面的附加功能。

（1）数据文件管理 对有限元模型数据文件、计算结果存储文件进行查询与管理，实现快速方便地调用原始数据文件，以防止在保存计算结果文件时无意中覆盖已存在的同名文件。

（2）数据校验功能 程序自动完成数据类型、数量匹配的检查，并反馈错误信息；图形显示模型网格、位移约束、载荷等内容，以辅助人工直观判断数值及有限元模型正误。

（3）后处理功能 图形显示变形图、内力图、应力云图等。

基于上述总体功能设定，将每类问题程序划分成8个功能模块，各功能模块的顺序如图3.8所示。对于初学编程的读者，建议重点演练基本功能函数，可将模型数据作为“Data_Model”函数，自行定义输出文件名。

2．平面线性三角形单元有限元程序功能

（1）预定的程序功能 采用线性三角形单元分析平面应力或平面应变问题的变形及应力时，外载荷包括节点集中力、自重与体力、面载荷、非零位移载荷等四种形式，其中体力可以有多组，面载荷为线性分布，面力作用形式则包括边界的法线方向、切线方向和倾斜方向三种。

程序通过数据文件输入模型数据，便于修改，具有较好的自动校验与模型显示辅助检查功能，计算结果以不同格式保存到与输入文件相对应的路径及文件名，显示可调整变形比例的变形图以及应力分布状态云图的后处理功能。

（2）平面线性三角形单元主函数 参照图3.8中的有限元分析程序的功能模块图，表3.5列出了平面三角形单元的主功能函数Plane_Triangular_Element调用的8个函数的一览表。

（3）平面三角形单元主程序函数Plane_Triangular_Element

1．function_Plane_Triangular_Element

2．％本程序采用线性三角形单元求解平面问题，计算在自重等体力、集中力以及

3．％线性分布面力作用下的变形和应力，并将结果存储到文件

4．％调用以下功能函数完成：读入有限元模型数据、模型图形显示、计算结构总刚、载荷矢量、

5．％求解有限元方程、应力分析、位移应力后处理等功能

提醒一点，函数式程序在调用其他函数时会进行数据交换，可用实参量与虚参量进行数据交换，也可将读入的变量设为全局变量。使用全局变量虽然表达形式简便，但编程时多个功能函数程序段中的变量符号应一致，还要事先将变量符号规划好，修改变量符号工作量较大，多人共同编写程序时更要注意。采用实参量与虚参量格式，编写程序时可灵活自由地选用变量符号；在调用功能函数时，需要使实参量与虚参量的数目、位置、格式保持一致。

MATLAB语言区分字符大小写。为突出程序可读性，本书编写的MATLAB程序，除了常用符号如弹性模量用大写E外，一般矩阵变量用大写首字母表示，简单变量用小写字母表示。

3.8.2 文件管理函数

1．文件管理函数预设功能

“文件管理函数”用于查找有限元模型数据输入文件、生成用于保存计算结果的输出数据文件名，其功能要求如下：

①查找并显示已保存在预定路径下的有限元模型数据文件，数据文件可以是文本（*.txt）或Excel表（*.xls或*.xlsx）类型；②系统会自动在模型数据文件的同路径下，生成原名与“RES”组合的文本格式（.txt）文件或Excel文件，并以两种格式保存有限元计算结果；③检测拟保存计算结果的文件名是否存在，以免误删或覆盖有用的数据文件，若文件名已存在，则提醒选择“覆盖或新建”输出文件名或文件路径。

实现上述功能有两种编程方式：字符变量法和图形用户界面（GUI）。字符变量法，将原始数据文件名定义为字符串变量，并通过键盘输入。若未查找到与“字符串变量”完全吻合的文件，说明输入的文件名不存在或路径错误，提醒重新输入直到正确为止。文件名称，即“字符串变量”必须包含原始数据文件的存储路径、文件名及扩展名，即使是同一文件名的不同版本的扩展名（如*.xls或*.xlsx），系统也会提醒输入错误。该方式编写的程序使用起来极不方便，文件名字符越长，输入出错概率越大。

图形用户界面（GUI）通过uiget-file函数显示相应的窗口，如图3.9所示，单击鼠标查找并选择文件图标，为了减少显示范围，快速浏览并选择有限元模型数据文件，可设定只优先显示文件夹及*.txt、*.xls或*.xlsx等格式的文件。

当检测到拟保存计算结果的文件名已存在时，通过询问函数questdlg，弹出如图3.10所示的询问对话框，通过鼠标选择是否覆盖。新建文件名时，可利用uiputfile函数的显示保存文件的对话框，用鼠标选择保存文件的路径，在对话框内输入文件名，输入文件名时可以不带扩展名。MATLAB的图形用户界面可以实现用户与计算机的直接交互操作，避免键盘输入所造成的文件名或路径错误，使用更方便。

2．文件管理函数流程图

文件管理函数程序 （详见链接文件：第3章 平面三角形单元/1．文件名管理函数.txt） GUI流程如图3.11所示。

3.8.3 有限元模型数据输入函数

有限元分析程序的变量赋值主要有键盘输入和数据文件输入两种方式。键盘输入操作简单，使用灵活，但每次运行程序都需要重复输入数据，输入工作量较大且输入错误不易更正，键盘输入法适用于数据量比较小、一次性输入等情况。

采用数据文件输入，数据直观、便于检查和修改，容易保存为不同的文件格式和预定的文件名，可反复调用、多次使用，对于数据量较大的情况优点更为突出，常用数据文件类型为文本（*.txt）或Excel表（*.xls或*.xlsx）等。采用文本（*.txt）格式数据输入，各类数据按约定顺序排列，程序简单，但数据可读性较差。采用Excel文件，可附加文字以增加数据的可读性，甚至可将各类数据填写在Excel任意单元格内。

1．模型数据输入函数预定功能

①从指定的文件中读取有限元分析所需的全部数据；②用图形显示有限元模型网格划分、节点及单元编号、位移约束、各种载荷作用，辅助人工检查原始数据。

2．数据填写顺序及格式

平面三角形单元的模型数据包括问题属性和材料属性及模型规模概况数据、节点坐标、单元数据、节点位移约束、节点集中载荷、体力、面力作用的边线、面力作用的节点等，共8类数据，各类数据填写顺序及格式要求见表3.6。

对表3.6的模型数据补充说明：

1）弹性问题属性和材料属性。问题属性代码，程序规定：“1”表示平面应力问题，“2”表示平面应变问题，其他数字无效。材料属性包括弹性模量、泊松比、重度、厚度4个数值。

2）单元数据为3列或4列。前3列表示单元的节点号，第4列为体力组别编号，比表3.2多1列，这是程序拓展功能，允许各单元体力不同。如果所有单元的体力相同，可略去第4列数据，采用3列数据。

3）关于线性分布面力的描述数据，采用两个矩阵组合表达方式：

其一，描述存在表面力作用边界的数据，为3列矩阵，如表3.7所示，数组 QMX （ mx ，3），其中 mx 表示有面力的单元边界数量，前两列分别为有表面力作用的单元边界两端节点的有序编号，要求物体始终处在从第1节点到第2节点连线的左侧，即逆时针方向包络物体；第三列为面力形式代码。

其二，描述存在表面力作用的节点处面力大小的数据，为3列或4列矩阵。如图3.5所示，面力有三种典型表达形式：法向力——表面力与边界面垂直，如正压力；切向力——表面力在边界面内，与其切线方向平行，如摩擦力；面力与边界不垂直或相切的倾斜状态，用面力在坐标轴上的分量表示。如果仅有法向力、切向力，则可采用3列矩阵；存在倾斜面力或面力混合编号时，则需采用4列矩阵。

描述三种形式面力的数据列数不同，但格式相似，统一采用四列表示方法：第1列为节点编号、第2列为面力类型、第3与第4列为对应节点处面力的大小，如表3.8所示，数组 QMD （ md ，4），其中 md 表示有面力的节点数量。为了程序统一与简便，填写表3.8中面力集度数据时，法向面力、切向面力的第4列填写“0”。当边界上有不同形式面力作用时，可以实现混合编号而不必再对划分类型分别编号，但同一节点或边界上存在多种形式面力时应分别描述，并分别编号和计数。典型面力代码及符号约定列于表3.9。

3．文本文件数据的格式及MATLAB中读取文本文件数据的基本语句

文本文件具有所占字节小、用MATLAB读取数据语句简单、快捷、添加删减数据灵活便利等优点，被很多程序设计者采用，但数据可读性差，各类数据格式、数量及排列顺序必须严格按程序要求填写，人工查询数据时定位困难。

为了增加数据的可读性以及利用循环读取的数据方式来简化程序步骤，建议有限元模式文本文件数据格式，采用“三句”表示法，即各类数据均包括 提示句 、 数据大小 、 数据主体 三部分：

（1） 提示句 用文字或字符串表示以下数据的类型，放在数据开始位置，程序要求只能占1行；（2） 数据大小 表示该类数据的 及 ，行数实质是该类数据的总体数目，列数与表3.6对应。采用两个整型数，应新起一行，不能与“提示句”处在同一行；

（3） 主体数据 模型计算使用数据。建议逐行填写，为便于数据检查，每行的第一列为编号，第二列至最后一列是数据有效部分；每行应写满要求的列数，不能缺位，序号不计位，否则会出现串位现象。

以下为节点坐标数据描述格式：

（ 提示句 ）：节点坐标【2列/3列】

注：斜体部分不含在数据文件中。

采用文本格式数据文件，应注意：

数据类型不能少，每类数据的结构还应一致。对于主体为零的类型数据，必须有“提示句”，第二项“数据大小”填写“00”，第三项“主体数据”不写。

各类数据必须完整填写预定的列数，“0”不能缺少。同行两个数字之间用空格分开，不能用逗号（，）分割，否则容易读取出错。添加或删减主体数据时，第二项“数据大小”应相应变更。

文本格式数据文件中增加了一些辅助查阅的文字或数字，在程序设计时，应将该部分占位无效的额外数据的读取舍弃，只保留用于计算的有效数据 （详见链接文件：第3章 平面三角形单元\2．平面三角形单元模型数据格式案例.txt） 。MATLAB读取文本文件数据的主要语句及其说明如下：

4．可填写在任意区域单元格内的Excel数据格式

这里提出一种新的数据读取方式，该方式的特点：①允许有限元模型中的各类数据没有顺序限制，可填写在Excel表中的任意区域单元格内；②允许插入文字提示、数据类型等内容，以增加数据的可读性；③按个人操作习惯编写数据，同类数据也可有不同的列数及格式。

5．有限元模型数据输入函数程序的流程图

图3.12表示具有数据匹配性检查、补齐默认项功能的有限元模型数据输入的流程图。采用文本数据文件时可以不进行匹配性检查和默认项补齐运算，读取文本数据文件的平面三角形单元模型数据输入函数为Plane_Tri_Model_txt。

3.8.4 显示有限元模型函数

显示有限元模型的目的是便于人工检查判断有限元模型数据的正误。其功能设定：①以多种颜色显示有限元模型单元网格，在单元中心处标注单元号，在节点处标注节点号；②以不同方式显示外载荷，集中力画在作用节点处；③面力以作用边中心为起点并指向面力方向，法向面力与边线垂直，切向面力画在边线上，倾斜面力以坐标分量表示；④位移约束用三角号表示，如图3.13所示；⑤若模型有误，则终止正在执行的功能函数并跳过后续所有功能函数，直接返回到MATLAB命令界面。

显示弹性平面问题有限元模型函数的流程如图3.14所示。

3.8.5 计算结构总刚度矩阵函数

1．功能要求

1）计算单元的刚度矩阵。建议按3.3.3节的步骤计算单元刚度矩阵，调用弹性矩阵函数Elastic_Matrix（pm，E，nv）生成弹性矩阵 （详见链接文件：第3章 平面三角形单元\3．弹性矩阵函数.txt） ，调用应变矩阵函数Plane_B3_Matrix（ie）计算单元的面积与应变矩阵 （详见链接文件：第3章 平面三角形单元\4．应变矩阵函数.txt） 。由于弹性矩阵与问题属性及材料性质有关，与单元无关，因此生成弹性矩阵的运算应放在单元循环体之前。

2）将单元刚度矩阵按节点编号分块后叠加到总刚度矩阵中，生成结构总的刚度矩阵。有的教材采用逐个元素添加赋值法编程，这样操作虽简单直观，但因单刚元素多而会造成程序冗长，单元自由度越多，程序就会越长，逐一赋值容易出错且不易检查。采用循环语句编程，程序简短、方法通用性强、准确无差错，但需要将单刚按节点分块并计算每个子块在总刚度矩阵中的位置。

2．单刚分块形成总刚度矩阵的功能函数

单元刚度矩阵大小是固定的，元素排列与对应的节点号相关联。单元刚度矩阵按节点分块的2×2阶子矩阵 K _rs （ r ， s = i ， j ， m ），在总刚度矩阵中的位置区域为［（2 r -1）：2 r ，（2 s -1）：2 s ］。集成总刚度矩阵，实质上就是将单刚子矩阵 K _rs （ r ， s = i ， j ， m ）累加到总刚与 r 、 s 相对应的［（2 r -1）：2 r ，（2 s -1）：2 s ］的位置上。

为了更好地理解单刚分块形成总刚度矩阵的过程，下面列出的程序是在已经确定单元刚度矩阵的前提下，实现分块与累加功能。三角形单元单刚集成总刚度矩阵函数为Assemble_Stiffness_Matrix（对单元ie）。

3．计算单刚集成结构总刚度矩阵函数的流程图

上述程序实现了单刚分块形成总刚度矩阵的过程，前提是已经确定单元刚度矩阵。在实际程序中，每计算一个单元刚度矩阵之后，将其直接叠加到总刚度矩阵中，计算单刚集成结构总刚度矩阵函数Plane_Tri_Stiff_Matrix （详见链接文件：第3章 平面三角形单元\5．总刚度矩阵函数.txt） 的流程如图3.15所示。

3.8.6 结构总载荷矢量函数

外力包括重力及体力、线性分布面力等效节点载荷、节点集中力这3种形式，程序的功能是计算单元等效节点载荷并集成有限元分析计算所需的总载荷矢量。该程序段分3块：

（1）体力 体力包含重力和分布体力两种：重力与材料密度有关，所有单元相同，默认等效载荷方向为纵向、取负值（- y 向）；一般分布体力由坐标分量表示，单元体力可不同。计算体力等效节点载荷的过程对单元循环。

（2）线性分布面力 包括法向面力、切向面力、分量表示的3种形式面力，按表3.9中规定的方向，计算面力等效节点载荷的过程是对有面力作用的边界循环。计算3种面力等效节点载荷的通用程序为Equivalent_Nodal_Force_Surface （详见链接文件：第3章 平面三角形单元\6．面力等效节点载荷的通用函数.txt） 。

（3）集中力 集中力则对集中力个数循环，直接按作用点实际编号叠加到总载荷矢量中。

计算单元等效节点载荷并集成结构总载荷矢量的函数为Plane_Tri_Load_Vector （详见链接文件：第3章 平面三角形单元\7．总载荷矢量函数.txt） ，其流程图如图3.16所示。

3.8.7 求解有限元方程函数

求解有限元方程可选用置“1”法、乘大数法或降阶法，按预定格式，将计算结果保存到有限元模型数据文件同目录下同名或自定义的文本文件（.txt）。拟定输出格式要求：标明问题属性、模型节点数、单元数等概况，每行显示节点编号及位移分量。

（1）置“1”法与乘大数法流程 置“1”法求解有限元方程 （详见链接文件：第3章 平面三角形单元\8．置“1”法求解有限方程函数.txt） 的流程如图3.17所示。乘大数法求解有限元方程的流程与置“1”法基本相同。在功能函数程序中，最关键的步骤是计算并确定每个位移约束所对应的总刚度矩阵主对角线元素的位置。设 BC （ r ，1）表示第 r 个约束位置的节点编号， BC （ r ，2）为位移约束方向代码，则与该约束对应的总刚度矩阵主对角线元素的位置为

i 0= kd ×［ BC （ r ，1）-1］＋ BC （ r ，2）

式中， kd 为节点自由度数，平面问题 kd =2，空间问题 kd =3。

（2）降阶法流程 只保留待定的节点位移作为未知量，压缩了线性方程组阶数，减小求解方程工程量，但打乱了节点顺序，总刚度矩阵及载荷矢量需要重新排列，因此在求解方程之前，要根据位移约束条件，进行有限元方程重组。降阶法的关键是有限元方程重构，解决方案是构成一个与节点位移条件相对应的顺序调整矩阵 T ₀₁ 。矩阵 T ₀₁ 与总刚度矩阵 K 同阶，在单位矩阵基础上通过调整矩阵行的排列而得到， T ₀₁ 是一个正交矩阵。调整后的关系为

式中，带波浪线的字母表示顺序调整后的量。

图3.18表示构建顺序调整矩阵功能函数的流程图，其中 ng 表示位移约束总数， nk 表示总自由度数。

3.8.8 计算应力的功能函数

预设功能：①计算单元应力分量；②计算单元的主应力和Mises应力；③统计各单元应力分量、主应力、Mises应力的极大值、极小值；④采用绕点加权平均法，计算节点的应力分量、主应力和Mises应力；⑤将应力结果以设定格式保存到预定的文件中。

计算应力的功能函数为Plane_Tri_Stress （详见链接文件：第3章 平面三角形单元\9．计算单元应力分量函数.txt） ，调用计算应变矩阵函数Plane_B3_Matrix、弹性矩阵函数Elas-tic_Matrix、计算主应力及Mises应力函数Main_Stress （详见链接文件：第3章 平面三角形单元\10．计算主应力函数.txt） 等3个函数，应力分析函数程序流程图如图3.19所示。

3.8.9 平面三角形单元后处理函数

后处理函数的功能包含显示变形图、显示应力分布云图两部分。

（1）显示变形图 显示未变形网格，根据图幅尺寸及变形量初步确定变形的放大系数，显示变形后网格，变形放大系数可调整。

（2）显示应力分布云图 显示不同应力分量、主应力、Mises应力的分布云图，图中有标明应力大小的颜色参考刻度条。

平面三角形单元后处理函数的流程图如图3.20所示。