《平面四边形八节点等参元程序》由会员分享，可在线阅读，更多相关《平面四边形八节点等参元程序（16页珍藏版）》请在金锄头文库上搜索。
1、 .平面四边形八节点 等参元程序=主程序=% 变量说明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
2、,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) ; %
3、结点力数组 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子块叠加到总刚中（对号入座
4、）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 endend% 由结点力生成总荷载向量列阵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
5、: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 endend% 计算结点位移DISP=HKFORCE;% 计算单元应力EDISP=zeros(16,1); %单元位移列向量清0STRES=zeros(3,NELEM);for i=1:NELEM %对单元个数进行循环 for j=1:8 %对结点进行循环 N1=LNODS(i,j); EDISP(j*2-1):(j*
6、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,%fn,DISP);fclose(FP1);fclose(FP2); % 关闭文件=ele_EK.m=% 计算单元刚度矩阵functio
7、n 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
8、(2,2)-J(1,2)*J(2,1); %求雅克比矩阵行列式的值 EK=EK+W(i)*W(j)*B*D*B*DJ*h; %高斯计分法求单元刚度矩阵 endendreturn=BM.m=% 求应变矩阵Bfunction 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,N
9、XY(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=% 求弹性矩阵Dfunction D=DM(E,mu)A1=mu; %平面应力的弹性系数A2=(1-mu)/2;A3=E/(1-mu*mu);%平面应变的弹性系
10、数%A1=POISS/(1-POISS);%A2=(1-2*POISS)/2/(1-POISS);%A3=E*(1-POISS)/(1+POISS)/(1-2*POISS);%弹性矩阵D=A3*1 A1 0; A1 1 0; 0 0 A2; return=NXY.m=% 形函数对x、y求导function NXY=nxy(ie,LNODS,COORD,KC,ETA)J=Jacob(ie,LNODS,COORD,KC,ETA); %调用函数求雅可比矩阵JJN=J(2,2),-J(2,1);. %求J矩阵的伴随矩阵 -J(1,2),J(1,1);DJ=J(1,1)*J(2,2)-J(1,2)*J(2,1); %求J矩阵行列式的值%借助复合函数求导法则求NXYNKC=nkc(KC,ETA); %调用函数求N对KC导数NETA=neta(KC,EYA); %调用函数求N对ETA导数NXY=1/DJ*JN*NKC;NETA;return=xhs.m=% 平面四节点等参元形函数function N=xhs(KC,ETA)N(1)=(1-KC)*(1-ETA)*(-KC-ETA-1)/4;N(2)=(1+KC)*(1-ETA)*(KC-ETA-1)/4;N(3)=(1+KC)*(1+ETA)*(KC+ETA-1)/4;
