单元刚度矩阵MATLAB编程.docx

资源描述

单元刚度矩阵MATLAB编程.docx

《单元刚度矩阵MATLAB编程.docx》由会员分享，可在线阅读，更多相关《单元刚度矩阵MATLAB编程.docx（16页珍藏版）》请在冰豆网上搜索。

单元刚度矩阵MATLAB编程.docx

单元刚度矩阵MATLAB编程

《有限元法》实验报告

专业班级力学（实验）1601

姓名田诗豪

学号

提交日期

实验编号

实验一

实验二

实验三

总分

得分

实验一（30分）

一、实验内容

编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB

函数文件[B3,S3,K3]=ele_mat_tri3（xy3,mat），其中：

输入变量xy3为结点坐标数组，mat为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。

（要求给

出3个不同算例进行验证，并绘制出单元形状和结点号）

、程序代码

通用函数

function[B3,S3,K3]=ele_mat_tri3（xy3,mat）

%fe成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数

%xy3结点坐标数组

%mat材料参数矩阵（弹性模量，泊松比，壁厚）

%B3应变矩阵

%S3应力矩阵

%K3单元刚度矩阵

°%*********************************xyh=[1,xy3（1,1）,xy3（1,2）;1,xy3（2,1）,xy3（2,2）;1,xy3（3,1）,xy3（3,2）];

A=*det（xyh）;

A=abs（A）;

D=mat

（1）/（1-mat

（2）A2）*[1,mat

（2）,0;mat

（2）,1,0;0,0,（1-mat

（2））/2];

b=zeros（1,3）;c=zeros（1,3）;

°%*********************************

fori=1:

ifi==1

j=2;

m=3;

elseifi==2

j=3;

m=1;

else

j=1;

m=2;

end

b（i）=xy3（j,2）-xy3（m,2）;

c（i）=xy3（m,1）-xy3（j,1）;

end

°%*********************************

B31=1/（2*A）*[b

（1）,0;0,c

（1）;c

（1）,b

（1）];

B32=1/（2*A）*[b

（2）,0;0,c

（2）;c

（2）,b

（2）];

B33=1/（2*A）*[b（3）,0;0,c（3）;c（3）,b（3）];

B3=[B31,B32,B33];

^%*********************************

S3=D*B3;

^%*********************************

K3=A*mat（3）*B3'*D*B3;

主程序

clear;clc;

xy3=[0,0;5,1;1,4];

mat=[3e6,,];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

[B3,S3,K3]=ele_mat_tri3（xy3,mat）

三、算例分析

算例1:

如图1所示三角形单元，结点坐标为1（0,0），2（5,2），3（1,4），弹性模量为

200GPa泊松比为、厚度为。

试求应变矩阵，应力矩阵和单元刚度矩阵。

图1算例1三角形单元

解：

根据如图1所示三角形单元及其几何和材料参数，编制主程序如下:

xy3=[0,0;5,2;1,4];

mat=[2e11,,];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

[B3,S3,K3]=ele_mat_tri3（xy3,mat）

运行程序，得到应变矩阵B3如下:

得到应力矩阵S3（Pa）如下:

+10

+09

+10

+09

+10

+09

+10

+09

+10

+09

得到单元刚度矩阵K3（Pa）如下:

+10

+09

+10

+09

+10

+09

+10

+09

+10

+09

+10

+09

+10

算例2:

如图2所示三角形单元，

结点坐标为1（0,0）

，2（3,0），3（0,5）

，弹性模量为

200GPa泊松比为、厚度为。

试求应变矩阵，应力矩阵和单元刚度矩阵。

图2算例2三角形单元

解：

根据如图2所示三角形单元及其几何和材料参数，编制主程序如下:

xy3=[0,0;3,0;5,0];

mat=[2e11,,];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

[B3,S3,K3]=ele_mat_tri3（xy3,mat）

运行程序，得到应变矩阵B3如下:

得到应力矩阵S3（Pa）如下:

+10

+00

+10

+00

+10

+00

+10

+00

得到单元刚度矩阵

K3（Pa）

如下：

+11

+10

+00

+10

+00

+10

+00

+10

+00

+10

+00

+10

+00

+10

算例3:

如图3所示三角形单元，结点坐标为1（0,0）,2（3,0）,3,丁3,弹性模量为200GPa泊松比为、厚度为。

试求应变矩阵，应力矩阵和单元刚度矩阵。

图3算例3三角形单元

解：

根据如图3所示三角形单元及其几何和材料参数，编制主程序如下:

clear;clc;

xy3=[0,0;3,0;,*sqrt（3）];

