有限元数值分析课程论文

《有限元数值分析》课程论文要求

1．简述平面问题有限元法（自选单元类型）。

2．写出计算程序框图。

3．自编一份平面问题有限元弹性应力计算程序。

4．以图1为算例，，分析计算结果。

图1

1.简述平面问题有限元法

1.1有限元法的基本思想：

1.1.1假想把连续系统分割成数目有限的单元，单元之间只在数目有限的指定点处（称为节点）相互连接，构成一个单元集合体来代替原来的连续系统。在节点上引进等效载荷，代替实际作用于系统上的外载荷。

1.1.2对每个单元由分块挖的思想，按一定的规则建立求解未知量与节点相互作用之间的关系。

1.1.3把所有单元的这种特性关系按一定的条件集合起来，引入边界条件，构成一组以节点变量为未知量的代数方程组，求解之就得到有限个节点处的待求变量。

1.2平面问题包括平面应力和平面应变问题

1.2.1平面应力问题

平面应力问题研究等厚度薄板状弹性体，厚度尺寸远远小于截面尺寸，

受力方向沿板面方向，不沿厚度变化。

1.2.2平面应变问题

平面应变问题处理面内受力但垂直于平面方向上不产生变形的二维受力问题，沿截面方向的截面形状和大小相同且厚度尺寸远远大于截面尺寸，载荷垂直于厚度方向且沿厚度均匀分布。

1.3结合所给模型说明平面问题的三角形单元求解

1.3.1 结构的离散化

将所给实际模型简化，并划分为如图所示的2个三节点三角形单元，并对其进行编码，节点编号和坐标、单元编号以及相关节点如表1、表2所示。

图2

表1 节点编码及其坐标

表2 单元编码及其相关节点

表3 节点及节点载荷

表4 节点及节点位移

节点载荷列阵(1)

节点位移列阵：（2）1.3.2 分别计算单元的刚度矩阵

将每个单元局部编码，都可看做如图4的形式，每个单元都有其对应的i、j、m。

图3

单元刚度矩阵（t为板料厚度）（3）

为几何矩阵，与节点坐标有关，其计算公式如下：

（4）其中

单元面积（5）

为弹性矩阵，与材料的弹性模量E以及泊松比μ有关，本问题按照平面应力问题求解，平面应力问题的弹性矩阵的表达式如式（7）所示：

（6）其中

在所划分的模型中，这2个单元的面积是相同的。因为本问题是简单的平面应力问题，在求解几何矩阵时，没有用到位移插值函数，而是直接进行几何矩阵的求解。首先求出节点坐标的差和单元面积，得到几何矩阵，并将弹性模量E和泊松比μ带入式（6）得到弹性矩阵，最后根据公式（3）求得单元的刚度矩阵。这里有2个单元，因此将有2个单元刚度矩阵。

1.3.3合成总体刚度矩阵

当单元进行局部编码时（如图4），每个单元的刚度矩阵都有如下形式;

（7）

按单元的位移自由度所对应的位置进行组装可以得到整体刚度矩阵

1.3.4 引入边界条件求解刚度方程

将所得到的总刚度矩阵、节点载荷列阵（1）式以及节点位移列阵（2）式代入整体刚度方程中，总体刚度方程如式（8）。

（8）

将位移边界条件

和力边界位移条件

代入总体刚度方程，对总体刚度矩阵

进行降阶，划去

的对应位移为0的行和列。

总体刚度矩阵方程降阶后为线性方程组，采用牛顿迭代法求解，得到节点位移阵列

{}δ，再将节点位移代入总体刚度矩阵方程中，可求得节点载荷{}F 。

1.3.5 求解单元应变和单元应力

将每个单元的位移和几何矩阵代入式（9），求得每个单元的应变 {}[]{}e

e

e

B δε⋅= （9）

其中{}{}m

m

j

j i i

e

v u v u v u

=δ

将每个单元的应变代入式（10），求得每个单元的应力 {}[]{}e

D εσ⋅= （10）

2、计算程序框图

3、平面应力问题有限元弹性应力计算程序

#include

#include "stdafx.h"

#include

#include "FYMatrixMN.h"

using namespace std;

#define d11(E,V) E/(1-V)/(1+V)

#define d12(E,V) E*V/(1-V)/(1+V)

#define d33(E,V) E/2/(1+V)

typedef FYNAMESPACE::FYMatrixMN FyMatrix;

int main(int argc, char* argv[])

{

cout<<"程序中采用的单位是MPa、mm、N"<

double E=210000, V=0.3,t=2;

double x1=0, y1=0, x2, y2=0, x3, y3, x4=0, y4;

cout<

cin>>x1>>y1>>x2>>y2>>x3>>y3>>x4>>y4;

//计算单元1的单元刚度矩阵,S11是面积

double S11;

FyMatrix danyuan1(6,6);

FyMatrix B1(3,6);

FyMatrix D1(3,3);

FyMatrix H1(3,6);

double x32,x21,x13;

double y12,y23,y31;

x32=x3-x2; x21=x2-x1; x13=x1-x3;

y12=y1-y2; y23=y2-y3; y31=y3-y1;

S11=((x2*y3-x3*y2)-(x1*y3-x3*y1)+(x1*y2-x2*y1))/2.0;

B1.AllZero();

B1(0,0)=y23;

B1(0,2)=y31;

B1(0,4)=y12;

B1(1,1)=x32;

B1(1,3)=x13;

B1(1,5)=x21;

B1(2,0)=x32;

B1(2,1)=y23;

B1(2,2)=x13;

B1(2,3)=y31;

B1(2,4)=x21;

B1(2,5)=y12;

D1.AllZero();

D1(0,0)=d11(E,V);

D1(0,1)=d12(E,V);

D1(1,0)=D1(0,1);

D1(1,1)=D1(0,0);

D1(2,2)=d33(E,V);