平面四节点等参单元matlab实现

计算力学报告平面四节点等参单元

学生姓名：**

学号：********

一、问题描述及分析

在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。

根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。

二、有限元划分描述

在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：

打开ansys，在单元类型中选择solid->Quad 4 node 182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements 选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。如下图：

图1

三、有限元程序及求解

程序求解使用了matlab语言。具体如下：

程序：

clc

clear

E=2e11; %弹性模量

NU=0.3; %泊松比

t=0.1; %厚度

X=xlsread('D:\data','nodes'); %读取节点坐标

elem=xlsread('D:\data','elements'); %读取单元编号

w=[1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]; %有位移

约束的节点

n=size(X,1); %节点数

m=size(elem,1); %单元数

K=zeros(2*n); %初始总体刚度矩阵

for i=1:m

syms Ks Et x y I1 I2 a b B; %定义可能存在的变量

e=[1,1;-1,1;-1,-1;1,-1];

for j=1:4 %形函数

N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);

end

x=0;y=0;

for j=1:4 %标准母单元映射到真实单元x=x+N(j)*X(elem(i,j),1);

y=y+N(j)*X(elem(i,j),2);

end

J1=jacobian([x;y],[Ks;Et]); %雅克比矩阵及其转置

J=J1';

for j=1:4

I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数

I2=diff(N(j),Et);

C=(J^-1)*[I1;I2];

a=C(1); %形函数对x，y的偏导数

b=C(2);

B(1,2*j-1)=a; %组成B阵

B(1,2*j)=0;

B(2,2*j-1)=0;

B(2,2*j)=b;

B(3,2*j-1)=b;

B(3,2*j)=a;

end

D=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵

k=zeros(8,8);

Kss=[-0.906179,-0.538469,0,0.538469,0.906179]; %5*5高斯积分点

Ett=[-0.906179,-0.538469,0,0.538469,0.906179];

H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数

for j=1:5 %高斯积分求单元刚度阵

for l=1:5

Ks=Kss(j);Et=Ett(l);

B=subs(B);

J=subs(J);

k=k+H(j)*H(l)*B'*D*B*det(J);

end

end

G=zeros(8,2*n); %初始总刚变换矩阵

G(1,2*elem(i,1)-1)=1; %总刚变换矩阵

G(2,2*elem(i,1))=1;

G(3,2*elem(i,2)-1)=1;

G(4,2*elem(i,2))=1;

G(5,2*elem(i,3)-1)=1;

G(6,2*elem(i,3))=1;

G(7,2*elem(i,4)-1)=1;

G(8,2*elem(i,4))=1;

K=K+G'*k*G; %总体刚度矩阵合成

end

KK=K;

b=size(w,1);

for i=1:b

K(2*w(i)-1,2*w(i)-1)=1e20;

K(2*w(i),2*w(i))=1e20;

end %置大数法

f=zeros(2*n,1); %初始载荷矩阵

f(10)=-10e3; %加载荷10kN

U=K\f; %节点位移

for i=1:m %将每个单元各个节点位移集合

u(:,i)=[U(2*elem(i,1)-1);U(2*elem(i,1));U(2*elem(i,2)-1);U(2*elem(i,2));U(2*ele m(i,3)-1);U(2*elem(i,3));U(2*elem(i,4)-1);U(2*elem(i,4))];

end

for i=1:m %求单元应力

syms Ks Et x y I1 I2 a b B;

e=[1,1;-1,1;-1,-1;1,-1];

for j=1:4

N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);

end

x=0;y=0;

for j=1:4

x=x+N(j)*X(elem(i,j),1);

y=y+N(j)*X(elem(i,j),2);