mat=[2e11,,];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

[B3,S3,K3]=ele_mat_tri3（xy3,mat）

运行程序，得到应变矩阵B3如下:

得到应力矩阵S3（Pa）如下:

+10

+00

+10

+00

+10

+00

得到单元刚度矩阵

K3（Pa）

如下：

+10

+08

+10

+08

+10

+08

+10

+08

+10

+00

+10

+00

+10

实验二（30分）

一、实验内容

编写一个计算平面4结点四边形等参元的刚度矩阵的MATLAB®数文件K4=ele_mat_quad4（xy4,mat），其中：

输入变量xy4为结点坐标数组，mat为材料参数矩阵；输出变量K4为单元刚度矩阵。

（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）

二、程序代码

通用函数

functionK4=ele_mat_quad4（xy4,mat）

毗成平面4结点四边形等参元的刚度矩阵的功能函数

主程序

xy4=[3,2;8,3;7,8;4,7];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

mat=[3e6,,];

K4=ele_mat_quad4（xy4,mat）

三、算例分析

算例1:

如图5所示四边形等参单元（图4为局部坐标系规则单元），已知4个结点整

体坐标系内的坐标为1（3,2）;2（8,3）;3（7,8）;4（4,7）。

弹性模量为200GPa泊松比

为、厚度为。

试求单元刚度矩阵。

图4局部坐标系规则单元

图5算例1四边形等参单元

解：

根据如图5所示四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

^%*********

输入结点坐标数组********

xy4=[3,2;8,3;7,8;4,7];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****mat=[2e11,,];

K4=ele_mat_quad4（xy4,mat）

运行程序，得到单元刚度矩阵

K4（Pa）如下

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+07

+09

+08

+09

+08

+09

+08

+09

+08

+09

+08

+09

+07

+08

+09

算例2:

如图

6所示四边形等参单元，已知

4个结点整体坐标系内的坐标为

1（1,1）;

2（4,1）;3（4,3）;4（1,4）。

弹性模量为180GPa泊松比为、厚度为。

试求单元刚度矩

阵。

图6算例2四边形等参单元

解：

根据如图6所示四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

^%*********

输入结点坐标数组********

xy4=[1,1;4,1;4,3;1,3];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚广***mat=[,,];

K4=ele_mat_quad4（xy4,mat）

运行程序，得到单元刚度矩阵K4（Pa）如下:

+09

+08

+09

+10

+09

+10

+09

+08

+09

+10

+09

+10

+09

+08

+09

+10

+09

+10

+09

+08

+09

+10

+09

+10

算例3:

如图

7所示四边形等参单兀，已知

4个结点整体坐标系内的坐标为

1（1,1）;

2（4,1）;3（5,3）;4（2,3）。

弹性模量为180GPa泊松比为、厚度为。

试求单元刚度矩

阵。

图7算例3四边形等参单元

解：

根据如图7所示四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

^%*********

输入结点坐标数组********

xy4=[1,1;4,1;5,3;2,3];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****mat=[,,];

K4=ele_mat_quad4（xy4,mat）

运行程序，

得到单元刚度矩阵

K4（Pa）如下：

+09

+08

+09

+08

+09

+08

+10

+09

+08

+10

+09

+08

+09

+10

+09

+10

+09

+10

+09

+08

+09

+08

+10

+09

+08

+09

+08

+09

+08

+10

+09

+10

+09

+10

+09

+10

实验三（40分）

1、实验内容

编写一个计算平面8结点四边形等参元刚度矩阵的MATLAB函数文件K8=ele_mat_quad8（xy8,mat），其中：

输入变量xy8为结点坐标数组，mat为材料参数矩阵；输出变量K8为单元刚度矩阵。

（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）

2、程序代码

通用函数

functionK8=ele_mat_quad8（xy8,mat）

毗成平面8结点四边形等参元刚度矩阵的功能函数

%xy8结点坐标数组

%mat材料参数矩阵（弹性模量，泊松比，壁厚）

%K8单元刚度矩阵

%D弹性矩阵

°%*********************************

PN=sym（zeros（2,8））;

D=mat

（1）/（1-mat

（2）A2）*[1,mat

（2）,0;mat

（2）,1,0;0,0,（1-mat

（2））/2];symsy1y2real

（1）=（1-y1）*（1-y2）*（-1-y1-y2）/4;N

（2）=（1+y1）*（1-y2）*（-1-y2+y1）/4;

N（3）=（1+y1）*（1+y2）*（-1+y2+y1）/4;N（4）=（1-y1）*（1+y2）*（-1+y2-y1）/4;

