有限元法课程设计题目.doc

《有限元法课程设计题目.doc》由会员分享，可在线阅读，更多相关《有限元法课程设计题目.doc（32页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、计算力学课程设计任务书1用三角元程序计算求如下各单元的应力。 （图中几何尺寸、力学量自定）2用三角元程序计算求如下各单元的应力。 （图中几何尺寸、力学量自定）3用三角元程序计算求如下各单元的应力。 （图中几何尺寸、力学量自定）4用平面桁架程序计算求如下各单元的应力。（图中几何尺寸、力学量自定）5用平面刚架程序计算求如下各单元的应力。6KN4.8KN/m2.5cm218KN2.5cm36利用有限元软件进行： 带孔平板的有限元建模与受力分析7利用任意有限元软件进行： 受内外压筒体的有限元建模与应力变形分析8、直杆受集中拉伸外载情形下圣维南原理的数值验证及分析9、自选题目 各组可自行选择题目, 但自

2、选题目需在经老师认可。报告格式：封皮（“计算力学课程设计”，班级，完成人员，指导教师）内容：（设计题目，内容；设计用到的有限元法基本内容，包括原理，公式等；程序，结果；如用软件做设计，写出软件操作过程，截图），A4输出，左侧面订。以班级为单位，6月21日前交到建业楼409。三角元程序4.3.1本节利用平面三角元，对结构体进行离散处理，通过程序运算，可得出单元应力与节点处位移。程序适用于平面应力问题，如果对形成弹性矩阵子程序MODPS作适当修改，就可用于计算平面应变问题。程序适用于各种边界条件及各种载荷。4.3.2平面三角元程序框图 单 元 循 环 单 元 循 环 开 始对数组进行定义输入已知量

3、：NPOINT，NNODE，NELEM，NLOAD，NZERO，ZX，ZY，JJS，LNODE，COORD，LZERO计 算 半 带 宽NHBW对右端项ALOAD，总刚ASTIF充零调用形成弹性矩阵D的子程序MODPS调用形成应变矩阵B的子程序BMATPS 调用形成应力矩阵S的子程序DBE 形成单元刚度矩阵ESTIF 形成半带总刚度矩阵ASTIF形成右端项ALOAD调用方程GUASS子程序，输出节点位移ALOAD调用形成单元矩阵D的子程序MODPS调用形成应变矩阵B的子程序BMATPS调用形成应力矩阵S的子程序DBE将整体位移ALOAD转换成单元节点位移DISP计算单元应力STRESS输出单元

4、应力STRESS 结 束 4.3.3子程序说明本程序包含四个子程序，其主要功能分别是：1.MODPS（DMATX，YOUNG，POISS） 该子程序是形成弹性矩阵D的，D阵的平面应力问题表达式为： (4-56) 如果要计算平面应变问题，可将上式中弹性模量、泊松比值分别用、值代替。2、BMATPS（INELE，COORD，LNODE，BMATX，AREA，NPOIN，NELEM）该子程序用来计算应变矩阵B，B的表达式为： (4-57)3、DBE（DMATX，BMATX，DMATX）该子程序用来形成应力矩阵S，S阵的表达式为 (4-58)4、GUASS（NZERO，LZERO，ASTIF，NHBW

5、，ALOAD，NP2，NPOIT）该子程序是用来进行支座处理和解方程的，最后输出节点处位移分量。支座处理的基本原理是如果位移受到约束，则将其对应刚度对角元素变为1，该对角元素所对应的行和列的其余元素变为0，并且将右端项约束位移对应元素充为0。解方程应用带消去法，解出位移后，将节点位移分量存放在原来存放右端项的标识符ALOAD中，最后将ALOAD（节点位移）输出。4.3.4数组及变量说明1、 数组说明COORD 节点坐标数组LNODE 单元节点数组ALOAD 荷载（右端项）数组，解出位移后，将节点处的位移充于其中ESTIF 单元刚度矩阵，为66阶方阵ASTIF 总刚度矩阵，带状存贮DBMAT 应

