平面四边形八节点等参元matlab程序

广州大学

《有限元方法与程序设计》

学院： 土木工程学院

专业： 结构工程

姓名： 曾一凡

学号： 21115160**

1

% 平面四边形八节点等参元MATLAB程序 % 变量说明&（2015级——结构工程——曾一凡） % YOUNG POISS THICK % 弹性模量 泊松比 厚度

% NPOIN NELEM NVFIX NFORCE % 总结点数 单元数 约束结点个数 受力结点数 % COORD LNODS FORCE % 结构节点整体坐标数组 单元定义数组 结点力数组 % ALLFORCE FIXED HK DISP

% 总体荷载向量 约束信息数组 总体刚度矩阵 结点位移向量 % 1本程序计算了节点位移和单元中心应力并输出到nonde8out.txt文本里 % 2在第四页举了一个实例 用MATLAB算出结果再用ANSYS算出的结果对比 %======================主程序===================== format short e %设定输出类型 clear %清除内存变量

FP1=fopen('nonde8.txt','rt'); %打开初始数据文件

FP2=fopen('nonde8out.txt','wt'); %打开文件 存放计算结果 NPOIN=fscanf(FP1,'%d',1); %结点数 NELEM=fscanf(FP1,'%d',1); %单元数

NFORCE=fscanf(FP1,'%d',1); %作用荷载的结点个数 NVFIX=fscanf(FP1,'%d',1); %约束数 YOUNG=fscanf(FP1,'%e',1); %弹性模量 POISS=fscanf(FP1,'%f',1); %泊松比 THICK=fscanf(FP1,'%f',1); %厚度

LNODS=fscanf(FP1,'%d',[8,NELEM])'; %单元结点号数组（逆时针） COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号 x,y坐标(整体坐标下) FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点力数组

% 节点力：结点号、X方向力(向右正)，Y方向力(向上正) FIXED=fscanf(FP1,'%d',NVFIX)';

% 约束信息：约束对应的位移编码（共计 NVFIX 组） EK=zeros(2*8,2*8); % 单元刚度矩阵并清零

HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零 X=zeros(1,8); %存放单元8个x方向坐标的向量并清零 Y=zeros(1,8); %存放单元8个y方向坐标的向量并清零 %--------------------------求总刚------------------------------- for i=1:NELEM % 对单元个数循环 for m=1:8% 对单元结点个数循环

X(m)=COORD(LNODS(i,m),1); %单元8个x方向坐标的向量 Y(m)=COORD(LNODS(i,m),2); %单元8个y方向坐标的向量 end

EK=eKe(X,Y,YOUNG,POISS,THICK); %调用单元刚度矩阵

2

a=LNODS(i,:); %临时向量,用来记录当前单元的结点编号 for j=1:8 %对行进行循环---按结点号循环 for k=1:8 %对列进行循环---按结点号循环

HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+... EK(j*2-1:j*2,k*2-1:k*2); % 单刚子块叠加到总刚中 end end end

ALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE); % 调用函数生成荷载向量 %-------------------------处理约束-------------------------------- for j=1:NVFIX % 对约束个数进行循环 N1=FIXED(j);

HK(1:2*NPOIN,N1)=0; HK(N1,1:2*NPOIN)=0; HK(N1,N1)=1; % 将零位移约束对应的行、列变成零，主元变成1 ALLFORCE(N1)=0; end

DISP=HK\\ALLFORCE; % 方程求解，HK先求逆再与力向量左乘得到位移 %-------------------------求应力--------------------------------- stress=zeros(3,NELEM); % 应力向量并清零

for i=1:NELEM % 对单元个数进行循环

D=(YOUNG/(1-POISS*POISS))*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]; %弹性矩阵 for k=1:8 % 对单元结点个数进行循环 N2=LNODS(i,k); % 取单元结点号

U(k*2-1:k*2)=DISP(N2*2-1:N2*2); %从总位移向量中取出当前单元的结点位移 B=eBe(X,Y,0,0); %调用单元中心的应变矩阵 end

stress(:,i)=D*B*U'; end

%===============计算单元刚度矩阵函数=================== function EK=eKe(X,Y,YOUNG,POISS,THICK) EK=zeros(16,16); % 张成16*16矩阵并清零

