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

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

====================主程序====================

% 变量说明2016

% NPOIN NELEM NVFIX % 总结点数 单元数 受约束结点数

% NFPOIN YOUNG POISS THICK % 结点力数 弹性模量 泊松比 厚度 % LNODS COORD FPOIN FIXED % 单元定义数组 结点坐标数组 结点力数组 约束信息数组 % HK FORCE DISP % 总刚度矩阵 总荷载向量 结点位移向量

clc;clear all; %清空变量 format short e % 设定输出类型

FP1=fopen('sj.txt','rt'); % 打开数据文件 FP1数据文件指针 FP2=fopen('jg.txt','wt'); % 写入结果的文件通道号

% 读入初始数据

NPOIN=fscanf(FP1,'%d',1); NELEM=fscanf(FP1,'%d',1); NFPOIN=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])' ; % 结点坐标数组 FPOIN=fscanf(FP1,'%f',[3,NFPOIN])' ; % 结点力数组 FIXED=fscanf(FP1,'%d',[3,NVFIX])' ; % 约束数组

HK=zeros(2*NPOIN,2*NPOIN); FORCE=zeros(2*NPOIN,1);

% 形成总刚度矩阵

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

EK=ele_EK(i,LNODS,COORD,YOUNG,POISS,THICK);

%调用函数生成单刚

%按2*2子块加入总刚中（共计64块）

for j=1:8 %对行进行循环---按结点号循环 N1=LNODS(i,j);

for k=1:8 %对列进行循环---按结点号循环 N2=LNODS(i,k);

% 单刚2 x 2子块叠加到总刚中（对号入座）

HK((N1*2-1):(N1*2),(N2*2-1):(N2*2))=HK((N1*2-1):(N1*2),(N2*2-1):(N2*2))+EK((j*2-1):(j*2),(k*2-1):(k*2)); end end end

% 由结点力生成总荷载向量列阵

for i=1:NFPOIN %对结点荷载个数进行循环 N1=FPOIN(i,1); %作用荷载的结点号

FORCE((2*N1-1):2*N1,1)=FORCE((2*N1-1):2*N1)+FPOIN(i,2:3)'; end

% 总刚、总荷载进行边界条件处理

for i=1:NVFIX %对约束个数进行循环 N1=FIXED(i,1);

for j=-1:0 if FIXED(i,j+3)==1

HK((N1*2+j),1:2*NPOIN)=0; HK(1:2*NPOIN,(N1*2+j))=0; HK((N1*2+j),(N1*2+j))=1; FORCE((N1*2+j),1)=0;

%将零位移约束对应行、列变成零，主元素为1。

end end end

% 计算结点位移 DISP=HK\\FORCE;

% 计算单元应力

EDISP=zeros(16,1); %单元位移列向量清0 STRES=zeros(3,NELEM);

for i=1:NELEM %对单元个数进行循环 for j=1:8 %对结点进行循环 N1=LNODS(i,j);

EDISP((j*2-1):(j*2),1)=DISP((2*N1-1):(2*N1),1); % end

for j=1:3 %对高斯积分点进行循环 for k=1:3

X=[-sqrt(0.6),0,sqrt(0.6)];

B=BM(i,LNODS,COORD,X(j),X(k));

%调用函数求应变矩阵B end end

D=DM(YOUNG,POISS); %调用函数求弹性矩阵D STRES(:,i)=D*B*EDISP; %根据S=D*B*EDISP求应力矩阵 end

fprintf(FP2,'%f\\n',DISP); fclose(FP1);

fclose(FP2); % 关闭文件

===================ele_EK.m===================

% 计算单元刚度矩阵

function EK=ele_EK(ie,LNODS,COORD,E,mu,h)

EK=zeros(16,16); %8节点的平面单元有8*2个自由度； % 3*3 高斯积分点和权系数

X=[-sqrt(0.6) 0 sqrt(0.6)]; %高斯积分点 W=[5/9,8/9,5/9]; %加系数权

D=DM(E,mu); %调用函数求D弹性矩阵 %生成单元刚度矩阵

for i=1:3 %对高斯积分点进行循环 for j=1:3

B=BM(ie,LNODS,COORD,X(i),X(j));

%调用函数求B应变矩阵 J=Jacob(ie,LNODS,COORD,X(i),X(j));

%调用函数求雅可比矩阵 DJ=J(1,1)*J(2,2)-J(1,2)*J(2,1);

%求雅克比矩阵行列式的值 EK=EK+W(i)*W(j)*B'*D*B*DJ*h;

%高斯计分法求单元刚度矩阵 end end return

===================BM.m===================

% 求应变矩阵B

function B=BM(ie,LNODS,COORD,KC,ETA) NXY=nxy(ie,LNODS,COORD,KC,ETA);

%形函数N对x、y分别求导。

B=[NXY(1,1),0,NXY(1,2),0,NXY(1,3),0,NXY(1,4),0,NXY(1,5),0,NXY(1,6),0,NXY(1,7),0,NXY(1,8),0;...

0,NXY(2,1),0,NXY(2,2),0,NXY(2,3),0,NXY(2,4),0,NXY(2,5),0,NXY(2,6),0,NXY(2,7),0,NXY(2,8);...

NXY(2,1),NXY(1,1),NXY(2,2),NXY(1,2),NXY(2,3),NXY(1,3),NXY(2,4),NXY(1,4),NXY(2,5),NXY(1,5),NXY(2,6),NXY(1,6),NXY(2,7),NXY(1,7),NXY(2,8),NXY(1,8)]; return

===================DM.m===================

% 求弹性矩阵D function D=DM(E,mu)

A1=mu; %平面应力的弹性系数 A2=(1-mu)/2; A3=E/(1-mu*mu);