6、力矩阵S阵，为36阶DMATX 弹性矩阵D阵，为33阶BMATX 应变矩阵S阵，为36阶STRES 单元应力矩阵，为31阶JJS 受载节点矩阵ZX 方向已知载荷向量ZY Y方向已知载荷向量DISP 单元节点位移矩阵 LZERO 约束位移矩阵2、 变量说明NPOIT 最大节点数NELEM 最大单元数NNODE 单元节点数（三角元中为3）NLOAD 最大受载节点数NZERO 最大约束位移个数YOUNG 弹性模量POISS 泊松比TE 板厚（假设t=1）NHBW 半带宽NP2 位移总数AREA 单元面积y232kg/cm2Mx14564-12模型图已知一对角受压的正方形薄板，厚度为1cm，荷载沿厚度

7、均匀分布，为2Kg/cm2,泊松比，求板内点的应力与位场。题中由于XZ 面及YZ面均为该板的对称面，所以只选取1/4部分作为计算模型如图4-12计算模型的输入数据有以下各量:NPOINNELEMNNODENLOADNZEROYOUNGPOISSTE643160.31.0JJS1ZX0ZY-1.0 ICOORD123456COORD(I,1)0.00.00.50.00.01.0COORD(I,2)1.00.50.50.00.00.0 编号数组123456LZERO13781012LNODEILNODE（I，1）1426LNODE（I，2）2553LNODE（I，3）32354.3.6平面三角元有

8、限元法程序c FOR FINITE ELEMENT DIMENSION COORD(6,2),LNODE(4,3),ALOAD(12),ESTIF(6,6), 1ASTIF(12,8),DBMAT(3,6),STRES(3),BMATX(3,6),DMATX(3,3), 1jjS(1),ZX(1),ZY(1),DISP(6),LZERO(6)c 输入已知数据 OPEN(1,FILE=D:NMXDATA.DAT) DATA NPOIN,NELEM,NNODE,NLOAD,NZERO/6,4,3,1,6/ DATA YOUNG,POISS,TE/.0,0.3,1.0/ DATA ZX/0.0/ D

