平面四边形八节点-等参元程序_burstlkd_平面8节点四边形单元的等参元程序_

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

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

% 变量说明 2016

% NPOIN NELEM NVFIX

% 总结点数 单元数 受约束结点数

% NFPOIN YOUNG POISS THICK

% 结点力数 弹性模量 泊松比 厚度

% LNODS COORD FPOIN FIXED

% 单元定义数组 结点坐标数组 结点力数组 约束信息数组

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);

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 块）

%

单刚

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))+E

K((j*2-1):(j*2),(k*2-1):(k*2));

end

end

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

for i=1:NFPOIN %

对结点荷载个数进行循环

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;

end

% 计算结点位移

DISP=HK\FORCE;

%

计算单元应力

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

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)];

end

fclose(FP2); % 关闭文件

%

计算单元刚度矩阵

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