D=(YOUNG/(1-POISS*POISS))*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]; %弹性矩阵 %高斯积分采用3*3个积分点

A1=5/9;A2=8/9;A3=5/9; %对应积分的加权系数 A=[A1 A2 A3]; r=(3/5)^(1/2);

x=[-r 0 r]; %积分点 for i=1:3 for j=1:3

B=eBe(X,Y,x(i),x(j)); %调用应变矩阵B J=Jacobi(X,Y,x(i),x(j)); %调用雅可比矩阵J EK=EK+A(i)*A(j)*B'*D*B*det(J)*THICK;

3

end end

%===============计算雅可比矩阵函数=================== function J=Jacobi(X,Y,s,t) [N_s,N_t]=DHS(s,t); x_s=0;y_s=0; x_t=0;y_t=0; for j=1:8

x_s=x_s+N_s(j)*X(j);y_s=y_s+N_s(j)*Y(j); x_t=x_t+N_t(j)*X(j);y_t=y_t+N_t(j)*Y(j); end

J=[x_s y_s;x_t y_t];

%===============计算应变矩阵函数=================== function B=eBe(X,Y,s,t) [N_s,N_t]=DHS(s,t); J=Jacobi(X,Y,s,t); B=zeros(3,16); for i=1:8

B1=J(2,2)*N_s(i)-J(1,2)*N_t(i); B2=-J(2,1)*N_s(i)+J(1,1)*N_t(i); B(1:3,2*i-1:2*i)=[B1 0;0 B2;B2 B1]; end

B=B/det(J);

%===============计算形函数的求导函数================ == function [N_s,N_t]=DHS(s,t)

N_s(1)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t); N_s(2)=1/2*(1+t)*(1-t);

N_s(3)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t); N_s(4)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t); N_s(5)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t); N_s(6)=-1/2*(1+t)*(1-t);

N_s(7)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t); N_s(8)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t); N_t(1)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t); N_t(2)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t); N_t(3)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t); N_t(4)=1/2*(1+s)*(1-s);

N_t(5)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t); N_t(6)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t); N_t(7)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t); N_t(8)=-1/2*(1+s)*(1-s); end

%===============计算总荷载矩阵函数===================

function ALLFORCE=FOECEXL(NPOIN,NFORCE,FORCE) % 本函数生成荷载向量 ALLFORCE=zeros(2*NPOIN,1); % 张成特定大小的向量，并赋值0 for i=1:NFORCE

4