9、ATA ZY/-1.0/ DATA JJS/1/ DATA LNODE/1,4,2,6,2,5,5,3,3,2,3,5/ DATA COORD/0.0,0.0,0.5,0.0,0.5,1.0,1.0,0.5,0.5,0.0,0.0,0.0/ DATA LZERO/1,3,7,8,10,12/C 计算半带宽 NHBW=0 DO 11 INELE=1,NELEM DO 11 I=1,NNODE DO 11 J=1,NNODE LN=IABS(LNODE(INELE,I)-LNODE(INELE,J) IF (LN.GT.NHBW) NHBW=LN11 CONTINUE NHBW=(NHBW+1)*

10、2 WRITE(1,*)半带宽 WRITE(1,12) NHBW12 FORMAT(1x,NHBW=,I3) NP2=2*NPOIN DO 50 I=1,NP2 ALOAD(I)=0.0 DO 50 J=1,NHBW ASTIF(I,J)=0.050 CONTINUEC 对单元循环 DO 70 INELE=1,NELEM CALL MODPS(DMATX,YOUNG,POISS) CALL BMATPS(INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM) CALL DBE(DMATX,BMATX,DBMAT) DO 30 I=1,6 DO 30 j=1,6 ES

11、TIF(I,J)=0.0 DO 30 K=1,3 ESTIF(I,J)=ESTIF(I,J)+DBMAT(K,I)*BMATX(K,J)*AREA*TE30 CONTINUE DO 40 ID=1,NNODE DO 40 II=1,2 IH=2*(ID-1)+II IDH=2*(LNODE(INELE,ID)-1)+II DO 35 JD=1,NNODE DO 35 JJ=1,2 IL=2*(JD-1)+JJ IDL=2*(LNODE(INELE,JD)-1)+JJ-IDH+1 IF (IDL.LE.0) GOTO 35 ASTIF(IDH,IDL)=ASTIF(IDH,IDL)+ESTIF(

12、IH,IL)35 CONTINUE40 CONTINUE70 CONTINUEC 求右端项 DO 90 I=1,NLOAD IL=JJS(I)*2 ALOAD(IL-1)=ALOAD(IL-1)+ZX(I) ALOAD(IL)=ALOAD(IL)+ZY(I)90 CONTINUEC 支座处理、解方程 CALL GAUSS(NZERO,LZERO,ASTIF,NHBW,ALOAD,NP2,NPOIN) WRITE(1,*)单元应力 WRITE(1,*)单元号 x y xy DO 400 INELE=1,NELEM CALL MODPS(DMATX,YOUNG,POISS) CALL BMATPS

13、(INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM) CALL DBE(DMATX,BMATX,DBMAT) DO 410 I=1,NNODE DO 410 J=1,2 LH=2*(I-1)+J MH=2*(LNODE(INELE,I)-1)+J DISP(LH)=ALOAD(MH)410 CONTINUE DO 420 I=1,NNODE STRES(I)=0.0 DO 420 J=1,6 STRES(I)=STRES(I)+DBMAT(I,J)*DISP(J)420 CONTINUE WRITE(1,430)INELE,(STRES(I1),I1=1,NNO

14、DE)430 FORMAT(1X,I5,1X,3F13.5,1X)400 CONTINUE STOP ENDC 求弹性矩阵D SUBROUTINE MODPS(DMATX,YOUNG,POISS) DIMENSION DMATX(3,3) DMATX(1,1)=YOUNG/(1.0-POISS*POISS) DMATX(1,2)=YOUNG*POISS/(1.0-POISS*POISS) DMATX(2,1)=DMATX(1,2) DMATX(2,2)=DMATX(1,1) DMATX(3,3)=YOUNG/(2.0*(1.0+POISS) RETURN ENDC 求应变矩阵B SUBROUT

15、INE BMATPS(INELE,COORD,LNODE,BMATX,AREA,NPOIN,NELEM) DIMENSION COORD(NPOIN,2),LNODE(NELEM,3),BMATX(3,6) IE=LNODE(INELE,1) JE=LNODE(INELE,2) ME=LNODE(INELE,3) BI=COORD(JE,2)-COORD(ME,2) BJ=COORD(ME,2)-COORD(IE,2) BM=COORD(IE,2)-COORD(JE,2) CI=COORD(ME,1)-COORD(JE,1) CJ=COORD(IE,1)-COORD(ME,1) CM=COOR

16、D(JE,1)-COORD(IE,1) AREA=(BJ*CM-BM*CJ)/2.0 DO 3 I=1,3 DO 3 J=1,63 BMATX(I,J)=0.0 CH=2.0*AREA BMATX(1,1)=BI/CH BMATX(1,3)=BJ/CH BMATX(1,5)=BM/CH BMATX(2,2)=CI/CH BMATX(2,4)=CJ/CH BMATX(2,6)=CM/CH BMATX(3,1)=BMATX(2,2) BMATX(3,2)=BMATX(1,1) BMATX(3,3)=BMATX(2,4) BMATX(3,4)=BMATX(1,3) BMATX(3,5)=BMATX(

17、2,6) BMATX(3,6)=BMATX(1,5) RETURN ENDC 求应力矩阵DB SUBROUTINE DBE(DMATX,BMATX,DBMAT) DIMENSION DBMAT(3,6),DMATX(3,3),BMATX(3,6) DO 3 I=1,3 DO 3 J=1,6 DBMAT(I,J)=0.0 DO 3 K=1,3 DBMAT(I,J)=DBMAT(I,J)+DMATX(I,K)*BMATX(K,J)3 CONTINUE ENDC 支座处理、解方程 SUBROUTINE GAUSS(NZERO,LZERO,ASTIF,NHBW,ALOAD,NP2,NPOIN) DIM

18、ENSION LZERO(NZERO),ASTIF(NP2,NHBW),ALOAD(NP2) DO 260 I=1,NZERO IZ=LZERO(I) ASTIF(IZ,1)=1.0 DO 210 J=2,NHBW ASTIF(IZ,J)=0.0210 CONTINUE J0=NHBW IF(IZ-NHBW.LE.0) J0=IZ DO 250 J=2,J0 M=IZ-J+1 ASTIF(M,J)=0.0250 CONTINUE ALOAD(IZ)=0.0260 CONTINUE KK=NP2-1 DO 290 K=1,KK IM=NP2 IF(NP2-K-NHBW+1.GT.0) IM=NH

19、BW+K-1 IK=K+1 DO 285 I=IK,IM L=I-K+1 C=ASTIF(K,L)/ASTIF(K,1) IN=NHBW-L+1 DO 280 J=1,IN M=J+I-K ASTIF(I,J)=ASTIF(I,J)-C*ASTIF(K,M)280 CONTINUE ALOAD(I)=ALOAD(I)-C*ALOAD(K)285 CONTINUE290 CONTINUE ALOAD(NP2)=ALOAD(NP2)/ASTIF(NP2,1) DO 315 IB=1,KK I=NP2-IB J0=NHBW IF(NHBW-NP2+I-1.GT.0) J0=NP2-I+1 DO 31

20、0 J=2,J0 IH=J+I-1 ALOAD(I)=ALOAD(I)-ASTIF(I,J)*ALOAD(IH)310 CONTINUE ALOAD(I)=ALOAD(I)/ASTIF(I,1)315 CONTINUE WRITE(1,*)节点位移 WRITE(1,*)节点号 x方向位移 y方向位移 DO 111 I=1,NPOIN II=2*I-1 IL=II+1 WRITE(1,580)I, ALOAD(II),ALOAD(IL)111 CONTINUE580 FORMAT(I5,7X,F10.7,7X,F10.7) RETURN END平面桁架程序2.3.1本程序适用于节点荷载作用下的桁

21、架分析问题，当有节间荷载存在时可按照静力等效原理将其转化为节点荷载。可求解平面桁架在静力荷载作用下的内力和位移。2.3.2桁架有限元程序框图开 始 输 入 数 据 数 组 定 义计算各杆截面面积和半带宽 调用形成单刚矩阵UNIT调用形成半带宽存贮的结构原始刚度矩阵TOTAL有节点荷载否输入节点荷载值，并将其送入相应的荷载列阵P（N）中考虑结构是否自重将杆自重引起的等效荷载叠加到P（N）中支座处理、解方程，并输出U（N）、V（N）调用UNIT，求各单元杆端内力 结 束单元循环没有有否是 2.3.3子程序说明（1）TOTAL（INELE，LNODE，ESTIF，ASTIF） 该子程序用于形成半带宽

22、存贮的结构原始刚度矩阵 （2）UNIT（INELE，YOUNG，COORD，LNODE，AA，ESTIF，LL，CX,CY） 该子程序利用式 （2平面桁架内力计算的标识符 NPOIN 最大节点数 NELEM 最大单元数 NLOAD 节点的荷载总数NZERO 节点的约束位移总数WT 为结构的自重EE 为材料的弹性模量LL 一维数组，用于存放单元杆件的长度AA 一维数组，用于存放单元杆件的面积COORD 节点坐标数组LNODE 单元节点数组BH 二维数组，用于存放单元截面尺寸NRES 二维数组，用于存放约束的位移值JP 二维数组，用于存放节点的荷载值ESTIF 四维数组，用于存放整体坐标系下的单元

23、刚度矩阵ASTIF 二十维数组，用于存放半带宽结构原始刚度矩阵P 用于存放节点的荷载列阵U 用于存放节点x方向的位移值V 用于存放节点y方向的位移值2.3.5例题及数据准备例如图2-6所示桁架，已知杆件材料的弹性模量E=2.1104KN/cm2,杆件截面高度H=10cm，截面的宽度为b=2cm，不计各杆的自重，求在荷载作用下，各杆的轴力。图2-6桁架例图计算模型的输入数据有以下各量：NPOIN NRESNELEMNLOADYOGWT3331210000.0 节点 单元12LNODE(1,I)12LNODE(2,I)23LNODE(3,I)13 编号数组123RES(I,1)0.00.00.0R

24、ES(I,2)124节点 I123COORD(I,1)0.06.00.0COORD(I，2)0.00.06.0JP(I,1)10.0JP(I,2)5 单元号I123BH(I,1)101010BH(I,2)2222.3.6平面桁架内力计算程序 INTEGER E,NELEM,Z,H REAL LL,ESTIF,ASTIF,JP DIMENSION COORD(3,2),LNODE(3,2),AA(200),BH(3,2),RES(3,2), &LL(200),ESTIF(4,4),ASTIF(400,20),JP(1,2),P(400),U(200), &V(200) OPEN(2,FILE=D

25、:NMXJIA.DAT,STATUS=NEW) C 输入已知数据DATA NPOIN,NELEM,NJP,NRES,EE,WT/3,3,1,3,21000,0/DATA COORD/0,6,0,0,0,6/DATA LNODE/1,2,1,2,3,3/DATA BH/3*2,3*10/DATA RES/3*0,1,2,4/C 计算各单元面积DO 200 E=1,NELEMAA(E)=BH(E,1)*BH(E,2) CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY)200 CONTINUE C计算半带宽 L2=2*NPOIN NHBW=0 DO 210 E

26、=1,NELEM M=ABS(LNODE(E,1)-LNODE(E,2) IF(NHBW.LT.M) NHBW=M210 CONTINUEWRITE(2,*) 半带宽 NHBW=2*(NHBW+1) WRITE(2,220) NHBW220 FORMAT(1X,NHBW=, I2)C单元循环 DO 300 I1=1,L2 DO 300 J1=1,NHBW300 ASTIF(I1,J1)=0.0 DO 400 E=1,NELEM CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) CALL TOTAL(E,LNODE,ESTIF,ASTIF) 400 CO

27、NTINUE DO 560 N=1,L2560 P(N)=0.0 IF(NJP.EQ.0) GOTO 650DATA JP/10,5/ DO 630 K1=1,NJP NN=JP(K1,2)+0.1 630 P(NN)=JP(K1,1) 650 IF(WT.LE.0.0) GOTO 750 DO 700 E=1,NELEM N1=LNODE(E,1) N2=LNODE(E,2) P(2*N1)=P(2*N1)-WT*AA(E)*LL(E)/2.0 P(2*N2)=P(2*N2)-WT*AA(E)*LL(E)/2.0700 CONTINUE WRITE(2,710)710 FORMAT(/4X,

28、荷载总数,8X,水平荷载,8X,铅垂荷载) DO 730 K=1,NO730 WRITE(2,740) K,P(2*K-1),P(2*K)740 FORMAT(4X,I2,8X,F8.3,8X,F8.3) 750 DO 800 I1=1,NRES Z=RES(I1,2)+1E-5 ASTIF(Z,1)=ASTIF(Z,1)*1E8 P(Z)=ASTIF(Z,1)*RES(I1,1)800 CONTINUE DO 850 K1=1,L2-1 IF(L2.GT.(K1+NHBW-1) THEN IM=K1+NHBW-1 ELSE IM=L2 ENDIF DO 850 I1=K1+1,IM L1=I

29、1-K1+1 C1=ASTIF(K1,L1)/ASTIF(K1,1) DO 830 J1=1,NHBW-L1+1 MM=J1+I1-K1 ASTIF(I1,J1)=ASTIF(I1,J1)-C1*ASTIF(K1,MM)830 CONTINUE P(I1)=P(I1)-C1*P(K1)850 CONTINUE P(L2)=P(L2)/ASTIF(L2,1) DO 900 I1=L2-1,1,-1 IF(NHBW.GT.(L2-I1+1) THEN JM=L2-I1+1 ELSE JM=NHBW ENDIF DO 880 J1=2,JM H=J1+I1-1 P(I1)=P(I1)-ASTIF(I

30、1,J1)*P(H)880 CONTINUE P(I1)=P(I1)/ASTIF(I1,1)900 CONTINUE WRITE(2,910)910 FORMAT(/10X,节点位移,10X,水平位移,10X,铅垂位移/) DO 930 N=1,NO U(N)=P(2*N-1) V(N)=P(2*N)930 WRITE(2,950) N,U(N),V(N)950 FORMAT(15X,I2,6X,F12.7,6X,F12.7) WRITE(2,970)970 FORMAT(/4X,单元号,8X,节点号,8X,N(KN),8X,Q(KN)/) DO 980 E=1,NELEM CALL UNIT

31、(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) N1=LNODE(E,1) N2=LNODE(E,2) ULNODE=U(N1)-U(N2) VLNODE=V(N1)-V(N2) D1=ESTIF(1,1)*ULNODE+ESTIF(1,2)*VLNODE D2=ESTIF(1,2)*ULNODE+ESTIF(2,2)*VLNODE FI=CX*D1+CY*D2 FJ=-FI TI=-CY*D1+CX*D2 TJ=-TI WRITE(2,990) E,N1,FI,TI,N2,FJ,TJ990 FORMAT(4X,I2,12X,I2,8X,F8.4,8X,F8.4/18

32、X,I2,8X,F8.4,8X,F8.4)980 CONTINUE WRITE(2,1000)1000FORMAT(/28X,结束,/15X,35(*)/) STOP END 子程序TOTAL形成总刚度矩阵SUBROUTINE TOTAL(E,LNODE,ESTIF,ASTIF) INTEGER E,DH,ZL,DL REAL ESTIF,ASTIF DIMENSION LNODE(3,2),ESTIF(4,4),ASTIF(400,20) DO 40 I1=1,2 DO 40 II=1,2 KH=2*(I1-1)+II DH=2*(LNODE(E,I1)-1)+II DO 40 J1=1,2

33、 DO 40 JJ=1,2 KL=2*(J1-1)+JJ ZL=2*(LNODE(E,J1)-1)+JJ DL=ZL-DH+1 IF(DL.GT.0) ASTIF(DH,DL)=ASTIF(DH,DL)+ESTIF(KH,KL)40 CONTINUE RETURN END 子程序UNIT形成单刚SUBROUTINE UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) INTEGER E REAL LL,ESTIF DIMENSION COORD(3,2),LNODE(3,2),AA(200),LL(200),ESTIF(4,4) N1=LNODE(E,1) N2=

34、LNODE(E,2) CX=COORD(N2,1)-COORD(N1,1) CY=COORD(N2,2)-COORD(N1,2) LL(E)=SQRT(CX*CX+CY*CY) CX=CX/LL(E) CY=CY/LL(E) EAL=EE*AA(E)/LL(E) ESTIF(1,1)=EAL*CX*CX ESTIF(1,2)=EAL*CX*CY ESTIF(2,2)=EAL*CY*CY ESTIF(2,1)=ESTIF(1,2) DO 10 I=1,2 DO 10 J=1,2 ESTIF(I,J+2)=-ESTIF(I,J) ESTIF(I+2,J)=-ESTIF(I,J)10 ESTIF(I

35、+2,J+2)=ESTIF(I,J) RETURN END平面刚架程序3.3.1平面刚架内力的计算程序是由一个主程序和五个子程序组成。3.3.2平面刚架程序框图开 始 输入已知量：NPOIN，NELEM，NJP，NPP，NRES，INF，YOUNG,WT 对数组进行定 义 计算半带宽NHBW 调用形成单刚子程序STIFF调用形成半带宽存贮的结构原始刚度矩阵TOTAL有节点荷载否输入节点荷载值，并将其送入相应的荷载列阵P（N）中单元循环没有有否 有非节点荷载否输入非节点荷载有考虑结构是否自重将杆自重引起的等效荷载叠加到P（N）中支座处理、解方程，并输出节点位移、水平位移调用FIXF，求各单元杆端

36、内力 结 束是 调用形成节点荷载FIXF否 调用形成单元常数子程序CXCY 3.3.3子程序说明1、 本子程序包含5个子程序，其主要功能是：（1） 子程序CXCY（E，L，IE，JE）的功能 用于计算各单元的杆件的长度及在结构坐标系中的倾角的三角函数cos、sin。（2） 子程序ELEA（E，L，IE，JE，A）的功能 用于计算结构坐标系下各单元刚度矩阵中的刚度系数。（3） 子程序TOTAL（E，LNODE，KE，KZ）的功能 用于形成半带宽的结构原始刚度矩阵。（4） 子程序FIFX（NP，FXI，FXJ，FYI，FYJ，MMI，MMJ，WL，IE，JE）的功能 它是用来计算各单元的等效节点荷

37、载。（5） 子程序STIFF（E，KE）的功能用于形成结构坐标系下的单元刚度矩阵 。 3.3.4标示符说明1、数组说明COORD 节点坐标数组LNODE 单元节点数组LL 单元各杆件的长度PP 存放杆上的荷载BH（100，2）存放单元杆件的截面宽度（面积）和截面高度（惯性矩）RES(20，2) 存放约束的位移JP（200，2）存放节点载荷PP（200，4）存放杆件载荷KE（6，6） 存放单元刚度矩阵KZ（200，20）存放半带宽存贮的结构原始刚度矩阵U（60）、V（60）、PHI（60）分别用于存放节点的位移1、 变量说明NPOIN 节点总数NELEM 单元总数NJP 节点荷载总数NPP 杆上荷载总数NRES 最大位移约束个数 YOUNG 杨氏模量WT 自重NP3 位移总数LL 单元杆长INF 特征量 当INF0（整数） 输入杆件截面宽度B和截面高度H,0 当INF0（整数） 输入杆件截面面积A和截面的惯性矩INHBW 最大的半带宽宽度3.3.5例题及数据准备：已知图示刚架，各杆的材料及截面均相同，杨氏模量E2108 kN/m2,A=0.01m2,惯性矩I=0.0032m4，q=30kN/m.