c语言计算平面桁架内力计算程序

桁架内力计算程序(TRUSS)一、程序功能及编制方法桁架内力计算程序（TRUSS），能计算任意平面和空间桁架（包括网架）在结点荷载作用下各结点的位移和各杆的轴力。

程序采用变带宽一维数组存储总刚度矩阵，先处理法引进支座条件。

计算结果输出各结点的位移和各杆的轴力。

二、程序使用方法使用方法与“APF”程序相同。

用文件编辑编辑器建立数据文件后即可运行。

计算结果将写在结果文件中。

三、数据文件填写格式数据文件填写格式大致与APF程序相似。

1．总信息：T，NJ，NE，NR，NB，NP，EO，DS其中：T——桁架类型，平面桁架 T=2，空间桁架 T=3。

NR——支座约束数。

其他变量与APF程序相同。

2．结点坐标数组XYZ（NJ， 3）每个结点填一行，每行三个数分别填写结点的x，y，z三个坐标数值,平面桁架只填x,y 值（单位：m）。

3．单元信息数组G（NE）采用紧缩存储方式，每个单元填一个数。

把单元的左端、右端结点号和杆的类型号三个数紧缩为一个数。

例如某单元左端结点号为15，在端结点号为8，类型号为3，则写成0.15083，一般格式为0.×××××。

4．单元截面信息数组AI(NB）填写各类单元的杆截面面积(m2）。

5．约束信息数组R(NR)采用紧缩存储方式，每个约束（支座链杆）填一个数。

把约束作用的作用点写在该数的整数部分，约束的方向写在小数部分。

x方向的约束为“l”，y方向的约束为“2”。

例如某支座链杆作用在 17号结点上，方向沿整体坐标 y方向，则写为 17.2，一般格式为××.×。

6．结点荷载信息数组F（NP，2）每个结点荷载填一行，每行两个数。

前一个数用紧缩方式填写荷载作用的结点号和作用方向，格式与约束信息的格式相同。

后一个数为荷载的数值。

单位为kN，与整体坐标方向一致者为正值，相反者为负值。

例如，作用在16号结点上，数值为183.5 kN，方向向下的力，则写成：16.2，-183.5（这里，假定坐标轴y轴向上）。

桁架结构内力计算方法
在计算桁架结构内力时，可以采用以下步骤：
1.给定载荷：首先确定桁架结构所受到的外部载荷，包括竖向荷载、
水平荷载和斜向荷载等。

这些载荷可以通过静力学分析或者实际测量得到。

2.确定支座反力：根据结构平衡条件，计算出桁架结构支座的反力。

支座反力是由桁架结构与支座之间的约束关系决定的。

3.确定节点平衡条件：桁架结构中的每个节点都应满足平衡条件，即
节点受力平衡。

根据节点的受力平衡条件，可以得到每个节点处的力平衡
方程。

4.建立杆件的受力方程：根据构件材料的力学性质和几何形状，建立
每根杆件的受力方程。

通常使用杆件受力平衡和伸缩力平衡方程。

5.解方程求解内力：将节点平衡条件和杆件受力方程组合起来，得到
一个线性方程组。

通过求解这个方程组，可以求解出各个构件的内力大小
和方向。

在具体计算过程中，可以采用不同的计算方法来求解桁架结构的内力。

以下是几种常用的计算方法：
1.切线法：切线法是一种基于几何形状的方法，通过假设桁架结构各
个构件处于弧形变形状态，利用切线关系计算出内力。

该方法适用于相对
简单的桁架结构。

2.牛顿-拉夫逊法：牛顿-拉夫逊法是一种基于力的平衡条件的方法，
通过迭代计算桁架结构内力。

该方法适用于复杂的桁架结构。

3.力法：力法是一种基于力平衡方程和几何条件的方法，通过逐个构件计算内力。

该方法适用于任意形状的桁架结构。

以上是桁架结构内力计算的基本方法和一些常用的计算方法。

在实际应用中，还可以根据具体情况选择适合的方法进行计算。

习题课3静定平面桁架的内力计算一、找出桁架的零杆(1)F P000000(2)8根零杆5根零杆F P000(3)12根零杆F PF P 00000000000F P 12根零杆(4)A 000000000006根零杆(5)a aaaS 1S 2F P 2F P F P F N 2F N 1AB C DⅡⅡⅠⅠ00000由于荷载反对称，该桁架除下部水平链杆AB 外，其余杆件受力反对称，故。

0=NCD F 1S F=∑I-I 右：20S F =∑II-II 右：12220222N P P F F F +⋅−⋅=10N F =22220222N P P F F F −+⋅+⋅=22N P F F =F PF P(7)(6)F P0附属部分6根零杆7根零杆F P0000000二、用简捷方法求桁架指定杆轴力150+II-II 左：(1)ABC DE FGHa /2a /2ⅡⅡⅠⅠ12解：CM=∑11 1.51.5()N P N P F a F a F F ⋅⋅==−压220/200.5Dy P y PMF a F a F F =+==−∑I-I 下：250.5 1.118()1N P P F F F =−⋅=−压简单桁架1.5F PF P F P F Pa a /2aa /2 1.5F P125(2)F NF yB =2F P dd dd F PAC DB dⅡⅡⅠⅠ1F yC =F PF yA =F PF xA =2F P 323F P2F P331)020()2)0()3)0()C N P P N P x N P y yC P M F d F d F d F F F F F F F F =+−======↑∑∑∑拉拉I-I 右：联合桁架解：F NF yB =2F P dd ddF PAC DB dⅡⅡⅠⅠ1F yC =F PF yA =F PF xA =2F P 323F P2F P1210()()032()D N P P yN P P P M F F d F dF F F F F ==−⋅=−==−=−∑∑压压II-II 左：整体平衡：10(322)()2ByA P P P P P MF F d F d F d F d F d==+⋅−−=↑∑(3)F P0A -F P-F P-F PF DEC 12F PF PPF 2F N F Pa /2aⅠⅠ00B 0EN MF ==∑1)I-I 右：02=N F 2)结点C ：1102()yy PN P FF F F F ===∑拉3)结点F ：aa /2aa /2联合桁架解：(4)F P 4m4mF P F P F P 4F P6.67F P6.67F P AB FC F N 2134ⅡⅠEⅡD GH 4m4m3m3mⅠ解：1220()0()33DN P xN P MF F FF F ====−∑∑拉压1) I-I 右：2) 结点E ：2222550()346xx PN x P FF F F F F ====∑拉简单桁架F PF N 2E 23PF F P 4m4mF P F P F P 4F P6.67F P6.67F P AB FC F N 2134ⅡⅠEⅡD G H4m4m3m3mⅠ553434343x Py P PN F F F F F F ==⋅==4) 结点D ：F PFN 3F N 5D2F P2F P /3410(48)2()6F N P P P M F F F F −==+=−∑压3) II-II 右：xF=∑(5)F PCA B dⅠ1 1.5F P 1.5F P02F P F Pd 0复杂桁架1)结点C：结构与荷载均对称，两斜杆轴力为零。

