《有限单元法FORTRAN程序设计作业.doc》由会员分享，可在线阅读，更多相关《有限单元法FORTRAN程序设计作业.doc（44页珍藏版）》请在三一文库上搜索。
1、结构工程 104811085 彭晓丽引言有限元分析是从二十世纪发展起来的一种新的结构分析方法，其分析问题的主要思想是把连续弹性体离散成为一群仅在节点处相互连接的有限单元的集合，通过把无限连续问题有限化进行进一步的分析研究。有限元最早是美国波音公司工程师特纳采用三角形和矩形单元，把结构力学中的位移法应用于飞机结构的分析中；后来，人们逐渐认识到，有限元是一种基于虚功原理的广义里茨法；到二十世纪七十年代后，随着计算机技术的迅猛发展，人们采用各种成熟的单元编制了一系列结构分析程序，例如威尔逊教授等编制的SAP系列、拜塞等编制的ADINA系列、ANSYS公司研制的ANSYS系列等。通过本学期课程的学习，
2、我们进一步学习了Fortran90编制空间桁架的基本方法和理念，并对以下程序进行了调试，结合实例进行了计算。程序1：平面三节点有限元分析程序（PLANE.FOR）试用程序PLANE.FOR计算图2所示悬臂梁结构，其受均布荷载q=1N/mm2 作用，E=2.1105 N/mm2 ，=0.3，厚度h=10mm，有限元法分析其位移及应力。图2 悬臂梁单元划分及荷载简图 梁视为平面应力状态，按图2尺寸划分为均匀的三角形网格，共80个单元，55个节点，坐标轴及单元与节点的编号如图。将均布荷载分配到相应的节点上，把有约束的节点51、52、53、54、55视作固定铰链，建立如图2所示的离散化计算模型。源程序
3、：PROGRAM MAINDIMENSION SK(300,30),EK(12,12),Q(300),MC(55),XY(3,100),XYE(3,4),& QE(12),NX(4,100)OPEN (7,FILE=INPUT.TXT)REWIND 7READ (7,*) NF,NE,NN,MB,ND,LE,LSREAD (7,*) E,UM,T10 FORMAT(7I5)12 FORMAT(3F15.2)WRITE (*,600) NF,NE,NN,MB,ND,LE,LS,E,UM,TME=NE*NFMS=NN*NFCALL INPUT (XY,Q,NX,MC,LS,NN,MS,NE,LE,
4、ND)WRITE (*,102) (XY(I,J),I=1,LS),J=1,NN)102 FORMAT (10X,XY/,(2X,6F12.3)WRITE (*,101) (Q(I),I=1,MS)101 FORMAT (10X,Q/,(2X,6F12.3)WRITE (*,500) (NX(I,J),I=1,NE),J=1,LE)500 FORMAT (10X,NX/,(2X,12I6) WRITE (*,400) (MC(I),I=1,ND)600 FORMAT (10X,NF NE NN MB ND LE LS E UM T/7(2X,I4),3(2X,F8.4)400 FORMAT (
5、10X,MC/,(2X,10I6)CALL STIFS (SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,& NN,LS,E,UM,T)CALL SOLVE (SK,Q,MS,MB)OPEN (9,FILE=OUT.DAT)REWIND 9WRITE (9,200)WRITE (9,250) (Q(I),I=1,MS)200 FORMAT (5X,DISPLACEMENT)250 FORMAT (2X, 6E14.5)CALL STRES (Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)STOP 1000ENDSUBR
6、OUTINE INPUT (XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND)READ (7,*) XYREAD (7,*) QREAD (7,*) NXREAD (7,*) MCCLOSE(7)10 FORMAT(6F11.2)20 FORMAT(12I5)RETURNENDSUBROUTINE STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,& NN,LS,E,UM,T)DIMENSION SK(MS,MB),EK(ME,ME),Q(MS),NX(NE
7、,LE),MC(ND),& XY(LS,NN),XYE(LS,NE)DO 35 I=1,MSDO 35 J=1,MB35 SK (I,J)=0DO 200 L=1,LEDO 40 J=1,NELJ=NX(J,L)DO 40 I=1,LS40 XYE(I,J)=XY(I,LJ)DO 50 I=1,MEDO 50 J=1,ME50 EK(I,J)=0.0CALL STIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)IF(L.EQ.1) WRITE(*,70)EK70 FORMAT (10X,EK/,(6E14.5)DO 200 I=1,NEDO 200 II=1,NFM=NF*(I-1
8、)+IIM1=NF*(NX(I,L)-1)+IIDO 200 J =1,NEDO 200 JJ=1,NF N=NF*(J-1)+JJN1=NF*(NX(J,L)-1)+JJMN=N1-M1+1IF(MN) 200,200,150150 SK(M1,MN)=SK(M1,MN)+EK(M,N)200 CONTINUEDO 220 I=1,NDM=MC(I)220 SK(M,1)=SK(M,1)*1E8RETURNENDSUBROUTINE SOLVE(SK,Q,MS,MB)DIMENSION SK(MS,MB),Q(MS)K1=MS-1DO 125 K=1,K1IF(K+MB-1-MS) 105,
9、106,106105 N=K+MB-1GOTO 110106 N=MS110 I1=K+1DO 125 I=I1,NL=I-K+1C=SK(K,L)/SK(K,1)J1=MB-L+1DO 122 J=1,J1M=J+I-K122 SK(I,J)=SK(I,J)-C*SK(K,M)125 Q(I)=Q(I)-C*Q(K)Q(MS)=Q(MS)/SK(MS,1)M=MS-1DO 145 I1=1,MI=MS-I1IF(MS-I+1-MB) 135,136,136135 N=MS-I+1GOTO 140136 N=MB140 DO 142 J=2,NL=J+I-1142 Q(I)=Q(I)-SK(I
10、,J)*Q(L)145 Q(I)=Q(I)/SK(I,1)WRITE(*,147)147 FORMAT (5X, DISPLACEMENT )WRITE( *, 150) (Q(I) ,I=1 ,MS)150 FORMAT(2X, 6E14.5)RETURNENDSUBROUTINE STRES (Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)DIMENSION Q(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS,NE)DO 400 L=1,LEDO 160 I=1,NEDO 160 J=1,NFN=NF*(I-1)+JN1=
11、NF*(NX(I,L)-1)+J160 QE(N)=Q(N1)WRITE (*,165) LWRITE (*,170) (QE(I),I=1,ME)165 FORMAT (4X,L=,I4)170 FORMAT (6E14.5)DO 200 J=1,NELJ=NX(J,L)DO 200 I=1,LS200 XYE(I,J)=XY(I,LJ)CALL STE (XYE,QE,NE,LS,ME,E,UM,T)400 CONTINUERETURNENDSUBROUTINE STIFE (EK,XYE,ME,NE,NF,LS,E,UM,T)DIMENSION EK(ME,ME),XYE(LS,NE),
12、B(3),C(3)B(1)=XYE(2,2)-XYE(2,3)B(2)=XYE(2,3)-XYE(2,1)B(3)=XYE(2,1)-XYE(2,2)C(1)=XYE (1,3)-XYE (1,2)C(2)=XYE(1,1)-XYE(1,3)C(3)=XYE(1,2)-XYE(1,1)AE=(B(2)*C(3)-B(3)*C(2)/2D=E*T/4/(1-UM*2)/AEDO 30 I=1,3DO 30 J=1,3M=2*(I-1)+1N=2*(J-1)+1EK(M,N)=D*(B(I)*B(J)+C(I)*C(J)*(1-UM)/2)EK(M,N+1)=D*(UM*B(I)*C(J)+C(I
13、)*B(J)*(1-UM)/2)EK(M+1,N)=D*(UM*C(1)*B(J)+B(I)*C(J)*(1-UM)/2)30 EK(M+1,N+1)=D*(C(I)*C(J)+B(I)*B(J)*(1-UM)/2)RETURNENDSUBROUTINE STE (XYE,QE,NE,LS,ME,E,UM,T)DIMENSION XYE(LS,NE),QE(ME),SG(3),B(3),C(3),S(3,6)B(1)=XYE(2,2)-XYE(2,3)B(2)=XYE(2,3)-XYE(2,1)B(3)=XYE(2,1)-XYE(2,2)C(1)=XYE(1,3)-XYE(1,2)C(2)=X
14、YE(1,1)-XYE(1,3)C(3)=XYE(1,2)-XYE(1,1)AE=(B(2)*C(3)-B(3)*C(2)/2D=E/2/(1-UM*2)/AEDO 220 I=1,3S(1,2*I-1)=D*B(I)S(2,2*I-1)=D*UM*B(I)S(3,2*I-1)=D*(1-UM)*C(I)/2S(1,2*I)=D*UM*C(I)S(2,2*I)=D*C(I)220 S(3,2*I)=D*(1-UM)*B(I)/2DO 250 I=1,3SG(I)=0DO 250 K=1,ME250 SG(I)=SG(I)+S(I,K)*QE(K)WRITE (9,300) SG300 FORM
15、AT (5X,SEGMA,3E20.4)RETURNEND输入文本：输出文本：程序2：平面刚架静力分析程序（PF.FOR） 试用程序PF.FOR计算图1所示结构。梁和柱均用45a号工字钢，截面面积A=102cm2，截面二次距I=32240cm4，弹性模量为E=2.1108kPa。图1 平面钢架荷载简图源程序：! PF.FOR (A program for the analysis of plane frame)! Version 6.3 2004! MAIN PROGRAM ! Reading and printing the control data, allocating the arra
16、ys ! and organizing the whole calculation by calling subroutines DIMENSION W (80000)CHARACTER IDFN*12,OUTFN*12,TITLE(5)*72WRITE(*,1)1 FORMAT (1X,Input Data File Name:)READ (*,(A12) IDFNOPEN (3,FILE=IDFN,STATUS=OLD)WRITE (*,2)2 FORMAT (1X,Output File Name:)READ (*,(A12)OUTFNOPEN (4,FILE=OUTFN,STATUS=
17、NEW)WRITE (4,3) IDFN3 FORMAT(1X,Input Data File Name:,1X,A12)WRITE (4,4) OUTFN4 FORMAT(5X,Output File Name:,1X,A12)READ (3,(A72) (TITLE(M),M=1,5)WRITE (4,5) (TITLE(M),M=1,5)5 FORMAT(1X,A72)WRITE(4,6)6 FORMAT(/13X,The Input Data/10X,The General& Information/6X,E,9X,NM,5X,NJ,5X,NS,5X,NLC)READ(3,*)E,NM
18、,NJ,NS,NLC WRITE(4,7)E,NM,NJ,NS,NLC7 FORMAT(1X,1PE10.3,4I7)WRITE(4,8)8 FORMAT(/10X,The Information of Members &9 /1X,member,2X,start,2X,end,9X,A,15X,I)L1=1L2=L1+NML3=L2+NML4=L3+NML11=L4+NML12=L11+NJL21=L12+NJL22=L21+NSL31=L22+NSL32=L31+3*NJL41=L32+6*NMCALL IOMJS (NM,NJ,NS,NLC,W(L1),W(L2),W(L3),W(L4)
19、, & W(L11),W(L12),W(L21),W(L22)CALL IDUM (NJ,N,W(L31),W(L11),W(L12)CALL LCVCT (NM,W(L1),W(L2),W(L32),NJ,W(L31)CALL LCDIA (NM,N,W(L32),W(L41),W(L41),W(L41),MAXBDW,NA)L51=L41+NL52=L51+72L53=L52+NA*2L54=L53L61=L54+N*2CALL FORMA (E,NM,NJ,N,NA,W(L1),W(L2),W(L3),W(L4), & W(L11),W(L12),W(L32),W(L51),W(L41)
20、,W(L52)CALL AS(NS,NJ,N,NA,W(L21),W(L31),W(L41),W(L52)CALL LDLT(N,NA,W(L41),W(L52),W(L53) DO 100 LC=1,NLC READ (3,*)NLJ L62=L61+NLJ L63=L62+NLJ L64=L63+NLJ L71=L61 L81=L71+6*NM*2 CALL B0(LC,N,NLJ,W(L54) IF(NLJ.NE.0) CALL IOLJB(N,NJ,NLJ,W(L31),W(L61),W(L62), & W(L63),W(L64),W(L54) READ(3,*)NLM L82=L81
21、+NLM L83=L82+NLM L84=L83+NLM CALL F0(NLM,NM,W(L71) IF(NLM.NE.0) CALL IOLMFB(NM,NJ,N,NLM,W(L81),W(L82), & W(L83),W(L84),W(L1),W(L2),W(L11),W(L12),W(L32),W(L71),W(L54) CALL BS(NS,NJ,N,W(L21),W(L22),W(L31),W(L54) CALL SLVEQ(N,NA,MAXBDW,W(L41),W(L52),W(L54) CALL OJD(NJ,N,W(L31),W(L54) CALL COTF(E,NM,NJ,
22、N,W(L1),W(L2),W(L3),W(L4),W(L11), & W(L12),W(L32),W(L54),W(L71)100 CONTINUEEND SUBROUTINE IOMJS(NM,NJ,NS,NLC,IST,IEN,AR,RI,X,Y,IS,VS)! Reading and printing the data of members, joints and supports DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ), & Y(NJ),IS(NS),VS(NS)READ(3,*) (IST(M),IEN(M),AR(M),RI(M
23、),M=1,NM)WRITE(4,1) (M,IST(M),IEN(M),AR(M),RI(M),M=1,NM)1FORMAT (1X,I4,I8,I6,1P2E16.6)WRITE(4,2)2FORMAT (/10X,The Joint Coordinates/1X,joint,8X,X,15X,Y)READ(3,*) (X(M),Y(M),M=1,NJ)WRITE(4,3) (M,X(M),Y(M),M=1,NJ)3FORMAT (1X,I4,2F16.6) WRITE (4,4)4FORMAT (/10X,The Information of Supports/4X,IS,8X,VS)R
24、EAD(3,*) (IS(M),VS(M),M=1,NS)WRITE(4,5) (IS(M),VS(M),M=1,NS)5FORMAT (1X,I5,F15.6)RETURNENDSUBROUTINE IDUN (NJ,N,IU,X,Y)! Determining the numbers of independent unknowns DIMENSION IU(3,NJ),X(NJ),Y(NJ)IU(1,1)=1IU(2,1)=2IU(3,1)=3N=3DO 200 J=2,NJ DO 100 I=J-1,1,-1 DX=ABS(X(J)-X(I) DY=ABS(Y(J)-Y(I) IF(DX
25、.LT.1E-4.AND.DY.LT.1E-4)THEN IU(1,J)=IU(1,I) IU(2,J)=IU(2,I) IU(3,J)=N+1 N=IU(3,J) GOTO 150 END IF100 CONTINUE IU(1,J)=N+1 IU(2,J)=N+2 IU(3,J)=N+3 N=IU(3,J)150 CONTINUE200 CONTINUERETURNENDSUBROUTINE LCVCT(NM,IST,IEN,LV,NJ,IU)! Determine the location vector of elementDIMENSION IST(NM),IEN(NM),LV(6,N
26、M),IU(3,NJ)DO 100 M=1,NM I=IST(M) J=IEN(M) LV(1,M)=IU(1,I) LV(2,M)=IU(2,I) LV(3,M)=IU(3,I) LV(4,M)=IU(1,J) LV(5,M)=IU(2,J) LV(6,M)=IU(3,J)100 CONTINUERETURNENDSUBROUTINE LCDIA(NM,N,LV,MIN,IBDW,LD,MAXBDW,NA)! Determine the location of diagonal elements of global stiffness! matrix ADIMENSION LV(6,NM),
27、MIN(N),IBDW(N),LD(N)DO 100 I=1,N MIN(I)=I100CONTINUEDO 400 M=1,NM MINLV=LV(1,M) DO 200 I=2,6 IF (LV(I,M).LT.MINLV) MINLV=LV(I,M)200 CONTINUE DO 300 I=1,6 IF(MINLV.LT.MIN(LV(I,M) MIN(LV(I,M)=MINLV300 CONTINUE400CONTINUEMAXBDW=0DO 500 I=1,N IBDW(I)=I-MIN(I)+1 IF(IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I)500CON
28、TINUELD(1)=IBDW(1)DO 600 I=2,N LD(I)=LD(I-1)+IBDW(I)600CONTINUENA=LD(N)RETURNENDSUBROUTINE LCS(M,NM,NJ,IST,IEN,X,Y,RL,C,S)! Calculating length, cosine and sine of memberDIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ)DOUBLE PRECISION RL,C,S,X1,Y1I=IST(M)J=IEN(M) X1=X(J)-X(I) Y1=Y(J)-Y(I)RL=SQRT(X1*X1+Y1*Y1)C=X
29、1/RL S=Y1/RLRETURNEND SUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4)! Calculating element stiffness matrix in local coordinatesDIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM)DOUBLE PRECISION RL,C,S,E1,E2,E3,E4CALL LCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)E1=E*AR(M)/RLE2=12.0*E*RI(M)/(RL*R
30、L*RL)E3=0.5*E2*RLE4=4.0*E*RI(M)/RLRETURNENDSUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)! Calculating element stiffness matrix in global coordinatesDIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ)DOUBLE PRECISION AE(6,6),C,S,E1,E2,E3,E4,A1,A2,A3,A4,A5,A6CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,
31、S,E1,E2,E3,E4)A1=E1*C*C+E2*S*SA2=(E1-E2)*C*SA3=E1*S*S+E2*C*CA4=E3*SA5=E3*CA6=E4AE(1,1)=A1AE(2,1)=A2AE(2,2)=A3AE(3,1)=-A4AE(3,2)=A5AE(3,3)=A6AE(4,1)=-A1AE(4,2)=-A2AE(4,3)=A4AE(4,4)=A1AE(5,1)=-A2AE(5,2)=-A3AE(5,3)=-A5AE(5,4)=A2AE(5,5)=A3AE(6,1)=-A4AE(6,2)=A5AE(6,3)=0.5*A6AE(6,4)=A4AE(6,5)=-A5AE(6,6)=A
32、6RETURNENDSUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,X,Y,LV,AE,LD,A)! Forming theglobal stiffness matrix ADIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ), & LV(6,NM),LD(N)DOUBLE PRECISION A(NA),AE(6,6)DO 300 M=1,NM CALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) DO 200 I=1,6 DO 100 J=1,I IF (LV(I,M).
33、GE.LV(J,M) THEN A(LD(LV(I,M)-LV(I,M)+LV(J,M) & =A(LD(LV(I,M)-LV(I,M)+LV(J,M)+AE(I,J) ELSE A(LD(LV(J,M)-LV(J,M)+LV(I,M) & =A(LD(LV(J,M)-LV(J,M)+LV(I,M)+AE(I,J)END IF100 CONTINUE200 CONTINUE300CONTINUERETURNENDSUBROUTINE AS (NS,NJ,N,NA,IS,IU,LD,A)! Introducing support conditions into the global stiffn
34、ess matrix ADIMENSION IS(NS),IU(3,NJ),LD(N)DOUBLE PRECISION A(NA)DO 100 M=1,NS J=IS(M)/10 I=MOD(IS(M),10) K=IU(I,J) A(LD(K)=1D22100CONTINUERETURNENDSUBROUTINE LDLT (N,NA,LD,A,T)! Solve equations (1) - decomposition of matrix ADIMENSION LD(N)DOUBLE PRECISION A(NA),T(N),SUMDO 300 I=2,N LDI=LD(I) I1=I-
35、LDI+LD(I-1)+1 DO 200 J=I1,I-1 LDJ=LD(J) J1=J-LDJ+LD(J-1)+1 IF(I1.GT.J1) J1=I1 SUM=0.0D0 DO 100 K=J1,J-1 SUM=SUM+T(K)*A(LDJ-J+K)100 CONTINUE T(J)=A(LDI-I+J)-SUM A(LDI-I+J)=T(J)/A(LDJ) A(LDI)=A(LDI)-T(J)*A(LDI-I+J)200 CONTINUE300CONTINUERETURNENDSUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)! Solve equations (
36、2) - forward and back substitutionDIMENSION LD(N)DOUBLE PRECISION A(NA),B(N)DO 200 I=2,N LDI=LD(I) I1=I-LDI+LD(I-1)+1 DO 100 J=I1,I-1 B(I)=B(I)-A(LDI-I+J)*B(J)100 CONTINUE200CONTINUEDO 300 I=1,N B(I)=B(I)/A(LD(I)300CONTINUEDO 500 I=N-1,1,-1 IMIN=I+MAXBDW IF(IMIN.GT.N) IMIN=N DO 400 J=I+1,IMIN LDJ=LD
37、(J) J1=J-LDJ+LD(J-1)+1 IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J)400 CONTINUE500CONTINUERETURNENDSUBROUTINE B0 (LC,N,NLJ,B)! Initializing joint load vector B to zeroDOUBLE PRECISION B(N) WRITE (4,1)LC,NLJ1FORMAT(/15X,Loading Case,I3/10X,The Loadings at Joints/17X,NLJ=,I4)DO 100 I=1,N B(I)=0.0D0100CONTINU
38、E RETURNENDSUBROUTINE IOLJB (N,NJ,NLJ,IU,LJ,FX,FY,FM,B)! Reading and printing the data of loadings at ! joints and forming joint load vector BDIMENSIONIU(3,NJ),LJ(NLJ),FX(NLJ),FY(NLJ),FM(NLJ)DOUBLE PRECISION B(N)WRITE (4,1)1 FORMAT(/1X,joint,8X,FX,14X,FY,16X,FM) READ (3,*)(LJ(M),FX(M),FY(M),FM(M),M=
39、1,NLJ) WRITE (4,2)(LJ(M),FX(M),FY(M),FM(M),M=1,NLJ)2 FORMAT(1X,I4,2F16.6,F18.6)DO 100 M=1,NLJ J=LJ(M) B(IU(1,J)=FX(M) B(IU(2,J)=FY(M) B(IU(3,J)=FM(M)100CONTINUERETURN ENDSUBROUTINE F0(NLM,NM,F)! Initializing the terminal forces of members to zeroDOUBLE PRECISION F(6,NM)WRITE (4,1)NLM1 FORMAT(/10X,Th
40、e Loadings at Members/17X,NLM=,I4)DO 200 J=1,NM DO 100 I=1,6 F(I,J)=0.0D0100 CONTINUE200CONTINUE RETURNENDSUBROUTINE IOLMFB(NM,NJ,N,NLM,LM,LT,VF,DST,IST,IEN,X,Y,LV,F,B)! Reading and printing the data of loadings at members, calculating fixed-end ! forces and equivalent joint loads, adding the latter
41、 to the joint load vector BDIMENSION LM(NLM),LT(NLM),VF(NLM),DST(NLM),IST(NM), & IEN(NM),X(NJ),Y(NJ),LV(6,NM)DOUBLE PRECISION RL,C,S,D1,D2,P1,P2,F1,F2,F3,F4, & F5,F6,G,B(N),F(6,NM)WRITE (4,1)1FORMAT(/1X,member,2X,type,9X,VF,14X,DST)READ (3,*) (LM(M),LT(M),VF(M),DST(M),M=1,NLM)WRITE (4,2)(LM(M),LT(M)
42、,VF(M),DST(M),M=1,NLM)2FORMAT(1X,I4,I7,1X,2F16.6)DO 100 M=1,NLM L=LM(M) CALL LCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S) D1=DST(M) D2=RL-D1 IF (LT(M).EQ.1.OR.LT(M).EQ.3) THEN P1=VF(M)*C P2=-VF(M)*S ENDIF IF(LT(M).EQ.2.OR.LT(M).EQ.4) THEN P1=VF(M)*S P2=VF(M)*C ENDIF IF(LT(M).EQ.1.OR.LT(M).EQ.2)THEN F1=-P1*D2/RL F4=-P1-F1 F2=-P2*D2*D2*(RL+2.0*D1)/(RL*RL*RL) F5=-P2-F2 F3=-P2*D1*D2*D2/(RL*RL) F6=P2*D1*D1*D2/(RL*RL) ENDIF IF(LT(M).EQ.3.OR.LT(M).EQ.4)THEN G=P2*D1*D1/(12.0*RL*RL) F3=-G*(6.0*RL-8.0*D1)*RL+3.0*D1*D1) F6=G*D1*(4.0*RL-3.0*D1) F5=-6.0*G*D1*(2.0-D1/RL) F2=-P2*D1-F5 F4=-P1*D1*D1/(2.0*RL)