N（5）=（1-y1A2）*（1-y2）/2;N（6）=（1-y2A2）*（1+y1）/2;

N（7）=（1-y1A2）*（1+y2）/2;N（8）=（1-y1）*（1-y2A2）/2;

°%*********************************i=1:

PN1（i）=diff（N（i）,y1）;PN2（i）=diff（N（i）,y2）;

fori=1:

PN（:

i）=[PN1（i）,PN2（i）]';end

J=PN*xy8;JN=inv（J）;

J1=JN（1,:

）;J2=JN（2,:

）;

B1=[J1*PN（:

1）,0;0,J2*PN（:

1）;J2*PN（:

1）,J1*PN（:

1）];

B2=[J1*PN（:

2）,0;0,J2*PN（:

2）;J2*PN（:

2）,J1*PN（:

2）];

B3=[J1*PN（:

3）,0;0,J2*PN（:

3）;J2*PN（:

3）,J1*PN（:

3）];

B4=[J1*PN（:

4）,0;0,J2*PN（:

4）;J2*PN（:

4）,J1*PN（:

4）];

B5=[J1*PN（:

5）,0;0,J2*PN（:

5）;J2*PN（:

5）,J1*PN（:

5）];

B6=[J1*PN（:

6）,0;0,J2*PN（:

6）;J2*PN（:

6）,J1*PN（:

6）];

B7=[J1*PN（:

7）,0;0,J2*PN（:

7）;J2*PN（:

7）,J1*PN（:

7）];

B8=[J1*PN（:

8）,0;0,J2*PN（:

8）;J2*PN（:

8）,J1*PN（:

8）];

B=[B1,B2,B3,B4,B5,B6,B7,B8];

°%*********************************

m=mat（3）*B'*D*B*det（J）;

%*********数值积分（Guass,n=4）****

（1）=;C

（2）=;

C（3）=;C（4）=;

（1）=;A

（2）=;

A（3）=;A（4）=;

sum=0;

fori=1:

forj=1:

sum=sum+A（i）*A（j）*subs（m,[y1,y2],[C（i）,C（j）]）;end

end

K8=vpa（sum,9）;

主程序

clc;clear;

xy8=[2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

mat=[3e6,,];

K8=ele_mat_quad8（xy8,mat）

3、算例分析

算例1:

如图9所示曲四边形等参单元（图8为局部坐标系规则单元），已知8个结点整体坐标系内的坐标为1（2,2）;2（7,2）;3（6,7）;4（3,6）;5（5,3）;6（6,4）;7（4,6）;

8（3,4）。

弹性模量为100GPa泊松比为、厚度为。

试求单元刚度矩阵。

图8局部坐标系规则单元

图9算例1曲四边形等参单元

解：

根据如图9所示曲四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

^%*********

输入结点坐标数组********

xy8=[2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

mat=[,,];

K8=ele_mat_quad8（xy8,mat）

1（-5,-4）;

弹性模量为

算例2:

如图10所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为

2（6,-12）;3（10,8）;4（-8,12）;5（0,-6）;6（5,-1）;7（0,4）;8（-4,2）

200GPa泊松比为、厚度为。

试求单元刚度矩阵。

图10算例2曲四边形等参单元

解：

根据如图10所示曲四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

xy8=[-5,-4;6,-12;10,8;-8,12;0,-6;5,-1;0,4;-4,2];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

mat=[2e11,,];

K8=ele_mat_quad8（xy8,mat）

算例3:

如图11所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为1（-6,-4）;

2（6,-4）;3（6,4）;4（-6,4）;5（0,-3）;6（3,0）;7（0,3）;8（-3,0）。

弹性模量为200GPa

泊松比为、厚度为。

试求单元刚度矩阵。

图11算例3曲四边形等参单元

解：

根据如图11所示曲四边形等参单元及其几何和材料参数，编制主程序如下:

clc;clear;

xy8=[-6,-4;6,-4;6,4;-6,4;0,-3;3,0;0,3;-3,0];

%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****

mat=[2e11,,];

K8=ele_mat_quad8（xy8,mat）

总结与备注

这次有限元编程实验给了我很多收获，不仅锻炼了Matlab编程的能力，还加强了Word排版的能力，最重要的是加深了对有限单元法的理解，加深了自己的力学功底。

本次实验的所有算例，不仅都通过我的程序得出结果，还经过宋戎妆同学的程序验证，故可以保证程序与算例的正确性。

在这里衷心地感谢周博教授的指导，同时感谢几位帮助过我的同学。

展开阅读全文