平面三角形单元常应变单元matlab程序的编制模板 下载本文

三角形常应变单元程序的编制与使用

有限元法是求解微分方程边值问题的一种通用数值方法,该方法是一种基于变分法(或变分里兹法)而发展起来的求解微分方程的数值计算方法,以计算机为手段,采用分片近似,进而逼近整体的研究思想求解物理问题。

有限元分析的基本步骤可归纳为三大步:结构离散、单元分析和整体分析。 对于平面问题,结构离散常用的网格形状有三角形、矩形、任意四边形,以三个顶点为节点的三角形单元是最简单的平面单元,它较矩形或四边形对曲边边界有更好的适应性,而矩形或四边形单元较三节点三角形有更高的计算精度。

Matlab语言是进行矩阵运算的强大工具,因此,用Matlab语言编写有限元中平面问题的程序有优越性。本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序,程序流程图见图1。

有限元法中三节点三角形分析结构的步骤如下:

1)整理原始数据,如材料性质、荷载条件、约

束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2)单元分析,建立单元刚度矩阵。 3)整体分析,建立总刚矩阵。

4)建立整体结构的等效节点荷载和总荷载矩

5)边界条件处理。 6)解方程,求出节点位移。 7)求出各单元的单元应力。

8)计算结果整理。计算结果整理包括位移和应

力两个方面;位移计算结果一般不需要特别的处理,利用计算出的节点位移分量,就可画出结构任意方向的位移云图;而应力解的 误差表现在单元内部不满足平衡方程,单元 与单元边界处应力一般不连续,在边界上应 力解一般与力的边界条件不相符合。

-

开始 输入初始数据 生成单刚集成总刚 施加约束信息 生成荷载向量 边界条件处理 计算结点位移 计算单元应力 计算结果整理 结束 图1 程序流程图

1

1.1 程序说明

%******************************************************************* % 三角形常应变单元求解结构主程序

%******************************************************************* ? 功能:运用有限元法中三角形常应变单元解平面问题的计算主程序。 ? 基本思想:单元结点按右手法则顺序编号。 ? 荷载类型:可计算结点荷载。

? 说明:主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算

法的全过程,即实现程序流程图的程序表达。

%----------------------------------------------------------------------------------------------------- 1 程序准备

format short e clear all clc

? 说明:

format short e - 设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法; clear all clc

- 清除所有已定义变量,目的是在本程序的运行过程中,不会- 清屏,使屏幕在本程序运行开始时

发生变量名相同等可能使计算出错的情况;

%----------------------------------------------------------------------------------------------------- 2 全局变量定义

global NNODE NPION NELEM NVFIX NFORCE COORD

LNODS YOUNG POISS THICK

global FORCE FIXED

global BMATX DMATX SMATX AREA global ASTIF ASLOD ASDISP

global FP1 ? 说明:

NNODE—单元结点数,NPION—总结点数, NELEM—单元数,NVFIX—受约束边界点数,NFORCE—结点力数,COORD—结构结点坐标数组,LNODS—单元定义数组,YOUNG—弹性模量,POISS—泊松比,THICK—厚度

-

%设定输出类型 %清除所有已定义变量 %清屏

2

FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向; FIXED— 约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0

BMATX—单元应变矩阵(3*6), DMATX—单元弹性矩阵(3*3),SMATX—单元应力矩阵(3*6),AREA—单元面积

ASTIF—总体刚度矩阵,ASLOD—总体荷载向量,ASDISP—结点位移向量

FP1—数据文件指针 3 打开文件

FP1=fopen('input.txt','rt'); ? 说明:

FP1=fopen('input.txt','rt'); - 打开已存在的输入数据文件input.txt,且设置其为只读格式,使程序在执行过程中不能改变输入文件中的数值,并用文件句柄FP1来执行

FP2=fopen('output.txt','wt'); -打开输出数据文件,该文件不存在时,通过此命令创建新文件,该文件存在时则将原有内容全部删除。该文件设置为可写格式,可在程序执行过程中向输出文件写入数据。

%----------------------------------------------------------------------------------------------------- 4 读入程序控制信息

NPION=fscanf(FP1,'%d',1) %结点个数(结点编码总数) NELEM=fscanf(FP1,'%d',1) %单元个数(单元编码总数) NFORCE=fscanf(FP1,'%d',1) %结点荷载个数 NVFIX=fscanf(FP1,'%d',1) %

YOUNG=fscanf(FP1,'%e',1) %弹性模量 POISS=fscanf(FP1,'%f',1) %泊松比 THICK=fscanf(FP1,'%d',1) %厚度

LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组(单元结点号) ? 说明:

建立LNODS矩阵,该矩阵指出了每一单元的连接信息。

矩阵的每一行针对每一单元,共计 NELEM;每一列相应为单元结点号(编码)、按逆时针顺序输入。

命令中,[3,NELEM]’表示读取NELEM行3列数据赋值给LNODS矩阵。 显然,LNODS(i,1:3)依次表示i单元的i,j,k结点号。 COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组 ? 说明:

-

%打开输入数据文件 存放初始数据

3