线弹性小变形空间板壳静力有限元计算程序Word文件下载.docx

1、 kna5=knode*1 if （kna5/2*2 .lt. kna5） kna5=kna5+1 knb4=kelem*1 if （knb4/2*2 .lt. knb4） knb4=knb4+1 knb2=kvar1*1 if （knb2/2*2 .lt. knb2） knb2=knb2+1 knb3=kvar1*1 if （knb3/2*2 .lt. knb3） knb3=knb3+1 kna0=1 kna1=kna1+kna0 kna2=kna2+kna1 kna3=kna3+kna2 kna4=kna4+kna3 kna5=kna5+kna4 if （kna5-1.gt.2500000

2、00） thenexceed memory of array iamemory of ia = 250000000memory needed = ,kna5, in prgram start stop 55555 endif knb0=1 knb1=knb1+knb0 knb2=knb2+knb1 knb3=knb3+knb2 knb4=knb4+knb3 if （knb4-1.gt.125000000） thenexceed memory of array ibmemory of ib = 125000000,knb4, call start（knode,kdgof,kcoor,kvar,

3、*kelem,maxt,kvar1,ia（kna0）,ia（kna1）,ia（kna2）, *ia（kna3）,ia（kna4）,ib（knb0）,ib（knb1）,ib（knb2）, *ib（knb3）, *filename） end subroutine start（knode,kdgof,kcoor,kvar, *kelem,maxt,kvar1,u0,u1,u2, *coor,inodvar,nodvar,numcol,lm,node, character*12 filename（20） DIMENSION NODVAR（KDGOF,KNODE）,COOR（KCOOR,KNODE）,R

4、（3）, * U0（KDGOF,KNODE）,U1（KDGOF,KNODE）,U2（KDGOF,KNODE）, * INODVAR（KNODE）,node（kelem） DIMENSION NUMCOL（KVAR1）,LM（KVAR1） CHARACTER*1 MATERIAL logical filflgC .C . KDGOF NUMBER OF D.O.FC . KNODE NUMBER OF NODESC . INODVAR ID DATAC . NODVAR DENOTE THE EQUATION NUMBER CORRESPONDING THE D.O.FC . U0 U1 U2

5、INITIAL VALUEC . COOR COORDINATESC . NODE ELEMENT NODAL CONNECTION6 FORMAT （1X, 15I4）7 FORMAT （1X,8F9.3）C.OPEN ID file OPEN （1,FILE=,FORM=UNFORMATTED,STATUS=OLD READ （1） NUMNOD,NODDOF,（NODVAR（I,J）,I=1,NODDOF）,J=1,NUMNOD） CLOSE （1） call chms（kdgof,knode,NODVAR）c WRITE（*,*） NUMNOD =,NUMNOD, NODDOF =,N

6、ODDOFc WRITE （*,*） ID =c WRITE （*,6） （NODVAR（I,J）,I=1,NODDOF）,J=1,NUMNOD）C. GET THE NATURAL NODAL ORDER DO 12 N=1,KNODE INODVAR（N）=N12 CONTINUEC. OPEN ORDER.NOD FILE AND READ THE NODAL ORDER IF THE FILE EXIST inquire（file=ORDER.NOD,exist=filflg） if （filflg） then READ （1） （INODVAR（I）,I=1,NUMNOD） CLOS

7、E（1） WRITE（*,*） NODORDER = WRITE（*,6） （INODVAR（I）,I=1,NUMNOD）C. GET NV BY ID NEQ=0 DO 20 JNOD=1,NUMNOD J=INODVAR（JNOD） DO 18 I=1,NODDOF IF （NODVAR（I,J）.NE.1） GOTO 18 NEQ = NEQ + 1 NODVAR（I,J） = NEQ18 CONTINUE20 CONTINUE DO 30 JNOD=1,NUMNOD DO 28 I=1,NODDOF IF （NODVAR（I,J）.GE.-1） GOTO 28 N = -NODVAR（