1题目结构如图所示： 杆的弹性模量E 为200000Mpa ，横截面面积A 为3250mm 2。

图 1 桁架示意图2实验材料PC 机一台，Microsoft Visual Studio 软件，Ansys 软件。

3实验原理（1）桁架结构特点桁架结构中的桁架指的是桁架梁，是格一种梁式结构。

桁架结构常用于大跨度的厂房、展览馆、体育馆和桥梁等公共建筑中。

由于大多用于建筑的屋盖结构，桁架通常也被称作屋架。

结构上由光滑铰链连接，载荷只作用于节点处，只约束线位移，各杆只有轴向拉伸和压缩。

（2）平面桁架有限元分析1、单元分析局部坐标系中的干单元如图所示：图 2 局部坐标系中的杆单元以下公式描述了整体位移和局部位移之间的关系：U=Tu 其中U=[ U ix U iy U jx U jy ]，T=[cos θ−sin θ00sin θcos θ0000cos θ−sin θ00sin θcos θ],u=[u ix u iy u jx u jy ]U 和u 分别代表整体坐标系和局部坐标系XY 系和局部坐标系xy 下节点i 和节点j 的位移。

T 是变形从局部坐标转换到整体坐标系下的变换阵，类似的局部力和整体力也有以下关系：F=Tf其中F=[ F ixF iy F jx F jy ] ，是整体坐标系下施加在节点i 和j 上的力的分量而且f=[ f ix f iy f jx f jy ]，代表局部坐标系下施加在节点i和j上的分量。

