有限元基础知识-结构作业

一、PF程序算例计算图1所示结构，两杆均为工字型钢，截面面积A=102cm2，截面二次矩输入文件 ASimpleSuppoted E2331AI12102E-32240E-23102E-32240E-XY004383000节点荷载 非节点荷载 12-224-4输出文件 ASimpleSuppoted TheInputDataTeral TheInformationof .102000E- .322400E- .102000E- .322400E-TheJoint TheInformationofSupports Loading TheLoadingsatJoints TheLoadingsatMembers - TheResultsofCalculationTheJointDisplacementsuv12.731849E--1.000000E--1.574435E-23.235083E--4.336787E--9.741742E- 3.235083E--1.000000E-1.871939E-TheM 12-- 2-3源程序!PF.FOR(Aprogramforysisofplane!Version !Mainprogramreadsthecontroldata anizesthe!calculationbycallingCHARACTERIDFN*12,OUTFN*12,TITLE(5)*72WRITE(*,1) FORMAT(1X,'InputDataFileName:')WRITE(*,2)FORMAT(1X,'OutputFileName:')FORMAT(1X,'InputDataFileName:',1X,A12)WRITE(4,4)OUTFNFORMAT(5X,'OutputFileName:',1X,A12)FORMAT(/13X,'TheInputData'//10X,'TeralFORMAT(/10X,'TheInformationofCALLSTART(E,NM,NJ,NS,NLC)SUBROUTINEDIMENSIONDIMENSIONIU(3,NJ),LV(6,NM),LD(90),LJ(50),FX(50),FY(50),FM(50)DIMENSIONLM(50),LT(50),VF(50),DST(50)DOUBLEPRECISIONAE(6,6),A(300),B(90),F(6,NM)CALLIOMJS(NM,NJ,NS,NLC,IST,IEN,AR,RI,X,Y,IS,VS)CALLIDUN(NJ,N,IU,X,Y)CALLLCVCTCALLLCDIACALLFORMA(E,NM,NJ,N,NA,IST,IEN,AR,RI,X,Y,LV,AE,LD,A)CALLAS(NS,NJ,N,NA,IS,IU,LD,A)CALLLDLT(N,NA,LD,A,B)DO100LC=1,NLCREADCALLB0(LIF(NLJ.NE.0)CALLIOLJB(N,NJ,NLJ,IU,LJ,FX,FY,FM,B)CALLIF(NLM.NE.0)CALLCALLBS(NS,NJ,N,IS,VS,IU,B)CALLSLVEQ(N,NA,MAXBDW,LD,A,B)CALLOJD(NJ,N,IU,B)CALL100SUBROUTINE Readdataofmembers,joints,supports&printDIMENSION FORMAT(/10X,'TheJointCoordinates'//1X,'joint',8X,'X',15X,'Y') FORMAT(/10X,'TheInformationofSupports'//4X,'IS',8X,'VS')READ(3,*)(IS(M),VS(M),M=1,NS) SUBROUTINE DeterminingthenumbersofindependentunknownsDIMENSIONIU(3,NJ),X(NJ),Y(NJ)DO200DO100I=J-1,1,-1DX=ABS(X(J)-DY=ABS(Y(J)-IF(DX.LT.1E-4.AND.DY.LT.1E-4)THENGOTO150ENDIF 150200SUBROUTINE DeterminelocationvectorofelementDIMENSIONIST(NM),IEN(NM),LV(6,NM),IU(3,NJ)DO100M=1,NM100SUBROUTINE Determinelocationofdiagonalelementsofglobal matrixDIMENSIONLV(6,NM),MIN(N),IBDW(N),LD(N)DO100I=1,N100DO400M=1,NMDO200I=2,6IF(LV(I,M).LT.MINLV) DO300IF(MINLV.LT.MIN(LV(I,M)))DO500IF(IBDW(I).GT.MAXBDW)MAXBDW=IBDW(I)500CONTINUEDO600I=2,N600CONTINUESUBROUTINE Calculaength,cosine&sineofmemberDIMENSIONIST(NM),IEN(NM),X(NJ),Y(NJ)DOUBLEPRECISIONRL,C,S,X1,Y1X1=X(J)-SUBROUTINEKEBAR!CalculateelementstiffnessmatrixalonglocalaxesDIMENSIONIST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM)DOUBLEPRECISIONC,S,E1,E2,E3,E4,RLCALLLCS(M,NM,NJ,IST,IEN,X,Y,RL,C,S)SUBROUTINEKE CalculateelementstiffnessmatrixalongglobalaxesDIMENSIONIST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ)DOUBLEPRECISIONAE(6,6),C,S,E1,E2,E3,E4,A1,A2,A3,A4,A5,A6CALLKEBAR(M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4)AE(4,1)=-AE(5,1)=-AE(5,2)=-SUBROUTINEFORMA FormglobalstiffnessmatrixDIMENSIONIST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&DOUBLEPRECISIONA(NA),AE(6,6)DO300CALLKE(M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)DO200I=1,6DO100IF(LV(I,M).GE.LV(J,M))A(LD(LV(I,M))-=A(LD(LV(I,M))-A(LD(LV(J,M))- 300CONTINUESUBROUTINEAS IntroducesupportconditionsintoglobalstiffnessmatrixADIMENSIONIS(NS),IU(3,NJ),LD(N)DOUBLEPRECISIONA(NA)DO100M=1,NS100CONTINUESUBROUTINELDLT Solveequations(1)- positionofmatrixADIMENSIONLD(N)DOUBLEPRECISIONA(NA),T(N),SUMDO300I=2,NDO200J=I1,I-1IF(J.GT.1)THENIF(I1.GT.J1)J1=I1DO100K=J1,J-1 T(J)=A(LDI-I+J)- 300CONTINUESUBROUTINESLVEQ Solveequations(2)-forward&backsubstitutionDIMENSIONLD(N)DOUBLEPRECISIONA(NA),B(N)DO200I=2,NDO100J=I1,I-1B(I)=B(I)-A(LDI- 200DO300I=1,N300DO500I=N-1,1,-1IF(IMIN.GT.N)IMIN=NDO400J=I+1,IMINJ1=J-LDJ+LD(J-IF(I.GE.J1)B(I)=B(I)-A(LDJ-400500CONTINUESUBROUTINEB0(L InitializejointloadvectorBDOUBLEPRECISIONB(N)WRITE FORMAT(/15X,'Loading//10X,'TheLoadingsatDO100I=1,N100CONTINUESUBROUTINEIOLJB Readdataofloadingatjoint&formjointloadvectorBDIMENSIONIU(3,NJ),LJ(NLJ),FX(NLJ),FY(NLJ),FM(NLJ)DOUBLEPRECISIONB(N) WRITE DO100M=1,NLJ100CONTINUESUBROUTINE InitializeterminalforcesofmembersDOUBLEPRECISIONF(6,NM)WRITE FORMAT(/10X,'TheLoadingsatDO200J=1,NMDO100I=1,6 200CONTINUESUBROUTINE Readdataofloadingatmember&calculatefixed-end addequivalentjointloadstovectorDIMENSIONLM(NLM),LT(NLM),VF(NLM),DST(NLM),IST(NM),&DOUBLEPRECISION WRITE DO100M=1,NLMCALLLCS(L,NM,NJ,IST,IEN,X,Y,RL,C,S)D2=RL-IF(LT(M).EQ.1.OR.LT(M).EQ.3)THENF4=-P1-B(LV(1,L))=B(LV(1,L))-B(LV(2,L))=B(LV(2,L))-F1*S-B(LV(3,L))=B(LV(3,L))-B(LV(4,L))=B(LV(4,L))-B(LV(5,L))=B(LV(5,L))-F4*S-B(LV(6,L))=B(LV(6,L))-100CONTINUESUBROUTINE IntroducesupportconditionsintojointloadvectorBDIMENSIONIS(NS),VS(NS),IU(3,NJ)DOUBLEPRECISIONB(N)DO100M=1,NS100CONTINUESUBROUTINE PrintjointdisplacementsDIMENSIONIU(3,NJ)DOUBLEPRECISIONB(N)WRITE(4,1)1FORMAT(/13X,'TheResultsof//10X,'TheJointWRITE(4,2)(J,B(IU(1,J)),B(IU(2,J)),B(IU(3,J)),J=1,NJ)2FORMAT(1X,I4,1P3E16.6)SUBROUTINE Calculate&printterminalforcesofmembersDIMENSIONIST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&DOUBLEPRECISIONWRITE1FORMAT(/10X,'TheTerminalDO100M=1,NMCALLKEBAR(M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4)F(2,M)=F(2,M)+E2*(U2-F(3,M)=F(3,M)+E3*(U2-F(4,M)=F(4,M)+E1*(U4-F(5,M)=F(5,M)+E2*(U5-U2)-WRITE(4,2)M,IST(M),F(1,M),F(2,M),F(3,M),& 100二、三角形单元程序算例h=10mm。现用有限元法分析其位移及应力。输入文件232ETXYXYXYXYXY00000000000-000000000-000000000-000000000-000000000-000000000-000000000-000000000-000000000-000000000-000000000-00000000单元节点约束节点的位移度序号输出文件***********TheinputET2322Q-----------172283394456***********Theresau----------3.0568E--------------------------------------------------------1.1911E-------------5.8304E--1.8891E-1.5068E-1.4920E--5.0731E-1.0190E--1.6415E-9.0267E--5.4597E--3.8697E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-2.0460E--5.7800E--8.6660E-----源程序 REWIND7WRITE(9,*)'***********Theinputdata***********'600FORMAT(10X,'NFNENNMBNDLELSEUMT'/7(2X,I4),3(2X,F8.4))CALLSUBROUTINE)CALLINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)CALLSOLVE(SK,Q,MS,MB)WRITE(9,*)'***********Theresauata***********'200FORMAT(5X,'DISPLACEMENT')250CALL SUBROUTINEINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSIONXY(LS,NN),Q(MS),NX(NE,LE),MC(ND)WRITE(9,600)LS,NN600FORMAT(10X,'LSNN'/2(2X,I4))DO35DO35 DO DO DO DO CALLSTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)IF(L.EQ.1)WRITE(*,70)EK DO DO200 DO200 DO200 200DO220 220SK(M,1)=SK(M,1)*1E8 DO125IF(K+MB-1-N=K+MB-GOTODO125 DO122 122SK(I,J)=SK(I,J)-125Q(I)=Q(I)-DO145 IF(MS-I+1-135N=MS-GOTO136140DO142L=J+I-142Q(I)=Q(I)-145Q(I)=Q(I)/SK(I,1)147FORMAT(5X,'DISPLACEMENT')150FORMAT(2X,6E14.5)SUBROUTINESTRES(E,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T) 160QE(N)=Q(N1)165170 200CALL400SUBROUTINESTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T) B(2)=XYE(2,3)-B(3)=XYE(2,1)-C(1)=XYE(1,3)-C(2)=XYE(1,1)-C(3)=XYE(1,2)- DO30J=1,3N=2*(J-EK(M,N)=D*(B(I)*B(J)+C(I)*C(J)*(1-EK(M,N+1)=D*(UM*B(I)*C(J)+C(I)*B(J)*(1-EK(M+1,N)=D*(UM*C(I)*B(J)+B(I)*C(J)*(1-30EK(M+1,N+1)=D*(C(I)*C(J)+B(I)*B(J)*(1-UM)/2) B(2)=XYE(2,3)-B(3)=XYE(2,1)-C(1)=XYE(1,3)-C(2)=XYE(1,1)-C(3)=XYE(1,2)- 220S(3,2*I)=D*(1- 250SG(I)=SG(I)+S(I,K)*QE(K)300FORMAT(5X,'SEGMA',3E20.4)三：矩形薄板弯曲的求3q的作用，已知长度a1000mmb500mmt20mmq1Nmm2，E2.1105Nmm2，泊松比0.3现利用有限单元法求解矩形w及弯矩Mx、My的分布曲线。3模OABC矩形网格，其节点及单元也如图6所示。建立矩形薄板的离散化计算模型。按照图示把矩形薄板划分为均匀的四节点的矩形薄板单元网格。模型一共有5630个节点，每个节点 度数目为90，任意节点i单元、节点的按照自右向左、自上向下的规则进行。AOy0；OCx0ABw0，切线转角y0，法线转角x0；BCw0，切线转角y0；根据上述的边界条件，在结构的总度中位移为零的度序号为：1、2、、、 、、 、、、1、32、33、44、46、47、48、59、、、、、、、、、82、83、85、86、88、8935个度被固定

Myk

Myl

Mym

MynZkZlZmZnqABMxk

1qAB23MxmMxnMxkMykM

1qA2Bm3mmMykMynMykmmem加 FORTRANPROGRAMMAINOPEN(7,FILE='INP.DAT')REWIND7READ(7,*)NF,NE,NN,MB,ND,LE,LSREAD(7,*)E,UM,T FORMAT FORMATWRITE(*,600)NF,NE,NN,MB,ND,LE,LS,E,UM,TCALLINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND) WRITE(*,101) FORMATWRITE(*,500)((NX(I,J),I=1,NE),J=1,LE)WRITE(*,400)(MC(I),I=1,ND)FORMAT(10X,'NFNENNMBNDLELSEUMCALLSOLVE(SK,Q,MS,MB)OPEN(9,FILE='OUT.DAT')REWIND9WRITEWRITE CALLSTRES(E,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)STOP1000SUBROUTINEINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSIONXY(LS,NN),Q(MS),NX(NE,LE),MC(ND) SUBROUTINESTIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,&DO35I=1,MSDO35 DO200DO40J=1,NEDO40I=1,LS DO50DO50 CALLSTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)IF(L.EQ.1)WRITE(*,70)EK DO200DO200II=1,NFDO200J=1,NEDO200JJ=1,NFIF(MN) DO220I=1,ND SUBROUTINESOLVE(SK,Q,MS,MB)DIMENSIONSK(MS,MB),Q(MS)DO125IF(K+MB-1-MS) N=K+MB-GOTO DO125I=I1,NDO122J=1,J1 SK(I,J)=SK(I,J)- Q(I)=Q(I)-DO145I1=1,MIF(MS-I+1-MB) N=MS-GOTO DO142L=J+I- Q(I)=Q(I)- SUBROUTINESTRES(E,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)DIMENSIONQ(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS,NE)DO400DO160DO160J=1,NF WRITE(*,165)WRITE(*,170) FORMAT FORMATDO200J=1,NEDO200 CALL SUBROUTINESTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)DIMENSIONEK(ME,ME),XYE(LS,NE)A=(XYE(1,2)-EK(6,5)=-EK(8,1)=-EK(8,4)=-EK(8,7)=-EK(9,8)=-EK(10,3)=-EK(10,6)=-EK(10,8)=-EK(10,9)=-DO10I=1,11DO10 SUBROUTINESTE(XYE,QE,NE,LS,ME,E,UM,T)DIMENSIONQE(ME),XYE(LS,NE)B=(XYE(2,4)-CALLSTR(QE,ME,A,B,E,UM,-A,-CALLSTR(QE,ME,A,B,E,UM,A,-CALLCALLSTR(QE,ME,A,B,E,UM,-A,B,T)SUBROUTINEDIMENSIONQE(ME),S(3,12),SG(3)C2=1.-A2=(1.-S(3,9)=-S(3,1