8、I,J）-1 NODVAR（I,J） = NODVAR（I,N）28 CONTINUE30 CONTINUEC. OPEN AND WRITE THE NV FILE OPEN（8,STATUS=unknown,FILE= ,FORM= WRITE（8） （NODVAR（I,J）,I=1,NODDOF）,J=1,NUMNOD） CLOSE（8）c WRITE（*,6） （NODVAR（I,J）,I=1,NODDOF）,J=1,NUMNOD）C. WRITE THE BOUNDAY CONDITION FILE BFD ACCORDING TO THE DISP0 FILEC.OPEN DISP

9、0 FILE OPEN（1,FILE= READ（1） NUMNOD,NODDOF,（U0（I,J）,I=1,NODDOF）,J=1,NUMNOD）C.OPEN BFD FILE WRITE（1） （U0（I,J）,I=1,NODDOF）,J=1,NUMNOD）C. GET THE INITIAL TIME FROM TIME0 FILEC.OPEN TIME0 FileFORMATTED READ（1,*） T0,TMAX,DT TIME = T0 IT = 0 TMAX,DT,TIME,IT =,TMAX,DT,TIME,ITC.OPEN TIME File WRITE（1） TMAX,D

10、T,TIME,ITC.OPEN COOR file READ （1） NUMNOD,NCOOR,（COOR（I,J）,I=1,NCOOR）,J=1,NUMNOD）COOR =c WRITE（*,7） （COOR（I,J）,I=1,NCOOR）,J=1,NUMNOD）C. GET THE INITIAL VALUE FROM THE DATA FILES BY PREPROCESSORdisp1 open（16,file=,status=old read（16） numnod,noddof,（U0（J,N）,J=1,NODDOF）,N=1,NUMNOD） close（16）disp2 read（

11、16） numnod,noddof,（U1（J,N）,J=1,NODDOF）,N=1,NUMNOD）disp3 read（16） numnod,noddof,（U2（J,N）,J=1,NODDOF）,N=1,NUMNOD） U0 = c WRITE（*,（6F13.3） （U0（J,N）,J=1,NODDOF）,N=1,NUMNOD）C WRITE（*,*） U1 = C WRITE（*,） （U1（J,N）,J=1,NODDOF）,N=1,NUMNOD）C. COMPUTE THE INITIAL VALUE BY BOUND.FOR zo = 0.0d0c DO 321 N=1,NUMNO

12、Dc DO 100 J=1,NCOORc100 R（J） = COOR（J,N）c DO 200 J=1,NODDOFc U0（J,N） = BOUND（R,zo,J）c U1（J,N） = BOUND1（R,zo,J）c U2（J,N） = BOUND2（R,zo,J）c200 CONTINUEc321 CONTINUEC.OPEN AND WRITE THE INITIAL VALUE FILE UNOD WRITE（1） （U0（I,J）,J=1,NUMNOD）,I=1,NODDOF）, * （U1（I,J）,J=1,NUMNOD）,I=1,NODDOF）, * （U2（I,J）,J=1

13、,NUMNOD）,I=1,NODDOF）, * （U0（I,J）,J=1,NUMNOD）,I=1,NODDOF）c. open IO file open（21,file=formatted read（21, （1a） material read（21,*） numtyp close（21） DO I=1,NEQ NUMCOL（i）=1 ENDDOC.OPEN ELEM0 file OPEN （3,FILE= NUMEL=0 KELEM=0 KEMATE=0 DO 2000 ITYP=1,NUMTYPC.INPUT ENODE READ （3） NUM,NNODE, * （NODE（I-1）*N

14、NODE+J）,J=1,NNODE）,I=1,NUM）cc WRITE（*,*） NUM =,NUM, NNODE =,NNODENODE =cc WRITE（*,6） （NODE（I-1）*NNODE+J）,J=1,NNODE）,I=1,NUM） IF （KELEM.LT.NUM*NNODE） KELEM = NUM*NNODE NNE = NNODE IF （MATERIAL.EQ.Y .OR. MATERIAL.EQ.y） THEN READ （3） MMATE,NMATE IF （KEMATE.LT.MMATE*NMATE） KEMATE = MMATE*NMATE NNE = NNE

