有限元计算结构力学大作业.docx

有限元计算结构力学大作业.docx

《有限元计算结构力学大作业.docx》由会员分享，可在线阅读，更多相关《有限元计算结构力学大作业.docx（22页珍藏版）》请在冰豆网上搜索。

有限元计算结构力学大作业

-标准化文件发布号：

（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

有限元-计算结构力学-大作业

SHANGHAIJIAOTONGUNIVERSITY

平面应力问题解的Matlab实现

姓名:

heiya168

学号:

帆哥

班级：

指导老师：

1绪论

有限元方法（finiteelementmethod），是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。

将它用于在科学研究中，可成为探究物质客观规律的先进手段。

将它应用于工程技术中，可成为工程设计和分析的可靠工具。

弹性体在载荷作用下，其基本方程可写成以下的三类方程和两种边界条件。

平衡方程——应力与外载荷的关系；几何方程——应变位移关系；物理方程——应力应变关系；力的边界条件；几何边界条件。

应用最小位能原理，并利用上述关系，最终建立由刚度方程，节点位移和等效节点载荷所构成的求解方程。

带入边界条件求解方程，就可以得出弹性力学问题的一般性解答。

本次大作业基于有限元方法的基本原理，使用Matlab这一平台，针对平面应力问题，采用四节点四边形单元编写了求解单元节点位移的程序。

主要内容包括：

1）介绍有限元的基本原理；2）编程基本思路及流程介绍；3）程序原理及说明；4）具体算例这四个部分。

2平面问题的四节点四边形单元

2.1单元的构造

（1）单元的几何和节点描述

平面4节点矩形单元如图4-6所示，单元的节点位移共有8个自由度。

节点的编号为1、2、3、4，各自的位置坐标为（xi,yi）,i=1,2,3,4，各个节点的位移（分别沿x方向和y方向）为（ui,vi）,i=1,2,3,4。

图2.1平面4节点矩形单元

若采用无量纲坐标

（2.1）

则单元4个节点的几何位置为

（2.2）

将所有节点上的位移组成一个列阵，记作

；同样，将所有节点上的各个力也组成一个列阵，记作

，那么

（2.3）

若该单元承受分布外载，可以将其等效到节点上，也可以表示为如式（2.3）所示的节点力。

利用函数插值、几何方程、物理方程以及势能计算公式，可以将单元的所有力学参量用节点位移列阵及相关的插值函数

来表示；下面进行具体的推导。

（2）单元位移场的表达

从图2-1可以看出，节点条件共有8个，即x方向4个

，y方向4个

，因此，x和y方向的位移场可以各有4个待定系数，即取以下多项式作为单元的位移场模式

（2.4）

它们是具有完全一次项的非完全二次项，以上两式中右端的第四项是考虑到x方向和y方向的对称性而取的，除此外xy项还有个重要特点，就是“双线性”，当x或y不变时，沿y或x方向位移函数呈线性变化，这与前面的线性项最为相容，而

或

项是二次曲线变化的。

因此，未选

或

项。

由节点条件，在x=xi，y=yi处，有

（2.5）

将式（2.5）代入式（2.4）中，可以求解出待定系数a0，…，a3和b0，…，b3，然后代回式（4-52）中，经整理后有

（2.6）

其中

（2.7）

如以无量纲坐标系（2.1）来表达，式（2.7）可以写成

（2.8）

将式（2.6）写成矩阵形式，有

（2.9）

其中，

为该单元的形状函数矩阵。

（3）单元应变场的表达

由弹性力学平面问题的几何方程（矩阵形式），有单元应变的表达

（2.10）

其中几何矩阵B（x,y）为

=

（2.11）

式（4-59）中的子矩阵

为

（2.12）

（4）单元应力场的表达

由弹性力学中平面问题的物理方程，可得到单元的应力表达式

（2.13）

（5）单元应力场的表达

以上已将单元的三大基本变量

用基于节点位移列阵来进行表达，见式（2.9）、式（2.10）及式（2.13）；将其代入单元的势能表达式中，有

（2.14）

其中

是4节点矩形单元的刚度矩阵。

将单元的势能对节点位移

取一阶极值，可得到单元的刚度方程

（2.15）

2.2等参变换

由前面的单元构造过程可以看出，一个单元的关键就是计算它的刚度矩阵，而由刚度矩阵的构成可知要实现两个坐标系中单元刚度矩阵的变换，必须计算两个坐标系之间的三种映射关系：

坐标映射

（2.16）

偏导数映射

（2.17）

面积映射

（2.18）

图2.2矩形单元映射为任意四边形单元

（1）两个坐标系之间的函数映射

设如图4-17所示的两个坐标系的坐标映射关系为

（2.19）

x和y方向上可以分别写出各包含有4个待定系数的多项式，即

（2.20）

其中待定系数a0，…，a3和b0，…，b3可由节点映射条件（2.19）来唯一确定。

对照前面4节点矩形单元的单元位移函数式（2.5），映射函数式（2.20）具有完全相同的形式，同样，将求出的待定系数再代回式（2.20）中，重写该式为

（2.21）

其中

（2.22）

（2.23）

这就可以实现两个坐标系间的映射。

（2）两个坐标系之间的偏导数映射

对物理坐标系（x,y）中的任意一个函数Φ（x,y），求它的偏导数，有

（2.24）

则偏导数的变换关系为

（2.25）

写成矩阵形式，有

（2.26）

其中

（2.27）

两个坐标系的偏导数映射关系

（2.28）

（3）两个坐标系之间的面积元映射

如图2-2所示，在物理坐标系（x,y）中，由dξ和dη所围成的微小平行四边形，其面积为

