单元刚度矩阵MATLAB编程.docx

资源描述

单元刚度矩阵MATLAB编程.docx

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

单元刚度矩阵MATLAB编程.docx

单元刚度矩阵MATLAB编程

单元刚度矩阵（等参元）MATLAB编程（总16页）

《有限元法》实验报告

专业班级力学（实验）1601

姓名田诗豪

学号10

提交日期

实验编号

实验一

实验二

实验三

总分

得分

实验一（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）

%生成平面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）^2）*[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所示三角形单元及其几何和材料参数，编制主程序如下：

clear;clc;

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

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所示三角形单元及其几何和材料参数，编制主程序如下：

clear;clc;

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

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,

，弹性模量为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------------------结点坐标数组

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

%K4-------------------单元刚度矩阵

%y1y2-----------------局部坐标系结点坐标

%D--------------------弹性矩阵

%*********************************

y1y2=[-1,-1;1,-1;1,1;-1,1];

D=mat

（1）/（1-mat

（2）^2）*[1,mat

（2）,0;mat

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

（2））/2];

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

（1）=;C

（2）=;

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

（1）=;A

（2）=;

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

sum=0;

fori=1:

forj=1:

y1=C;y2=C;

k=1:

%*********************************

PN1

（1）=*（y2（j）-1）;PN2

（1）=*（y1（i）-1）;

PN1

（2）=*（y2（j）-1）;PN2

（2）=*（y1（i）+1）;

PN1（3）=*（y2（j）+1）;PN2（3）=*（y1（i）+1）;

PN1（4）=*（y2（j）+1）;PN2（4）=*（y1（i）-1）;

fork=1:

PN（:

k）=[PN1（k）,PN2（k）]';

end

J=PN*xy4;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）];

B=[B1,B2,B3,B4];

%*********************************

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

sum=sum+A（i）*A（j）*m;

end

K4=vpa（sum,9）;

主程序

clc;clear;

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

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）如下：

+09

+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）^2）*[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-y1^2）*（1-y2）/2;N（6）=（1-y2^2）*（1+y1）/2;

N（7）=（1-y1^2）*（1+y2）/2;N（8）=（1-y1）*（1-y2^2）/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

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）

运行程序，得到单元刚度矩阵K8（GPa）如下：

算例2：

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

展开阅读全文