+09 +08
+09 +09 +10
+09 +09 +09
+10 +09 +09
+09 +09 +09
+10 +09 +09
+09 +09 +09
+09 +09 +10
1(1,1)； 2(4,1);
+09
算例3:如图7所示四边形等参单元，已知4个结点整体坐标系内的坐标为
3(5,3); 4(2,3)。弹性模量为180GPa,泊松比为、厚度为。试求单元刚度矩阵。
3
Q.5
0
2
3 6
图7算例3四边形等参单元
解：根据如图7所示四边形等参单元及其几何和材料参数，编制主程序如下: clc;clear;
%***输入结点坐标数组** xy4=[1,1;4,1;5,3;2,3];
%****输入材料参数矩阵(弹性模量 mat=[,,];
K4=ele_mat_quad4(xy4,mat)
泊松比，壁厚
)****
运行程序， 得到单元刚度矩阵
+09 +09 +09 +09 +09 +09 +08
K4( Pa)
如下：
+09 +09 +08 +09 +09 +09 +09
+09 +09 +08 +09 +09 +09 +08 +10
+09 +08 +10 +09 +08 +08 +09 +09
+09 +09 +09 +10 +09 +10 +09 +10
+09 +09 +08 +10 +09 +09 +08 +09
+08 +08 +09 +09 +09 +08 +10 +09
+09 +10 +09 +10 +09 +09 +09 +10
+09
+09
实验三（40分）
1、 实验内容
编写一个计算平面8结点四边形等参元刚度矩阵的 MATLAB函数文件K8 = ele_mat_quad8(xy8,mat),其中：输入变量xy8为结点坐标数组，mat为材料参数 矩阵；输出变量K8为单元刚度矩阵。(要求给出3个不同算例进行验证，并绘制 出单
元形状和结点号) 2、 程序代码
通用函数
function K8 = 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))/ syms y1 y2 real
N(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:8;
PN1(i)=diff(N(i),y1);PN2(i)=diff(N(i),y2); for i=1:8
PN(:,i)=[PN1(i),PN2(i)]'; end
J=PN*xy8;JN=i nv(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)**** C(1)=; C(2)=;
2];
C(3)=; C(4)=; A(1)=; A(2)=; A(3)=; A(4)=; sum=0; for i = 1:4
for j = 1:4
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曲四边形等参单元
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个结点整体坐标系内的坐标为 解：根据如图9所示曲四边形等参单元及其几何和材料参数，编制主程序如下:
1(-5,-4);
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)
运行程序，得到单元刚度矩阵 K8( GPa)如下:
单元刚度矩阵(等参元)MATLAB编程