（2.29）

由于dξ和dη在物理坐标系（x,y）中的分量为

（2.30）

其中i和j分别为物理坐标系（x,y）中的x方向和y方向的单位向量。

由（2.29）式，则有面积微元的变换计算

（2.31）

2.3边界条件的处理——置“1”法

设边界条件为第r个自由度的位移为零，即

；可置整体刚度矩阵中所对应对角元素位置的

，而该行和该列的其它元素为零，即

，同时也置对应的载荷元素则

，则

（2.32）

进行以上设置后，这时方程（2.32）应等价于原方程加上了边界条件

，下面考察这种等价性，就（2.32）中的第r行，有

（2.33）

由于置

，则有

即为所要求的位移边界条件。

而式（5-53）中除第r行外，其它各行在对应于r列的位置上都置了“0”，这相当于考虑了

的影响，除此之外其余各项的影响不变；这恰好就是原方程加上了边界条件

的影响。

3有限元分析流程

3.1程序原理和流程

该程序的特点如下：

问题类型：

可用于弹性力学平面应力问题和平面应变问题的分析

单元类型：

四节点四边形单元

材料性质：

单一的均匀的弹性材料

节点信息：

节点信息文件node.txt

单元信息：

单元信息文件element.txt

载荷信息：

等效节点力文件force.txt

约束信息：

节点约束信息文件constrain.txt

结果文件：

输出节点位移文件node_displace.txt

程序流程图如下：

图3.1程序流程图

3.2使用的函数

Stiffness（E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp）：

计算单元的刚度矩阵（具体说明见附录）

Assembly（KK,k,i,j,m）：

进行单元刚度矩阵的组装（具体说明见附录）

3.3文件管理

主程序源文件：

FEM.m

输入参数文件：

node.txt（节点信息文件）

element.txt（单元信息文件）

constrain.txt（位移边界条件文件）

force.txt（载荷状况文件）

输出位移文件：

node_displace.txt（节点位移文件）

3.4数据文件格式

4算例——开方孔的矩形板拉伸分析

4.1问题的具体参数与载荷

（a）对边拉伸的带孔矩形薄板（b）1/4模型的单元划分

图4-1受均匀荷载的作用的带孔矩形薄板薄板及有限元分析模型

如图4-1（a）所示的矩形带孔薄板对边受均匀荷载的作用，该结构在边界上受正向分布拉力0.3MPa。

若取板厚0.02m，弹性模量2.1E11，泊松比μ=0.3，采用FEM.m程序，和ANSYS分别分析该问题，最后给出各节点位移值。

4.2Matlab程序计算

（1）结构的离散化与编号

该薄板的荷载和几何形状关于x轴和轴对称，故可只取结构的1/4作为计算模型。

将此模型化分为四个全等的直角三角形单元，单元编号和节点编号如图4-1（b）所示。

（2）输入节点信息node.txt

（3）输入单元信息element.txt

（4）输入载荷信息force.txt

（5）输入载荷信息constrain.txt

（6）运行主函数FEM.m，输入E=2.1E11，NU-0.3，h=0.02。

（7）输出各节点位移

4.3ANSYS建模计算

选用plane42板单元，建立如下图所示模型，在2、3节点加载水平方向的载荷10000N.

图4-2受均匀荷载的作用的带孔矩形薄板薄板有限元分析模型

分析计算所得的最后结果为：

4.4误差分析

节点位移

2X

2Y

3X

3Y

4X

4Y

2SUM

3SUM

4SUM

Matlab

7.65E-06

2.53E-07

7.33E-06

-1.43E-06

1.87E-06

2.025E-07

7.66E-06

7.469E-06

1.88E-06

ANSYS

7.54E-06

-8.00E-07

7.81E-06

-2.58E-06

1.68E-06

2.45E-07

7.58E-06

8.22E-06

1.70E-06

误差

1.02%

9.18%

9.94%

表4-1Matlab程序与Ansys结果对比表

误差范围在10%以内，还是比较大的，导致误差的原因可能有以下几个方面：

1）刚度矩阵计算时，使用的积分方法不同，本程序中使用的是两点高斯积分。

2）刚度矩阵求解是，使用的算法不同，本程序使用的是高斯消去法。

3）约束处理时使用的方法不同，本程序使用的是置一法。

5总结

用了两周的时间，完成了这次大作业。

在做这次大作业之前，我对于有限元算法的了解还很肤浅，只是模模糊糊的能够记得刚度阵建立的过程，矩阵表示的公式，对于其中具体的物理意义，各个量之间相互关系的表示不太清楚，在做完这次大作业后，对于四边形单元刚度阵的建立过程有了清晰深入的认识和掌握。

关于等参变化的认识，是这次大作业中的难点之一。

对于和位移差值函数相同的坐标差值函数，如何理解其中“坐标差值”物理意义，在差值完成，刚度矩阵建立后，刚度矩阵的位移模式是局部坐标系中的，还是整体坐标系中的位移，这个问题思考了很久，也和很多同学进行了争论和探讨。

最终得到了一致的结论，收获很多。

在完成大作业之前，我对于Matlab编程还是一窍不通，在编程的过程中，很多好的想法无法用计算机语言来实现是一件很痛苦的事情，但是在同学的帮助下和参考资料的借鉴下，完成了这次编程。

其中，刚度矩阵和组装矩阵是自己独立完成的，主程序借鉴了清华大学曾攀教授编写的《有限元分析基础教程》中平面三节点单元的主程序。

写这个总结的时候，也就是大作业完成的时候，其中还有很多不足的地方有待改进，比如两点高斯积分精确度较差的问题，比如无法直接