ALLFORCE((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);

%FORCE(i,1)为作用点,FORCE(i,2:3)为x,y方向的结点力 end

%-------------------------输出节点位移和单元中心应力------------------------------- for i=1:NPOIN

fprintf(FP2,'x%d=%d\\n',i, DISP(2*i-1)); %输出结点x方向的位移 fprintf(FP2,'y%d=%d\\n',i, DISP(2*i)); %输出结点y方向的位移 end

for j=1:NELEM

fprintf(FP2,'%d x=%f\\n',j,stress(1,j)); %输出单元x方向的应力 fprintf(FP2,'%d y=%f\\n',j,stress(2,j)); %输出单元y方向的应力 fprintf(FP2,'%d xy=%f\\n',j,stress(3,j)); %输出单元切应力 end

%------------------------实例计算并用ANSYS进行对比结果----------------------------

如图所示一个4m*1m悬臂梁，在3节点作用1*105N的竖向力，参数如下：弹性模量2.0*108，泊松比0.3，划分成四个单元，每个单元八个节点，单元尺寸是1m*1m。

1.用MATLAB进行计算

在MATLAB当前工作目录下存入初始数据文件，nonde.txt文本数据内容在表第1列，计算结果（位移和应力）在表第2、3、列，第4列为ANSYS计算的结点竖向位移 23 4 1 6 2E08 0.3 1 x1=-2.335725e-02 x16=5.504682e-07 ANSYS 竖向位移 1 2 3 4 5 6 7 8 y1=-1.298092e-01 y16=-1.135171e-02 1 -0.12951 7 6 5 9 10 11 12 13 x2=-8.552411e-05 x17=-1.009374e-02 2 -0.13063 12 11 10 14 15 16 17 18 y2=-1.304313e-01 y17=-1.224831e-02 3 -0.13216 17 16 15 19 20 21 22 23 x3=2.414924e-02 x18=-1.428159e-02 4 -0.10634 4 0 y3=-1.314740e-01 y18=-2.498899e-02 5 -0.83270E-01 4 0.5 x4=2.339276e-02 x19=5.239181e-03 6 -0.82963E-01 4 1 y4=-1.063855e-01 y19=-3.600960e-03 7 -0.83027E-01 3.5 1 x5=2.212796e-02 x20=0 8 -0.10684 3 1 y5=-8.301568e-02 y20=0 9 -0.61320E-01 3 0.5 x6=1.768304e-05 x21=0 10 -0.41589E-01 3 0 y6=-8.281498e-02 y21=0 11 -0.41216E-01 3.5 0 x7=-2.217334e-02 x22=0 12 -0.41644E-01 5

2.5 1 2 1 2 0.5 2 0 2.5 0 1.5 1 1 1 1 0.5 1 0 1.5 0 0.5 1 0 1 0 0.5 0 0 0.5 0 3 0 -1E05 39 40 41 42 43 44 y7=-8.298558e-02 x8=-2.314206e-02 y8=-1.064554e-01 x9=2.026594e-02 y9=-6.115975e-02 x10=1.770574e-02 y10=-4.153420e-02 x11=-3.209328e-06 y11=-4.112391e-02 x12=-1.769400e-02 y12=-4.152824e-02 x13=-2.028457e-02 y13=-6.114672e-02 x14=1.428595e-02 y14=-2.498661e-02 x15=1.009174e-02 y15=-1.224757e-02 y22=0 13 -0.61222E-01 x23=-5.240129e-03 14 -0.25027E-01 y23=-3.600612e-03 15 -0.12282E-01 ----------应力--------- 16 -0.11371E-01 17 -0.12285E-01 1x=-18069.777552 18 -0.25067E-01 1y=8572.177205 19 -0.36203E-02 1xy=-83189.411527 20 0.0000 2x=3732.743571 21 0.0000 2y=-1485.769440 22 0.0000 2xy=-87734.904796 23 -0.36116E-02 3x=-669.435180 3y=275.080059 3xy=-92666.357985 4x=98.015042 4y=-40.262019 4xy=-67107.620819 2.用ANSYS软件分析对比 四边形八结点采用的是solid 8node 183单元，输入参数，在ANSYS中输出结果

对比结果：在ANSYS中编号3（右上方）的Y方向的最大位移为-0.132165与 MATLAB中编号3的Y方向的最大位移为-1.314740e-01基本吻合。

6

2.5 1 2 1 2 0.5 2 0 2.5 0 1.5 1 1 1 1 0.5 1 0 1.5 0 0.5 1 0 1 0 0.5 0 0 0.5 0 3 0 -1E05 39 40 41 42 43 44 y7=-8.298558e-02 x8=-2.314206e-02 y8=-1.064554e-01 x9=2.026594e-02 y9=-6.115975e-02 x10=1.770574e-02 y10=-4.153420e-02 x11=-3.209328e-06 y11=-4.112391e-02 x12=-1.769400e-02 y12=-4.152824e-02 x13=-2.028457e-02 y13=-6.114672e-02 x14=1.428595e-02 y14=-2.498661e-02 x15=1.009174e-02 y15=-1.224757e-02 y22=0 13 -0.61222E-01 x23=-5.240129e-03 14 -0.25027E-01 y23=-3.600612e-03 15 -0.12282E-01 ----------应力--------- 16 -0.11371E-01 17 -0.12285E-01 1x=-18069.777552 18 -0.25067E-01 1y=8572.177205 19 -0.36203E-02 1xy=-83189.411527 20 0.0000 2x=3732.743571 21 0.0000 2y=-1485.769440 22 0.0000 2xy=-87734.904796 23 -0.36116E-02 3x=-669.435180 3y=275.080059 3xy=-92666.357985 4x=98.015042 4y=-40.262019 4xy=-67107.620819 2.用ANSYS软件分析对比 四边形八结点采用的是solid 8node 183单元，输入参数，在ANSYS中输出结果

对比结果：在ANSYS中编号3（右上方）的Y方向的最大位移为-0.132165与 MATLAB中编号3的Y方向的最大位移为-1.314740e-01基本吻合。

6

本文来源：https://www.bwwdw.com/article/efka.html