平面杆系结构有限元编制程序

1平面杆系结构有限元编制程序——斜腿刚构计算一、程序编制说明本程序全部采用VisualBasic语言编写，可以进行简单杆系结构的受力分析，能够精确计算出结构体系的结点位移和内力，并采用文本方式输入和输出数据，方便实用。例子是针对斜腿刚构的计算，全桥共划分为6个单元，计算简图如下。二、子程序划分DATA：数据的准备与输入ZG：总体刚体矩阵DG：单元刚度矩阵TRANSLATE：坐标转换矩阵RESULT：杆端位移与杆端力计算MULT：矩阵相乘MULT1：矩阵与向量相乘GAS：高斯消元法BOUNDARY：边界条件LOAD：外荷载三、标识符说明：1、整数型变量NN：结点总个数ELE：单元总个数GD：固定支座个数JZ：铰支座个数ZC：支承个数MS:系统总自由度数目（总刚阶数）PN：荷载作用的结点数（结点荷载个数）2、一维数组变量EI：杆件抗弯刚度EA：杆件抗拉刚度Cn：cos值Sn：sin值LN：单元左端结点标号数组RN：单元右端结点标号数组FX：水平方向结点外力FY：垂直方向结点外力MOMENT：结点弯矩Z：总体坐标系下结点位移存储P：总体坐标系下结点力存储3、二维数组变量BKE：局部坐标单元刚度矩阵存储TKE：体系总体刚度矩阵存储T：坐标转换矩阵存储四、源程序代码DimI,JAsDouble2DimNNAsInteger,ELEAsInteger,MSAsIntegerDimGDAsInteger,JZAsInteger,ZCAsInteger,PNAsIntegerDimLRN(1To12,1To12)AsDouble,L(1To7)AsDouble,G(1To21,1To21)AsDoubleDimEA(1To7)AsDouble,EI(1To7)AsDouble,CnAsDouble,SnAsDoubleDimX(1To7)AsDouble,Y(1To7)AsDouble,P(1To21)AsDoubleDimTKE(1To21,1To21)AsDouble,TKE1(1To21,1To21)AsDoubleDimBKE(1To6,1To6)AsDouble,T(1To6,1To6)AsDoubleDimTZ(1To6,1To6)AsDouble,Z(1To21)AsDoubleDimLN(1To6)AsInteger,RN(1To6)AsIntegerDimFX(1To7)AsDouble,FY(1To7)AsDouble,MOMENT(1To7)AsDoublePrivateSubForm_Click()Opene:\homework\gginput.txtForInputAs#1Opene:\homework\ggoutput.txtForOutputAs#2CallData'输入数据CallZG'总刚度矩阵CallBoundary'输入边界条件Callload'输入荷载CallResult'计算结果（杆端力）Close#1Close#2EndSubPublicSubData()'数据的准备与输入Print#2,平面杆系有限元程序——斜腿刚构计算Print#2,=========================================Print#2,原始数据的输入Print#2,Print#2,结点数;Spc(3);单元数;Spc(3);固定支座数;Spc(3);铰支座数;Spc(3);支承数;Spc(3);结点荷载数;Spc(3);系统总自由度数Input#1,NN,GD,JZ,ZC,PNMS=3*NN:ELE=NN-1Print#2,NN;Spc(6);ELE;Spc(6);GD;Spc(6);JZ;Spc(6);ZC;Spc(6);PN;Spc(6);MSPrint#2,Print#2,=========================================Print#2,结点坐标的输入Print#2,Print#2,结点;Spc(3);X;Spc(6);Y;Spc(3);单位：米ForI=1ToNNInput#1,I,X(I),Y(I)Print#2,I;Spc(3);X(I);Spc(6);Y(I)NextIPrint#2,Print#2,=========================================Print#2,杆件相关信息的输入Print#2,Print#2,杆件号;Spc(3);左结点号;Spc(3);右结点号;Spc(3);抗弯刚度;Spc(3);抗拉刚度;Spc(3);杆件长度;Spc(3);E=2.06E+11N/m23ForI=1ToELEInput#1,I,LN(I),RN(I),EI(I),EA(I)NextIForI=1ToELEL(I)=Sqr((Y(RN(I))-Y(LN(I)))^2+(X(RN(I))-X(LN(I)))^2)Print#2,I;Spc(8);LN(I);Spc(8);RN(I);Spc(8);Format(EI(I),0.000E+00);Spc(8);Format(EA(I),0.000E+00);Spc(8);(Format(L(I),0.000))NextI'输入结点编号LRN(1,1)=3#:LRN(1,2)=3#:LRN(1,3)=4#LRN(2,1)=1#:LRN(2,2)=2#:LRN(2,3)=3#LRN(1,4)=5#:LRN(1,5)=7#:LRN(1,6)=6#LRN(2,4)=4#:LRN(2,5)=5#:LRN(2,6)=5#Print#2,Print#2,=========================================Print#2,输入结点荷载Print#2,Print#2,结点号;Spc(3);FX;Spc(4);FY;Spc(4);MOMENT;Spc(4);单位：NForI=1ToNNInput#1,I,FX(I),FY(I),MOMENT(I)NextIForI=1ToNNPrint#2,I;Spc(3);Format(FX(I),0.0E+);Spc(4);Format(FY(I),0.0E+);Spc(4);Format(MOMENT(I),0.0E+)NextIPrint#2,Print#2,=========================================Print#2,输入约束条件EndSubPublicSubDG(K)'形成单元刚度矩阵'ForK=1ToELEForI=1To6ForJ=1To6BKE(I,J)=0#NextJNextIBKE(1,1)=EA(K)/L(K)BKE(2,1)=0#BKE(2,2)=12*EI(K)/L(K)*L(K)*L(K)BKE(3,1)=0#BKE(3,2)=6*EI(K)/L(K)*L(K)BKE(3,3)=4*EI(K)/L(K)BKE(4,1)=-EA(K)/L(K)BKE(4,2)=0#BKE(4,3)=0#BKE(4,4)=EA(K)/L(K)4BKE(5,1)=0#BKE(5,2)=-12*EI(K)/L(K)*L(K)*L(K)BKE(5,3)=-6*EI(K)/L(K)*L(K)BKE(5,4)=0#BKE(5,5)=12*EI(K)/L(K)*L(K)*L(K)BKE(6,1)=0#BKE(6,2)=6*EI(K)/L(K)*L(K)BKE(6,3)=2*EI(K)/L(K)BKE(6,4)=0#BKE(6,5)=-6*EI(K)/L(K)*L(K)BKE(6,6)=4*EI(K)/L(K)ForI=1To6ForJ=1To6BKE(I,J)=BKE(J,I)NextJNextI'Print#2,K,单刚'ForI=1ToELE'ForJ=1ToELE'Print#2,(Format(BKE(I,J),0.000));Spc(2);'NextJ'Print#2,'NextI'NextKEndSubPublicSubTranslate(K)'形成坐标变换矩阵'ForK=1ToELECn=(X(RN(K))-X(LN(K)))/L(K)Sn=(Y(RN(K))-Y(LN(K)))/L(K)ForI=1To6ForJ=1To6T(I,J)=0#NextJNextIT(1,1)=CnT(1,2)=SnT(1,3)=0#T(2,1)=-SnT(2,2)=CnT(2,3)=0#T(3,1)=0#T(3,2)=0#T(3,3)=1#ForI=1To3ForJ=1To3T(I+3,J+3)=T(I,J)5NextJNextI'Print#2,K;单元,坐标变换矩阵'ForI=1ToELE'ForJ=1ToELE'Print#2,(Format(T(I,J),0.000));Spc(2);'NextJ'Print#2,'NextI'NextK'形成坐标变换矩阵的转置矩阵(逆矩阵）ForI=1To6ForJ=1To6TZ(I,J)=T(J,I)NextJNextIEndSubPublicSubZG()'形成总体刚度矩阵DimMAsInteger,M1AsInteger,HAsInteger,H1AsIntegerDimIAsInteger,JAsInteger,WAsInteger,IIAsInteger,JJAsIntegerForI=1ToMSForJ=1ToMSTKE(I,J)=0#NextJNextIForW=1ToELE'形成总体坐标的单元刚度矩阵CallTranslate(W)CallDG(W)CallMULT(BKE,T,G,6,6,6)CallMULT(TZ,G,TKE1,6,6,6)'形成总体刚度矩阵ForI=1To2ForII=1To3M=3*(I-1)+IIM1=3*(LRN(I,W)-1)+IForJ=1To2ForJJ=1To3H=3*(J-1)+JJH1=3*(LRN(J,W)-1)+JJTKE(M1,H1)=TKE(M1,H1)+TKE1(M,H)NextJJNextJNextIINextINextW6Print#2,总刚ForI=1ToMSForJ=1ToMSPrint#2,(Format(TKE(I,J),0.0E+00));Spc(2);NextJPrint#2,NextIEndSubPublicSubMULT(A,B,C,N1,N2,N3)'矩阵相乘运算DimA(N1,N2)AsDouble,B(N2,N3)AsDouble,C(N1,N3)AsDoubleForI=1ToN1ForJ=1ToN3C(I,J)=0#NextJNextIForI=1ToN1ForJ=1ToN3ForK=1ToN2C(I,J)=C(I,J)+A(I,K)*B(K,J)NextKNextJNextIEndSubPublicSubGAS(A,B,N,X)'高斯消元法DimKAsDouble,HAsDoubleForK=1ToN-1ForJ=K+1ToNA(K,J)=A(K,J)/A(K,K)NextJB(K)=B(K)/A(K,K)ForI=K+1ToNForJ=K+1ToNA(I,J)=A(I,J)-A(I,K)*A(K,J)NextJB(I)=B(I)-A(I,K)*B(K)NextINextKX(N)=B(N)/A(N,N)ForI=N-1To1Step-1H=0ForJ=I+1ToNH=H+A(I,J)*X(J)NextJX(I)=B(I)-HNextIEndSub7PublicSubload()'形成结点荷载向量ForI=1ToMSP(I)=0#NextIP(1)=0#P(2)=0#P(3)=0#P(4)=0#P(5)=0#P(6)=0#P(7)=0#P(8)=0#P(9)=0#P(10)=0#P(11)=-100000P(12)=0#P(13)=35000P(14)=0#P(15)=0#P(16)=0#P(17)=0#P(18)=0#P(19)=0#P(20)=0#