黄科4节点程序

1四边形四节点等参元matlab计算程序姓名：黄科学号：2111316108%***变量说明***%NNODE：单元节点数NFORCE：受力结点数%LNODS：单元定义数组COORD：结构节点整体坐标数组%ASLOAD：总体荷载向量DISP：结点位移向量%HK：总体刚度矩阵Ke：单元刚度矩阵%TJLOAD：体积力%FIXED：约束信息数组，(n,1)：约束点号；(n,2)与(n,3)：分别为约束点x%方向和y方向的约束情况,受约束为1否则为0；n：受约束节点数目；%FORCE：结点力数组；(i,1)为作用点；(i,2)为x方向的节点力；(i,3)为y方向%的节点力；n：节点力数目；%E：弹性模量NPOIN：总结点数%h：厚度NELEM：单元数%v：泊松比NVFIX：约束结点个数%—————————————读取初始数据——————————————clearclcformatshorteFP1=fopen('4jd.txt','rt');%打开数据文件%***读入控制数据***NPOIN=fscanf(FP1,'%d',1);%总结点数NELEM=fscanf(FP1,'%d',1);%单元数NFORCE=fscanf(FP1,'%d',1);%结点力数NVFIX=fscanf(FP1,'%d',1);%受约束结点数NNODE=fscanf(FP1,'%d',1);%单元结点数E=fscanf(FP1,'%e',1);%弹性模量v=fscanf(FP1,'%f',1);%泊松比h=fscanf(FP1,'%f',1);%厚度%***读取结构信息***LNODS=fscanf(FP1,'%f',[NNODE,NELEM])';%单元定义：单元结点编码（逆时针）；共[NELEM行,NFPOIN列]COORD=fscanf(FP1,'%f',[2,NPOIN])';%结点号x坐标，y坐标(整体坐标下)；共[NPOIN行,2列]2FORCE=fscanf(FP1,'%f',[3,NFPOIN])';%节点力：结点号、X方向力(向右为正)，Y方向力(向上为正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组：受约束节点号，x，y方向受约束情况（受约束为1）；共[NVFIX%行,3列]%—————————————平面应力问题求解—————————————%形成刚度矩阵%计算刚度矩阵，并对约束条件进行处理HK=zeros(2*NPOIN,2*NPOIN);%张成总刚度矩阵并清零%调用子程序生成单元刚度矩阵forie=1:NELEM%ie为单元号Ke=StiffnessMatrix(ie,E,v,COORD,LNODS,h);%调用单元刚度矩阵a=LNODS(ie,:);%临时向量，用来记录当前单元的节点编号forj=1:4%对行按结点号进行循环fork=1:4%对得按结点号进行循环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)+…Ke(j*2-1:j*2,k*2-1:k*2);%单刚集合成总刚endendend%—————————————对荷载向量进行处理—————————————ASLOAD=zeros(2*NPOIN,1);%张成总荷载向量并清0fori=1:NFORCE%对每个节点力进行循环B1=FORCE(i,1)*2-1;B2=FORCE(i,1)*2;%FORCE（i,1）为作用点ASLOAD(B1)=FORCE(i,2);%FORCE（i,2）为x方向的节点力ASLOAD(B2)=FORCE(i,3);%FORCE（i,3）为y方向的节点力endTJLOAD=ele_LOAD(NELEM,COORD,LNODS,R,h);GTJLOAD(2*NPOIN,1)=0;%张成整体等效节点荷载向量并清0forie=1:NELEM%对单元的个数进行循环forj=1:NNODE%对单元的结点个数进行循环M1=LNODS(ie,j)*2;%对节点的第二个位移进行编号GTJLOAD((M1-1):M1,1)=GTJLOAD((M1-1):M1,1)+TJLOAD(j*2-1:j*2,1);%集成总体的等效节点荷载向量endendASLOAD=ASLOAD+GTJLOAD;%求出结构的节点荷载向量%———————————总刚、总荷载进行边界条件处理————————————%将约束信息加入总刚,总荷载（边界条件处理）fori=1:NVFIXifFIXED(i,2)==1N1=FIXED(i,1)*2-1;HK(N1,;)=0;%将一约束号处的总刚列向量清03HK(:,N1)=0;%将一约束号处的总刚行向量清0HK(N1,N1)=1;%将行列交叉处的元素置为1ASLOAD(N1)=0;endifFIXED(i,3)==1N2=FIXED(i,1)*2;HK(N2,;)=0;%将一约束号处的总刚列向量清0HK(:,N2)=0;%将一约束号处的总刚行向量清0HK(N2,N2)=1;%将行列交叉处的元素置为1ASLOAD(N2)=0;endend%—————————————————————————————————DISP=HK\ASLOAD%计算节点位移向量%—————————————求解单元应力——————————————stress=zeros(3,NELEM);%张成应力矩阵赋空值strain=zeros(3,NELEM);%张成应变矩阵赋空值forie=1:NELEM%对每个单元进行循环EDISP(8,1)=0;%张成单元的位移矩阵并清0d=LNODS(ie,:);forj=1:NNODE%对每个单元节点进行循环EDISP(j*2-1:j*2)=DISP(d(j)*2-1:d(j)*2);%从总位移向量中取出当前单元的节点位移endD=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵DB=MatrixB(ie,1,1,COORD,LNODS);%调用计算应变矩阵B函数,赋值（1，1）处S=D*B;%计算应力矩阵Sstress(:,ie)=S*EDISP;%将每个单元的应力矩阵集合在一起strain(:,ie)=B*EDISP;%将每个单元的应变矩阵集合在一起endstress%输出应力矩阵strain%输出应变矩阵%———————————————子程序———————————————functionKe=StiffnessMatrix(ie,E,v,COORD,LNODS,h)%计算平面应变等参数单元的刚度矩阵%说明：用高斯积分法求解平面等参数单元的刚度矩阵Ke=zeros(8,8);%张成单元刚度矩阵并清零D=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵D%3x3高斯积分点和权系数xi=[-0.774596669241483,0,0.774596669241483];%高斯点w=[0.555555555555555,0.8888888888888889,0.555555555555555];%权系数fori=1:length(xi)forj=1:length(xi)B=MatrixB(ie,xi(i),xi(j),COORD,LNODS);4J=Jacobi(ie,xi(i),xi(j),COORD,LNODS);Ke=Ke+w(i)*w(j)*transpose(B)*D*B*det(J)*h;endendreturn%—————————————————————————————————functionTJLOAD=ele_LOAD(ie,COORD,LNODS,R,h)%计算单元等效节点体积力%输入参数：ie——单元号;R——密度;h——厚度%返回值：TJLOAD——结构的等效节点荷载向量TJLOAD(8,1)=0;x=COORD(LNODS(ie,1:4),1);%取出每个单元结点在整体坐标系下的x方向的坐标y=COORD(LNODS(ie,1:4),2);%取出每个单元结点在整体坐标系下的y方向的坐标%利用三次高斯积分集成单元的体积力的等效节点力矩阵%3x3高斯积分点和权系数xi=[-0.774596669241483,0,0.774596669241483];%取高斯点w=[0.555555555555555,0.8888888888888889,0.555555555555555];%权系数fori=1:length(xi)forj=1:length(xi)[N_kc,N_eta]=N_kceta(ie,xi(i),xi(j));J=Jacobi(ie,xi(i),xi(j),COORD,LNODS);N=Shape(xi(i),xi(j));TJLOAD(1:2:7)=TJLOAD(1:2:7)+N’*N*x*det(J)*w(i)*w(j)*h;TJLOAD(2:2:8)=TJLOAD(2:2:8)+N’*N*y*det(J)*w(i)*w(j)*h;%单元的体积力等效成节点的等效节点力endendTJLOAD=TJLOAD*R;return%—————————————————————————————————functionN=Shape(kc,eta)%计算形函数的值%输入参数：ie——单元号;kc,eta——单元内局部坐标%返回值：N——形函数的值N1=(1-kc)*(1-eta)/4;N2=(1+kc)*(1-eta)/4;N3=(1+kc)*(1+eta)/4;N4=(1-kc)*(1+eta)/4;N=[N1N2N3N4];return%—————————————————————————————————functionB=MatrixB(ie,kc,eta,COORD,LNODS)%计算单元的应变矩阵B%输入参数：ie——单元号;kc,eta——局部坐标5%返回值：B——在局部坐标处的应变矩阵B[N_kc,N_eta]=N_kceta(ie,kc,eta);%调用形函数对局部坐标求导的函数J=Jacobi(ie,kc,eta,COORD,LNODS);%调用雅克比矩阵函数J1=[J(2,2)-J(1,2);-J(1,2)J(1,1)]/[J(1,1)*J(2,2)-J(1,2)*J(2,1)];B=zeros(3,8);%张成应变矩阵并清0fori=1:4B1=J1(1,1)*N_kc(i)+J1(1,2)*N_eta(i);B2=J1(2,1)*N_kc(i)+J1(2,2)*N_eta(i);B(1:3,(2*i-1):2*i)=[B10;0B2;B2B1];endreturn%—————————————————————————————————function[N_kc,N_eta]=N_kceta(ie,kc,eta)%计算形函数对局部坐标的导数%输入参数：ie——单元号;kc,eta——局部坐标%返回值：N_kc——在局部坐标处的形函数对kc坐标的导数%N_eta——在局部坐标处的形函数对eta坐标的导数N_kc=zeros(1,4);N_eta=zeros(1,4);%张成形函数对局部坐标导数的矩阵并清0N_kc=zeros(1,4);N_eta=zeros(1,4);N_kc(1)=-(1-eta)/4;N_eta(1)=-(1-kc)/4;N_kc(2)=(1-eta)/4;N_eta(2)=-(1+kc)/4;N_kc(3)=(1+eta)/4;N_eta(3)=(1+kc)/4;N_kc(4)=-(1+eta)/4;N_eta(4)=(1-kc)/4;return%—————————————————————————————————functionJ=Jacobi(ie,kc,eta,COORD,LNODS)%计算雅克比矩阵%输入参数：ie——单元号;kc,eta——局部坐标%返回值：J——在局