15、-1 ENDIFMMATE =,MMATE, NMATE =,NMATE DO 1000 NE=1,NUM L=0 DO 700 INOD=1,NNE NODI=NODE（NE-1）*NNODE+INOD） DO 600 IDGF=1,KDGOF INV=NODVAR（IDGF,NODI） IF （INV.LE.0） GOTO 600 L=L+1 LM（L）=INV600 CONTINUE700 CONTINUE NUMEL=NUMEL+1C WRITE （*,*） L,LM =,LC WRITE （*,（1X,15I5） （LM（I）,I=1,L） if （l.gt.0） call ACLH

16、（NEQ,NUMCOL,l,lm）1000 continue2000 CONTINUEc CLOSE（1） CLOSE（3） call BCLH（NEQ,NUMCOL） MAXA=NUMCOL（NEQ）C.OPEN SYS File OPEN （2,FILE= WRITE（2） NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE CLOSE （2） OPEN（2,FILE= write（2） （NUMCOL（I）,I=1,NEQ） CLOSE（2）c write（*,*） NEQ,NUMCOL=,NEQc write（*,6） （NUMCOL（i）,i=1,NEQ） END

17、subroutine chms（kdgof,knode,id） dimension id（kdgof,knode）,ms（1000）,is（1000） do 1000 k=1,kdgof m = 0 do 800 n=1,knode if （id（k,n）.le.-1） id（k,n）=-1 if （id（k,n）.le.1） goto 800 j=id（k,n） j0=0 if （m.gt.0） then do i=1,m if （j.eq.ms（i） j0=is（i） enddo if （j0.eq.0） then m=m+1 ms（m）=j is（m）=n id（k,n）=1 else

18、id（k,n）=-j0-1800 continue return SUBROUTINE ACLH（NEQ,NUMCOL,ND,LM） DIMENSION LM（ND）,NUMCOL（NEQ） LS=LM（1）+1 DO 100 I=1,ND110 IF（LM（I）-LS） 120,100,100120 LS=LM（I）100 CONTINUE DO 200 I=1,ND II=LM（I） ME=II-LS IF（ME.GT.NUMCOL（II） NUMCOL（II）=ME200 CONTINUE RETURN SUBROUTINE BCLH（NEQ,NUMCOL） DIMENSION NUMC

19、OL（NEQ）C NUMCOL（1） = 1 DO 490 I=2,NEQ490 NUMCOL（I） = NUMCOL（I） + NUMCOL（I-1） + 13. eshell3da.for，Galerkin法求解位移场的主程序 common /cc/ ic（62500000） MAXT=250000000/2/2 read（2） NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE IF （MAXA.GT.MAXT） THENMATRIX A EXCEED CORE MEMERY . ,MAXAREQUIRED CORE MEMERY . ,MAXT STOP 0000 K

20、VAR=KNODE*KDGOF KCOOR=3C KELEM=31250000KNODE,KDGOF,KVAR,KCOOR,KELEM = WRITE（*,（1X,6I7） KNODE,KDGOF,KVAR,KCOOR,KELEM kna1=kdgof*knode*1 if （kna1/2*2 .lt. kna1） kna1=kna1+1 knc1=kdgof*knode*2 knc2=kcoor*knode*2 knc7=kdgof*knode*2 knc3=neq*2 knb1=maxa*2 knb2=maxa*2 kna2=neq*1 if （kna2/2*2 .lt. kna2） kna2=kna2+1 knc6=kemate*2 kna3=kelem*1 if （kna3/2*2 .lt. kna3） kna3=kna3+1 knc8=100000*2 knc5=neq*2 knc4=kdgof*knode*2 if （kna3-1.gt.125000000） then,kna3, in prgram eshell3da if （knb2-1.gt.250000000） thenexc