在假设的二力杆条件下，杆只能沿着局部坐标系的x方向变形，内力也总是沿着局部坐标系x的方向，因此将y方向的位移设置为0，局部坐标系下内力和位移通过刚度矩阵有如下关系：[f ixf iyf jxf jy]=|k0−k00000−k0k00000|=[U ixU iyU jxU jy]这里k=k eq=AE/L，写成矩阵形式有：f=Ku将f和u替换成F和U有：T-1F=KT-1U将方程两边乘以T得到：F=TKT-1U其中T-1是变换矩阵T的逆矩阵，替换方程中的TKT-1和U矩阵的值，相乘后得到：[F ixF iy F jx F jy]= k[cos2θsinθcosθ−cos2θ−sinθcosθsinθcosθsin2θ−sinθcosθ−sin2θ−cos2θ−sinθcosθcos2θsinθcosθ−sinθcosθ−sin2θsinθcosθsin2θ][U ixU iyU jxU jy]上述方程代表了施加外力、单元刚度矩阵和任意单元节点的整体位移之间的关系。

实验三平面刚架程序设计一．平面刚架内力和位移计算总框图二．平面桁架静力分析源程序（PFSAP.FOR）C ANALYSIS PROGRAM FOR PLANE FRAMEREAL K(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100), & EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6),& FE(6),D(100),ADE(6),DE(6),RT(6,6),AFE(6),F(3)INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),M(6),& JEAI(100),NOOPEN (6,FILE='PFSAP.IN',STA TUS='OLD')OPEN (8,FILE='PFSAP.OUT',STA TUS='NEW')1 READ (6,*)NOIF(NO.EQ.0)STOPWRITE (8,'(/A5,I3,A1)')'(NO.=',NO,')'CALL READ(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)DO 5 I=1,NP(I)=0.DO 5 J=1,N5 K(I,J)=0DO 10 IE=1,NECALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALL MR(R,IE,JE,X,Y)CALL MAKE(KE,R,AKE)CALL CALM(M,IE,JN,JE)CALL MK(K,AKE,M)10 CONTINUEDO 20 IP=1,NPFCALL MR(R,JPF(1,IP),JE,X,Y)CALL TRAN(R,RT)CALL PE(FE,IP,JPF,PF,AL)CALL MULV6 (RT,FE,AFE)CALL CALM(M,JPF(1,IP),JN,JE)CALL MF(P,AFE,M)20 CONTINUEDO 30 I=1,NPJ30 P(JPJ(I))=P(JPJ(I))+PJ(I)CALL SLOV(K,P,D,N)WRITE(8,'(/2(26(1H*),A))') 'RESULTS OF CALCULATION'WRITE(8,40)40 FORMAT(/5X,'NO.N',4X,'X-DISPLACEMENT',2X,& 'Y-DISPLACEMENT',3X,'ANG.ROT.(RAD)')DO 60 KK=1,NJDO 50 II=1,3F(II)=0.I1=JN(II,KK)50 IF(I1.GT.0)F(II)=D(I1)60 WRITE(8,70)KK,F(1),F(2),F(3)70 FORMAT(I8,2X,3G16.5)WRITE(8,80)80 FORMAT(/'NO.E',5X,'N(1)',8X,'Q(1)',8X,'M(1)',& 8X,'N(2)',8X,'Q(2)',8X,'M(2)')DO 130 IE=1,NECALL MADE(IE,JN,JE,D,ADE)CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALL MR(R,IE,JE,X,Y)CALL MULV6(R,ADE,DE)CALL MULV6(KE,DE,FF)DO 100 IP=1,NPFIF (JPF(1,IP).EQ.IE) THENCALL PE(FE,IP,JPF,PF,AL)DO 90 I=1,690 FF(I)=FF(I)-FE(I)ENDIF100 CONTINUEWRITE(8,110) IE,(FF(I),I=1,6)110 FORMA T(I5,2X,6G12.5)130 CONTINUEGOTO 1ENDSUBROUTINE READ(NJ,N,NE,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI, & JPJ,PJ,JPF,PF)REAL X(100),Y(100),EAI(3,100),PJ(100),PF(2,100)INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100), & TITLE(20)READ(6,'(20A4)') (TITLE(I),I=1,20)WRITE(8,'(/7X,20A4)')TITLEREAD (6,*) NJ,N,NE,NM,NPJ,NPFWRITE(8,'(/3(5X,A4,1H:,I2))')'NJ=',NJ,& 'N=',N,'NE=',NE,'NM=',NM,'NPJ=',NPJ,'NPF=',NPFWRITE (8,5)5 FORMAT(/4X,'NO. (1) (2) (3)',10X,'X ',8X,'Y')READ (6,10)((JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)10 FORMAT(2(3I5,2G16.4))DO 20 I=1,NJ20 WRITE (8,'(2X,1H(,I2,1H),3I6,4X,2F10.3)') I,JN(1,I),JN(2,I),& JN(3,I),X(I),Y(I)WRITE (8,30)30 FORMAT(/10X,'ELEMENT NO. NODE-1 NODE-2 MA TERIALS')READ (6,40)(JE(1,I),JE(2,I),JEAI(I),I=1,NE)40 FORMAT(5(3I5))DO 50 I=1,NE50 WRITE (8,'(14X,I2,3(7X,I3))') I,JE(1,I),JE(2,I),JEAI(I)READ(6,*)((EAI(I,J),I=1,3),J=1,NM)WRITE(8,60)(J,(EAI(I,J),I=1,3),J=1,NM)60 FORMAT(/3X,'NO.MAT',6X,'ELASTIK MODULUS',8X,& 'AREA',5X,'MOMENT OF INERTIA'/(I6,9X,3G16.6))IF(NPJ.EQ.0) GOTO 90WRITE(8,'(/20X,16H NODEL LOADS )')WRITE(8,'(16XA)')' NO.DISP. V ALUE'READ (6,70) (JPJ(I),PJ(I),I=1,NPJ)70 FORMAT (5(I5,G16.4))DO 80 I=1,NPJ80 WRITE(8,'(14X,I7,F16.3)') JPJ(I),PJ(I)90 CONTINUEIF(NPF.EQ.0) GOTO 130WRITE(8,'(/20X,16HNON-NODEL LOADS )')WRITE(8,'(11X,A,8X,A,9X,A)')'NO.E NO.LOAD.MODEL','A','C'READ (6,100) (JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)100 FORMAT (2(2I5,2G16.4))DO 110 I=1,NPF110 WRITE(8,120) (JPF(J,I),J=1,2),PF(1,I),PF(2,I)120 FORMAT(6X,2I8,10X,2F10.3)130 CONTINUERETURNENDSUBROUTINE MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)REAL KE(6,6),X(100),Y(100),EAI(3,100),AL(100),LINTEGER JE(2,100),JEAI(100)II=JE(1,IE)JJ=JE(2,IE)MT=JEAI(IE)L=SQRT((X(JJ)-X(II))**2+(Y(JJ)-Y(II))**2)AL(IE)=LA1=EAI(1,MT)*EAI(2,MT)/LA2=EAI(1,MT)*EAI(3,MT)/L**3A3=EAI(1,MT)*EAI(3,MT)/L**2A4=EAI(1,MT)*EAI(3,MT)/LKE(1,1)=A1KE(1,4)=-A1KE(2,2)=12*A2KE(2,3)=6*A3KE(2,5)=-12*A2KE(2,6)=6*A3KE(3,3)=4*A4KE(3,5)=-6*A3KE(3,6)=2*A4KE(4,4)=A1KE(5,5)=12*A2KE(5,6)=-6*A3KE(6,6)=4*A4DO 10 I=1 ,6DO 10 K=I ,610 KE(K,I)=KE(I,K)RETURNENDSUBROUTINE MR(R,IE,JE,X,Y)REAL R(6,6),X(100),Y(100),L,CX,CYINTEGER JE(2,100)I=JE(1,IE)J=JE(2,IE)L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)CX=(X(J)-X(I))/LCY=(Y(J)-Y(I))/LDO 10 J=1,6DO 10 I=1,610 R(I,J)=0.DO 20 I=1,4,3R(I,I)=CXR(I,I+1)=CYR(I+1,I)=-CYR(I+1,I+1)=CX20 R(I+2,I+2)=1.RETURNENDSUBROUTINE MAKE(KE,R,AKE)REAL KE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)CALL TRAN(R,RT)CALL MULV(RT,KE,TMP)CALL MULV(TMP,R,AKE)RETURNENDSUBROUTINE CALM(M,IE,JN,JE)INTEGER M(6),JN(3,100),JE(2,100),IEDO 10 I=1,3M(I)=JN(I,JE(1,IE))10 M(I+3)=JN(I,JE(2,IE))RETURNENDSUBROUTINE MK(K,AKE,M)REAL K(200,200),AKE(6,6)INTEGER M(6)DO 10 I=1,6DO 10 J=1,6IF(M(I).NE.0.AND.M(J).NE.0)& K(M(I),M(J))=K(M(I),M(J))+AKE(I,J)10 CONTINUERETURNENDSUBROUTINE PE(FE,IP,JPF,PF,AL)REAL FE(6),PF(2,100),AL(100),LINTEGER JPF(2,100)A=PF(1,IP)C=PF(2,IP)L=AL(JPF(1,IP))IND=JPF(2,IP)DO 5 I=1,65 FE(I)=0.GOTO(10,20,30,40,50,60),IND10 FE(2)=(7*A/20+3*C/20)*LFE(3)=(A/20+C/30)*L**2FE(5)=(3*A/20+7*C/20)*LFE(6)=-(A/30+C/20)*L**2RETURN20 FE(5)=A*C**3*(2*L-C)/2/L**3FE(2)=A*C-FE(5)FE(3)=A*C**2*(6*L*L-8*C*L+3*C*C)/12/L/LFE(6)=-A*C**3*(4*L-3*C)/12/L/LRETURN30 FE(2)=A*(L-C)**2*(L+2*C)/L**3FE(3)=A*C*(L-C)**2/L**2FE(5)=A-FE(2)FE(6)=-A*C**2*(L-C)/L**2RETURN40 FE(2)=-6*A*C*(L-C)/L**3FE(3)=A*(L-C)*(L-3*C)/L**2FE(5)=-FE(2)FE(6)=A*C*(3*C-2*L)/L**2RETURN50 FE(1)=A*(1-C/L)FE(4)=A*C/LRETURN60 FE(1)=C*L/2.FE(4)=FE(1)RETURNENDSUBROUTINE MULV6(A,B,C)REAL C(6),A(6,6),B(6)DO 10 I=1,6C(I)=0.DO 10 J=1,610 C(I)=C(I)+A(I,J)*B(J)RETURNENDSUBROUTINE MF(P,AFE,M)REAL P(100),AFE(6)INTEGER M(6)DO 10 I=1,6IF(M(I).NE.0)P(M(I))=AFE(I)+P(M(I))10 CONTINUERETURNENDSUBROUTINE SLOV(AK,P,D,N)REAL AK(200,200),P(100),D(100)DO 5 I=1,1005 D(I)=P(I)DO 10 K=1,N-1DO 10 I=K+1,NC=-AK(K,I)/AK(K,K)DO 20 J=I,N20 AK(I,J)=AK(I,J)+C*AK(K,J)10 D(I)=D(I)+C*D(K)D(N)=D(N)/AK(N,N)DO 40 I=N-1,1,-1DO 30 J=I+1,N30 D(I)=D(I)-AK(I,J)*D(J)40 D(I)=D(I)/AK(I,I)RETURNENDSUBROUTINE MADE(IE,JN,JE,D,ADE)REAL ADE(6),D(100)INTEGER IE,JN(3,100),JE(2,100)DO 3 I=1,63 ADE(I)=0DO 10 I=1,3IF (JN(I,JE(1,IE)).NE.0) ADE(I)=D(JN(I,JE(1,IE)))IF (JN(I,JE(2,IE)).NE.0) ADE(I+3)=D(JN(I,JE(2,IE)))10 CONTINUERETURNENDSUBROUTINE TRAN(R,R T)REAL R(6,6),RT(6,6)DO 10 I=1,6DO 10 J=1,610 RT(I,J)=R(J,I)RETURNENDSUBROUTINE MULV(A,B,C)REAL A(6,6),B(6,6),C(6,6)DO 10 I=1 , 6DO 10 J=1 , 6C(I,J)=0.DO 10 K=1 , 610 C(I,J)=C(I,J)+A(I,K)*B(K,J)RETURNEND三．课后习题（1）准备原始数据○1确定结点、划分单元，建立整体坐标系与局部坐标系如下图所示。