资源描述：
《黄科4节点程序.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、四边形四节点等参元matlab计算程序姓名：黄科学号：%***变量说明***%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方
2、向的节点力；(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
3、,'%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
4、,'%f',[2,NPOIN])';%结点号x坐标，y坐标(整体坐标下)；共[NPOIN行,2列]FORCE=fscanf(FP1,'%f',[3,NFPOIN])';%节点力：结点号、X方向力(向右为正)，Y方向力(向上为正)FIXED=fscanf(FP1,'%d',[3,NVFIX])';%约束信息数组：受约束节点号，x，y方向受约束情况（受约束为1）；共[NVFIX%行,3列]%—————————————平面应力问题求解—————————————%形成刚度矩阵%计算刚度矩阵，并对约束条件进行处理
5、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,(
6、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
7、)=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);%
8、集成总体的等效节点荷载向量endendASLOAD=ASLOAD+GTJLOAD;%求出结构的节点荷载向量%———————————总刚、总荷载进行边界条件处理————————————%将约束信息加入总刚,总荷载（边界条件处理）fori=1:NVFIXifFIXED(i,2)==1N1=FIXED(i,1)*2-1;HK(N1,;)=0;%将一约束号处的总刚列向量清0HK(:,N1)=0;%将一约束号处的总刚行向量清0HK(N1,N1)=1;
