$LARGE C C <<< ＪＳＴＰ基礎講座用プログラム　−　メイン >>> C C @ ２次元問題（平面ひずみ／軸対称） C @ 弾性／ラグランジェ乗数法剛塑性ＦＥＭ C @ フルマトリックスソルバー C C * Ver.1.1 C Coding started in Apr. 27, 1991 C Coding ended in Aug. 27, 1991 C C * Max mesh size C Division in lateral direction Nx : 4 C thickness direction Ny : 4 C C 1.本プログラムは，ＪＳＴＰ基礎講座受講生の実習のために， C ＪＳＴＰ基礎講座委員が共同で作成したものです． C 桑原利彦（東京農工大学） C 高橋　進（日産自動車） C 湯川伸樹（名古屋大学） C 柳本　潤（東京大学） C 2.ＪＳＴＰ基礎講座受講生以外の使用を原則として禁止します． C 3.使用を希望する場合には，委員まで問い合わせて下さい． C PROGRAM JSTP IMPLICIT REAL*8(A-H,P-Z) CHARACTER*1 INCRED CHARACTER*20 FILEOUT C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLOK01/ NX,NY COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/BLKRED/ INCRED(66) COMMON/BLKPRS/ VDPRES COMMON/BLKFRC/ FRCOEF COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM C C . 出力ファイル名 WRITE(6,9950) READ(5,9951) FILEOUT OPEN(16,FILE=FILEOUT,STATUS='NEW',FORM='FORMATTED',ERR=9952) C C . 配列データの大きさ C ((( 節点数 ))) MINDX1= 25 C ((( 要素数 ))) MINDX2= 16 C C . データの入力 CALL READIN(NODNUM,NELNUM,ICAL,YMAX,REXLIM,REXFAC,NSTMAX,IERROR) IF(IERROR.GE.1) STOP C C . 収束条件，緩和係数，各ステップの圧下率 DTSTEP= YMAX/100.0/DABS(VDPRES) C C . 摩擦力のスムージングのための係数 USCRIT= DABS(VDPRES)/1000.0 C C . 繰り返し回数のカウント C ( トータル ) ITRTOT= 0 C ( ステップ ) NSTSTP= 0 C C . 各ステップの計算の開始 1 NSTSTP= NSTSTP+1 C C .. 印刷 WRITE( 6,7000) NSTSTP WRITE(16,7000) NSTSTP C C .. マトリックスの作成 CALL SVEVAL(NELNUM,ICAL) C C .. 収束計算の開始 ITRSTP= 0 2 ITRTOT= ITRTOT+1 ITRSTP= ITRSTP+1 C C ... 内部変形エネルギ，摩擦エネルギの初期化 POWERI= 0.0 POWERF= 0.0 C C ... ひずみ速度の計算 CALL STRCAL(NELNUM) C C ... 節点力の計算 IF(ICAL.LE.1) CALL BFORCE(NELNUM,ICAL,USCRIT,POWERF) C C ... 大連立方程式の組立て CALL ASMBLY(NODNUM,NELNUM,ICAL,ITRSTP,POWERI) C C ... 大連立方程式の求解，節点速度の修正 CALL SOLVER(NODNUM,NELNUM,REXFAC,REXNRM) C C ... 収束の判定 IF(MOD(ITRSTP-2,10).EQ.0) THEN WRITE( 6,7001) NSTSTP,ITRTOT,ITRSTP,REXNRM WRITE(16,7001) NSTSTP,ITRTOT,ITRSTP,REXNRM ELSE WRITE( 6,7002) NSTSTP,ITRTOT,ITRSTP,REXNRM WRITE(16,7002) NSTSTP,ITRTOT,ITRSTP,REXNRM END IF IF(REXNRM.LE.REXLIM) GO TO 3 IF(ICAL .GE.2 ) GO TO 3 C C ... 次の収束計算の開始 GO TO 2 C C .. 結果の印刷，節点座標の更新 3 CALL PRTOUT(NODNUM,NELNUM,ICAL,POWERI,POWERF,DTSTEP) C C .. 解析の終了 IF(NSTSTP.LT.NSTMAX) GO TO 1 C C . メッセージの出力 WRITE( 6,7003) NSTSTP,ITRTOT WRITE(16,7003) NSTSTP,ITRTOT C C . フォーマット文 9950 FORMAT(//' ... 出力ファイル名を入力して下さい') 9951 FORMAT(A20) 7000 FORMAT(//' *** ',I3,' ステップの計算を開始します ***') 7001 FORMAT( ' ... NSTSTP=',I3,' ITRTOT=',I3, * ' ITRSTP=',I3,' REXNRM=',D12.4) 7002 FORMAT(1H+,' ... NSTSTP=',I3,' ITRTOT=',I3, * ' ITRSTP=',I3,' REXNRM=',D12.4) 7003 FORMAT(//' *** 計算終了 ***', * ' ステップ数=',I3,' 繰り返し数=',I3) 9952 STOP END C C SUB. READIN SUBROUTINE READIN(NODNUM,NELNUM,ICAL,YMAX, * REXLIM,REXFAC,NSTMAX,IERROR) IMPLICIT REAL*8 (A-H,P-Z) CHARACTER*1 INCRED C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLOK01/ NX,NY COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/BLKRED/ INCRED(66) COMMON/BLKPRS/ VDPRES COMMON/BLKFRC/ FRCOEF COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM COMMON/SRBLK2/ EQHIST(16) C C ... 加工条件の読み込み WRITE(6, 999) READ(5,*) ICAL WRITE(6,1000) READ(5,*) XMAX,XMIN,YMAX WRITE(6,998) READ(5,*) PLXMAX WRITE(6,1001) READ(5,*) NX,NY IF(ICAL.LE.1) THEN WRITE(6,1002) READ(5,*) S0,SA,SB,PH,PR WRITE(6,1003) READ(5,*) FRCOEF WRITE(6,1004) READ(5,*) VDPRES WRITE(6,1005) READ(5,*) REXLIM,REXFAC,NSTMAX ELSE WRITE(6,1011) READ(5,*) YNGMOD,POISNM WRITE(6,1012) READ(5,*) VDPRES FRCOEF= 0.0 REXLIM= 0.0 REXFAC= 1.0 NSTMAX= 0 END IF C C ... 配列サイズのチェック IERROR= 0 LARAY1=(NX+1)*(NY+1) LARAY2= NX * NY IF(LARAY1.GT.MINDX1) THEN WRITE( 6,9000) LARAY1 WRITE(16,9000) LARAY1 IERROR= 1 END IF IF(LARAY2.GT.MINDX2) THEN WRITE( 6,9001) LARAY2 WRITE(16,9001) LARAY2 IERROR= 1 END IF IF(IERROR.GE.1) RETURN C C ... 初期要素分割の作成 NODNUM= 0 DO 100 IY=1,NY+1 YCOTMP= YMAX *DFLOAT(IY-1)/DFLOAT(NY) DO 100 IX=1,NX+1 XCOTMP= (XMAX-XMIN)*DFLOAT(IX-1)/DFLOAT(NX)+XMIN NODNUM= NODNUM+1 X(NODNUM)= XCOTMP Y(NODNUM)= YCOTMP 100 CONTINUE C C ... 境界条件の設定 NODNUM= 0 DO 102 IY=1,NY+1 DO 102 IX=1,NX+1 NODNUM= NODNUM+1 IF(IY.EQ.1) THEN IF(IX.EQ.1) THEN INCRED(2*NODNUM-1)= 'C' INCRED(2*NODNUM )= 'C' ELSE INCRED(2*NODNUM-1)= ' ' INCRED(2*NODNUM )= 'C' END IF END IF IF(IY.GT.1.AND.IY.LT.NY+1) THEN IF(IX.EQ.1) THEN INCRED(2*NODNUM-1)= 'C' INCRED(2*NODNUM )= ' ' ELSE INCRED(2*NODNUM-1)= ' ' INCRED(2*NODNUM )= ' ' END IF END IF IF(IY.EQ.NY+1) THEN IF(IX.EQ.1) THEN INCRED(2*NODNUM-1)= 'C' INCRED(2*NODNUM )= 'P' ELSE INCRED(2*NODNUM-1)= ' ' INCRED(2*NODNUM )= 'P' END IF IF(X(NODNUM).GT.PLXMAX) INCRED(2*NODNUM )= ' ' END IF 102 CONTINUE IF(XMIN.GT.0.0) THEN NODNUM= 0 DO 103 IY=1,NY+1 DO 103 IX=1,NX+1 NODNUM= NODNUM+1 IF(IX.EQ.1) INCRED(2*NODNUM-1)= ' ' 103 CONTINUE END IF C C ... 初期速度分布の作成 NODNUM= 0 DO 101 IY=1,NY+1 YVETMP= VDPRES*DFLOAT(IY-1)/DFLOAT(NY) DO 101 IX=1,NX+1 XVETMP= -VDPRES*0.1 NODNUM= NODNUM+1 IF(ICAL.LE.1) THEN VX(NODNUM)= XVETMP VY(NODNUM)= YVETMP ELSE VX(NODNUM)= 0.0 VY(NODNUM)= 0.0 IF(INCRED(2*NODNUM).EQ.'P') VY(NODNUM)= VDPRES END IF 101 CONTINUE C C ... インシデンスの計算 NELNUM= 0 DO 200 IY=1,NY DO 200 IX=1,NX NELNUM= NELNUM+1 NLTMP1= (IY-1)*(NX+1)+IX NLTMP2= (IY-1)*(NX+1)+IX+1 NLTMP3= IY *(NX+1)+IX+1 NLTMP4= IY *(NX+1)+IX INCDNS(NELNUM,1)= NLTMP1 INCDNS(NELNUM,2)= NLTMP2 INCDNS(NELNUM,3)= NLTMP3 INCDNS(NELNUM,4)= NLTMP4 200 CONTINUE DO 201 NE=1,NELNUM EQHIST( NE)= 0.0 IF(ICAL.LE.1) THEN INCRED(2*NODNUM+NE)= ' ' ELSE INCRED(2*NODNUM+NE)= '*' END IF 201 CONTINUE C C ... フォーマット文 999 FORMAT(//' *** ＪＴＳＰ基礎講座プログラムを開始します ***', * //' なお，受講生以外の使用は禁止されています', * ///' ... 解析モードを指定して下さい', * /' 0 - 剛塑性変形，平面ひずみ', * /' 1 - 剛塑性変形，軸対称', * /' 2 - 弾性変形，平面ひずみ', * /' 3 - 弾性変形，軸対称') 1000 FORMAT(//' ... ビレットの 1/2外径， 1/2内径， 1/2初期高さを', * '入力して下さい') 998 FORMAT(//' ... 工具の幅を入力して下さい') 1001 FORMAT(//' ... 半径方向，高さ方向分割数を入力して下さい') 1002 FORMAT(//' ... 材料特性値 S0,a,b,n,m を入力して下さい') 1003 FORMAT(//' ... 摩擦係数を入力して下さい（＜０：クーロン摩擦）') 1004 FORMAT(//' ... 工具の移動速度を入力して下さい') 1005 FORMAT(//' ... 収束判定値，緩和係数，最大ステップ数を', * '入力して下さい') 1011 FORMAT(//' ... ヤング率，ポアソン比を入力して下さい') 1012 FORMAT(//' ... 工具の変位を入力して下さい') 9000 FORMAT(//' @@@ 節点数が大き過ぎます @@@ 節点数＝',I3) 9001 FORMAT(//' @@@ 要素数が大き過ぎます @@@ 要素数＝',I3) C RETURN END C C SUBR. SVEVAL SUBROUTINE SVEVAL(NELNUM,ICAL) IMPLICIT REAL*8(A-H,P-Z) DIMENSION GZAI(4),EATA(4),XE(4),YE(4),SS(4,4), * SW(4,4),SL(4),FAIX(4),FAIY(4),A(2,2), * SHPFN0(4),SHPFN1(4),SHPFN2(4),SHPF12(4), * BMGZAI(4,8),BMEATA(4,8),WKARYC(8,4),BMATMP(4,8) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/MATRX0/ BMATRX(16,4,8) COMMON/MATRX1/ CMATRX(16,8,8) COMMON/MATRX2/ PMATRX(16, 8) COMMON/MATRX3/ DETVOL(16) COMMON/MATRX4/ DM(4,4) C C . サブマトリックスの初期化 DO 1 I=1,4 DO 1 J=1,4 1 SS(I,J)= 0.0 DO 2 I=1,2 2 SS(I,I)= 1.0 SS(3,3)= 2.0 IF(ICAL.EQ.1.OR.ICAL.EQ.3) SS(4,4)= 1.0 C C . ［Ｄ］マトリックスの計算 C .. 剛塑性解析 DO 3 I=1,4 DO 3 J=1,4 DM(I,J)= 0.0 3 SL(I )= 0.0 IF(ICAL.LE.1) THEN DM(1,1)= 1.0 DM(2,2)= 1.0 DM(3,3)= 1.0 SL(1 )= 1.0 SL(2 )= 1.0 SL(3 )= 0.0 IF(ICAL.EQ.1) THEN DM(4,4)= 1.0 SL(4 )= 1.0 END IF C .. 弾性解析 ELSE CFTEMP= YNGMOD*(1.0-POISNM)/(1.0+POISNM)/(1.0-2.0*POISNM) DO 4 I=1,2 4 DM(I,I)= CFTEMP* 1.0 DM(1,2)= CFTEMP* POISNM /(1.0-POISNM) DM(2,1)= CFTEMP* POISNM /(1.0-POISNM) DM(3,3)= CFTEMP*(1.0-2.0*POISNM)/(1.0-POISNM) IF(ICAL.EQ.3) THEN DM(4,4)= CFTEMP*1.0 DM(1,4)= CFTEMP*POISNM/(1.0-POISNM) DM(4,1)= CFTEMP*POISNM/(1.0-POISNM) DM(2,4)= CFTEMP*POISNM/(1.0-POISNM) DM(4,2)= CFTEMP*POISNM/(1.0-POISNM) END IF SL(1)= 0.0 SL(2)= 0.0 SL(3)= 0.0 SL(4)= 0.0 END IF C .. [SW]=[SS][DM] DO 5 I=1,4 DO 5 J=1,4 SUMTMP= 0.0 DO 6 K=1,4 6 SUMTMP= SUMTMP+ SS(I,K)*DM(K,J) 5 SW(I,J)= SUMTMP C C . ［Ｂ］マトリックスの計算 GZAI(1)= -1.0 EATA(1)= -1.0 GZAI(2)= 1.0 EATA(2)= -1.0 GZAI(3)= 1.0 EATA(3)= 1.0 GZAI(4)= -1.0 EATA(4)= 1.0 XINTPN = 0.0 YINTPN = 0.0 DO 1000 NE=1,NELNUM DO 20 JNUM=1,4 KE= INCDNS(NE,JNUM) XE(JNUM)= X(KE) 20 YE(JNUM)= Y(KE) CALL BMACAL(ICAL,XE,YE,XINTPN,YINTPN,GZAI,EATA, * SHPFN0,SHPFN1,SHPFN2,BMATMP,DETTMP) DO 21 INUM=1,4 DO 21 JNUM=1,8 21 BMATRX(NE,INUM,JNUM)= BMATMP(INUM,JNUM) DETVOL(NE )= DETTMP C C . 安定化［Ｂ］マトリックスの計算（４節点双１次要素） C .. 形状関数の微係数 DO 23 JNUM=1,4 23 SHPF12(JNUM)= GZAI(JNUM)*EATA(JNUM)/4.0 SX0 = 0.0 SX1 = 0.0 SX2 = 0.0 SY1 = 0.0 SY2 = 0.0 SX12= 0.0 SY12= 0.0 DO 24 JNUM=1,4 XX= XE(JNUM) YY= YE(JNUM) SX0 = SX0 + SHPFN0(JNUM)*XX SX1 = SX1 + SHPFN1(JNUM)*XX SX2 = SX2 + SHPFN2(JNUM)*XX SY1 = SY1 + SHPFN1(JNUM)*YY SY2 = SY2 + SHPFN2(JNUM)*YY SX12= SX12+ SHPF12(JNUM)*XX SY12= SY12+ SHPF12(JNUM)*YY 24 CONTINUE A(1,1)= SX1 A(1,2)= SY1 A(2,1)= SX2 A(2,2)= SY2 CALL INVERS(A,RDUMMY) DO 22 JNUM=1,4 DS1= A(1,1)*SHPFN1(JNUM)+ A(1,2)*SHPFN2(JNUM) DS2= A(2,1)*SHPFN1(JNUM)+ A(2,2)*SHPFN2(JNUM) FAIX(JNUM)= DS1 22 FAIY(JNUM)= DS2 C C .. 安定化［Ｂ］マトリックスの成分 DO 26 JNUM=1,4 ESG1= A(1,1)* 0.0 * + A(1,2)*(SHPF12(JNUM)-SX12*FAIX(JNUM)-SY12*FAIY(JNUM)) ESG2= A(2,1)* 0.0 * + A(2,2)*(SHPF12(JNUM)-SX12*FAIX(JNUM)-SY12*FAIY(JNUM)) ESE1= A(1,1)*(SHPF12(JNUM)-SX12*FAIX(JNUM)-SY12*FAIY(JNUM)) * + A(1,2)* 0.0 ESE2= A(2,1)*(SHPF12(JNUM)-SX12*FAIX(JNUM)-SY12*FAIY(JNUM)) * + A(2,2)* 0.0 BMGZAI(1,2*JNUM-1)= ESG1 BMGZAI(1,2*JNUM )= 0.0 BMGZAI(2,2*JNUM-1)= 0.0 BMGZAI(2,2*JNUM )= ESG2 BMGZAI(3,2*JNUM-1)= ESG2/ 2.0 BMGZAI(3,2*JNUM )= ESG1/ 2.0 BMGZAI(4,2*JNUM-1)= 0.0 BMGZAI(4,2*JNUM )= 0.0 BMEATA(1,2*JNUM-1)= ESE1 BMEATA(1,2*JNUM )= 0.0 BMEATA(2,2*JNUM-1)= 0.0 BMEATA(2,2*JNUM )= ESE2 BMEATA(3,2*JNUM-1)= ESE2/ 2.0 BMEATA(3,2*JNUM )= ESE1/ 2.0 BMEATA(4,2*JNUM-1)= 0.0 BMEATA(4,2*JNUM )= 0.0 IF(ICAL.EQ.1) THEN BMGZAI(4,2*JNUM-1)=(SHPFN1(JNUM)*SX0-SHPFN0(JNUM)*SX1)/SX0**2 BMEATA(4,2*JNUM-1)=(SHPFN2(JNUM)*SX0-SHPFN0(JNUM)*SX2)/SX0**2 END IF 26 CONTINUE C C . ［Ｃ］，［Ｐ］マトリックスの計算 C .. [C]=[B]t[SW][B] DO 30 I=1,8 DO 30 J=1,8 30 CMATRX(NE,I,J)= 0.0 DO 31 I=1,8 DO 31 J=1,4 SMTMPC= 0.0 SMTMPE= 0.0 DO 32 K=1,4 32 SMTMPC= SMTMPC+ BMATRX(NE,K,I)*SW(K,J) 31 WKARYC(I,J)= SMTMPC DO 33 I=1,8 DO 33 J=1,8 SMTMPC= 0.0 DO 34 K=1,4 34 SMTMPC= SMTMPC+ WKARYC(I,K)*BMATRX(NE,K,J) 33 CMATRX(NE,I,J)= SMTMPC C .. [P]=[SL][B] DO 35 J=1,8 SUMTMP= 0.0 DO 36 K=1,4 36 SUMTMP= SUMTMP+ SL(K)*BMATRX(NE,K,J) 35 PMATRX(NE,J)= SUMTMP C ... [C]=[C]+[B]gzai-t[SW][B]gzai DO 41 I=1,8 DO 41 J=1,4 SMTMPC= 0.0 SMTMPE= 0.0 DO 42 K=1,4 42 SMTMPC= SMTMPC+ BMGZAI(K,I)*SW(K,J) 41 WKARYC(I,J)= SMTMPC DO 43 I=1,8 DO 43 J=1,8 SMTMPC= 0.0 DO 44 K=1,4 44 SMTMPC= SMTMPC+ WKARYC(I,K)*BMGZAI(K,J) 43 CMATRX(NE,I,J)= CMATRX(NE,I,J)+SMTMPC/3.0 C ... [C]=[C]+[B]eata-t[SW][B]eata DO 46 I=1,8 DO 46 J=1,4 SMTMPC= 0.0 SMTMPE= 0.0 DO 47 K=1,4 47 SMTMPC= SMTMPC+ BMEATA(K,I)*SW(K,J) 46 WKARYC(I,J)= SMTMPC DO 48 I=1,8 DO 48 J=1,8 SMTMPC= 0.0 DO 49 K=1,4 49 SMTMPC= SMTMPC+ WKARYC(I,K)*BMEATA(K,J) 48 CMATRX(NE,I,J)= CMATRX(NE,I,J)+SMTMPC/3.0 1000 CONTINUE C RETURN END C C SUBR. STRCAL SUBROUTINE STRCAL(NELNUM) IMPLICIT REAL*8(A-H,P-Z) DIMENSION VENODE(8) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/MATRX0/ BMATRX(16,4,8) COMMON/MATRX1/ CMATRX(16,8,8) COMMON/MATRX2/ PMATRX(16, 8) COMMON/MATRX3/ DETVOL(16) COMMON/SRBLK1/ EQRATE(16),STRATE(16,4) COMMON/SRBLK2/ EQHIST(16) C C . 相当ひずみ速度の計算（剛塑性解析） DO 1000 NE=1,NELNUM EXXDOT= 0.0 EYYDOT= 0.0 EXYDOT= 0.0 ESSDOT= 0.0 DO 10 JNUM=1,4 KENODE= INCDNS(NE,JNUM) VENODE(2*JNUM-1)= VX(KENODE) 10 VENODE(2*JNUM )= VY(KENODE) DO 20 LNUM=1,8 EXXDOT= EXXDOT+ BMATRX(NE,1,LNUM)*VENODE(LNUM) EYYDOT= EYYDOT+ BMATRX(NE,2,LNUM)*VENODE(LNUM) EXYDOT= EXYDOT+ BMATRX(NE,3,LNUM)*VENODE(LNUM) 20 ESSDOT= ESSDOT+ BMATRX(NE,4,LNUM)*VENODE(LNUM) STRATE(NE,1)= EXXDOT STRATE(NE,2)= EYYDOT STRATE(NE,3)= EXYDOT STRATE(NE,4)= ESSDOT EQTEMP= EXXDOT**2+EYYDOT**2+2.0*EXYDOT**2+ESSDOT**2 1000 EQRATE(NE)= DSQRT(2.0D+00*EQTEMP/3.0D+00) C RETURN END C C SUB. BFORCE SUBROUTINE BFORCE(NELNUM,ICAL,USCRIT,POWERF) IMPLICIT REAL*8(A-H,P-Z) CHARACTER*1 INCRED DIMENSION XINT(5),GZAI(2),FAISHP(5,2),WEIGHT(5),NF(2), * FAIWRK(2,4),FORSTF(4,4) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLOK01/ NX,NY COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/BLKRED/ INCRED(66) COMMON/BLKFRC/ FRCOEF COMMON/MATRX8/ DMRIGH(66) COMMON/MATRX9/ DMLEFT(66,66) COMMON/BLOPRS/ PRESUR(16) COMMON/SRBLK1/ EQRATE(16),STRATE(16,4) COMMON/SRBLK2/ EQHIST(16) C C . 初期化 NMINTE= 5 PAI= 3.1415926535897932D+00 DO 1 I=1,2*MINDX1+MINDX2 DMRIGH(I )= 0.0 DO 2 J=1,2*MINDX1+MINDX2 2 DMLEFT(I,J)= 0.0 1 CONTINUE IF(ICAL.GE.2) RETURN C C . 形状関数 GZAI(1)= -1.0 GZAI(2)= 1.0 XINT(1)= -0.906179845938664D+00 XINT(2)= -0.538469310105683D+00 XINT(3)= 0.0 XINT(4)= 0.538469310105683D+00 XINT(5)= 0.906179845938664D+00 WEIGHT(1)= 0.236926885056189D+00 WEIGHT(2)= 0.478628670499366D+00 WEIGHT(3)= 0.568888888888889D+00 WEIGHT(4)= 0.478628670499366D+00 WEIGHT(5)= 0.236926885056189D+00 DO 3 INTE=1,NMINTE DO 3 JNUM=1,2 3 FAISHP(INTE,JNUM)= (1+XINT(INTE)*GZAI(JNUM))/2.0 C C . 節点力の計算 DO 1000 NE=(NY-1)*NX+1,NELNUM NF(1)= INCDNS(NE,4) NF(2)= INCDNS(NE,3) DO 4 I=1,4 DO 4 J=1,4 4 FORSTF(I,J)= 0.0 IF(INCRED(2*NF(1)).EQ.'P'.AND.INCRED(2*NF(2)).EQ.'P') THEN C .. 線素の長さ DL= X(NF(2))-X(NF(1)) IF(ICAL.EQ.1) DL= DL*PAI*(X(NF(1))+X(NF(2))) C .. 接触圧力のチェック ER= EQRATE(NE) EH= EQHIST(NE) SM= YESTRS(EH,ER) ST= YESTRS(EH,ER)/DSQRT(3.0D+00) CF= 2.0*SM/3.0/ER SYY= CF*STRATE(NE,2)+PRESUR(NE) IF(FRCOEF.GE.0.0) THEN FRS= FRCOEF*ST ELSE FRS= FRCOEF*SYY IF(DABS(FRS).GE.ST) FRS= ST END IF IF(FRS.LE.0.0) FRS= 0.0 C .. 節点力の計算 DO 100 INTE=1,NMINTE USX= 0.0 DO 21 JNUM=1,2 21 USX= USX+ FAISHP(INTE,JNUM)*VX(NF(JNUM)) USN= DABS(USX) FRC= -2.0/PAI*DATAN(USN/USCRIT) FRX= USX/USN*FRC*FRS DO 22 JNUM=1,2 LN= NF(JNUM) DMRIGH(2*LN-1)= DMRIGH(2*LN-1) * + FRX*FAISHP(INTE,JNUM)*WEIGHT(INTE)*DL/2.0 22 CONTINUE C .. 摩擦エネルギーの計算 POWERF= POWERF- FRX*USX*WEIGHT(INTE)*DL/2.0 C .. 係数マトリックスの計算 DO 31 IWRK=1,2 DO 31 JWRK=1,4 31 FAIWRK(IWRK,JWRK)= 0.0 FAIWRK(1,1)= FAISHP(INTE,1) FAIWRK(2,2)= FAISHP(INTE,1) FAIWRK(1,3)= FAISHP(INTE,2) FAIWRK(2,4)= FAISHP(INTE,2) FRM= 2.0/PAI*DATAN(USN/USCRIT)*FRS/USN DO 32 ITMP=1,4 DO 32 JTMP=1,4 SUMTMP= 0.0 DO 33 KTMP=1,2 33 SUMTMP= SUMTMP+ FAIWRK(KTMP,ITMP)*FAIWRK(KTMP,JTMP) 32 FORSTF(ITMP,JTMP)= FORSTF(ITMP,JTMP) * + FRM*SUMTMP*WEIGHT(INTE)*DL/2.0 100 CONTINUE C .. アセンブリー DO 41 INUM=1,2 KI= NF(INUM) I1= 2*INUM-1 I2= 2*INUM DO 41 JNUM=1,2 KJ= NF(JNUM) J1= 2*JNUM-1 J2= 2*JNUM DMLEFT(2*KI-1,2*KJ-1)= DMLEFT(2*KI-1,2*KJ-1)+ FORSTF(I1,J1) DMLEFT(2*KI-1,2*KJ )= DMLEFT(2*KI-1,2*KJ )+ FORSTF(I1,J2) DMLEFT(2*KI ,2*KJ-1)= DMLEFT(2*KI ,2*KJ-1)+ FORSTF(I2,J1) DMLEFT(2*KI ,2*KJ )= DMLEFT(2*KI ,2*KJ )+ FORSTF(I2,J2) 41 CONTINUE END IF 1000 CONTINUE C RETURN END C C SUB. ASMBLY SUBROUTINE ASMBLY(NODNUM,NELNUM,ICAL,ITRSTP,POWERI) IMPLICIT REAL*8(A-H,P-Z) CHARACTER*1 MAP(66) DIMENSION RPSTF1(8,8),RPSTF2(8,8),VENODE(8),WK(8) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/MATRX0/ BMATRX(16,4,8) COMMON/MATRX1/ CMATRX(16,8,8) COMMON/MATRX2/ PMATRX(16, 8) COMMON/MATRX3/ DETVOL(16) COMMON/MATRX8/ DMRIGH(66) COMMON/MATRX9/ DMLEFT(66,66) COMMON/SRBLK1/ EQRATE(16),STRATE(16,4) COMMON/SRBLK2/ EQHIST(16) C C . 剛性方程式両辺の計算 PAI= 3.14159265D+00 DO 1000 NE=1,NELNUM SX0= 0.0 DO 11 JNUM=1,4 KENODE= INCDNS(NE,JNUM) SX0= SX0+ X(KENODE)/4.0 VENODE(2*JNUM-1)= VX(KENODE) 11 VENODE(2*JNUM )= VY(KENODE) VOLTMP= 4.0*DETVOL(NE) IF(ICAL.EQ.1.OR.ICAL.EQ.3) VOLTMP= VOLTMP*2.0*PAI*SX0 IF(ICAL.LE.1) THEN ER= EQRATE(NE) EH= EQHIST(NE) SM= YESTRS(EH,ER) SD= YEDERV(EH,ER) C1= SM *2.0D+00/3.0D+00/ER C2= (ER*SD-SM)*4.0D+00/9.0D+00/ER**3 C1= C1* VOLTMP C2= C2* VOLTMP ELSE C1= VOLTMP C2= 0.0 END IF DO 21 I=1,8 WK(I)= 0.0 DO 21 J=1,8 21 WK(I)= WK(I)+ CMATRX(NE,I,J)*VENODE(J) DO 22 I=1,8 DO 22 J=1,8 RPSTF1(I,J)= C1*CMATRX(NE,I,J)+ C2*WK(I)*WK(J) 22 RPSTF2(I,J)= C1*CMATRX(NE,I,J) SUMTMP= 0.0 DO 23 I=1,8 23 SUMTMP= SUMTMP+ WK(I)*VENODE(I) POWERI= POWERI+ C1*SUMTMP C ... 右辺 DO 51 JNUM=1,4 KENODE= INCDNS(NE,JNUM) VR1= 0.0 VR2= 0.0 DO 52 J=1,8 VR1= VR1+ RPSTF2(2*JNUM-1,J)*VENODE(J) 52 VR2= VR2+ RPSTF2(2*JNUM ,J)*VENODE(J) IN1= 2*KENODE-1 IN2= 2*KENODE DMRIGH(IN1)= DMRIGH(IN1)- VR1 51 DMRIGH(IN2)= DMRIGH(IN2)- VR2 VR3= 0.0 DO 53 J=1,8 53 VR3= VR3+ PMATRX(NE,J)*VENODE(J) DMRIGH(2*NODNUM+NE)= DMRIGH(2*NODNUM+NE)- VR3*VOLTMP C ... 左辺のアセンブリー DO 61 I=1,4 KI = INCDNS(NE,I) I1 = 2* I-1 I2 = 2* I LI1= 2*KI-1 LI2= 2*KI DO 61 J=1,4 KJ = INCDNS(NE,J) J1 = 2* J-1 J2 = 2* J LJ1= 2*KJ-1 LJ2= 2*KJ DMLEFT(LI1,LJ1)= DMLEFT(LI1,LJ1)+ RPSTF1(I1,J1) DMLEFT(LI1,LJ2)= DMLEFT(LI1,LJ2)+ RPSTF1(I1,J2) DMLEFT(LI2,LJ1)= DMLEFT(LI2,LJ1)+ RPSTF1(I2,J1) DMLEFT(LI2,LJ2)= DMLEFT(LI2,LJ2)+ RPSTF1(I2,J2) 61 CONTINUE DO 71 J=1,4 KJ = INCDNS(NE,J) J1 = 2* J-1 J2 = 2* J LJ1= 2*KJ-1 LJ2= 2*KJ DMLEFT(2*NODNUM+NE,LJ1) * = DMLEFT(2*NODNUM+NE,LJ1)+ VOLTMP*PMATRX(NE,J1) DMLEFT(2*NODNUM+NE,LJ2) * = DMLEFT(2*NODNUM+NE,LJ2)+ VOLTMP*PMATRX(NE,J2) DMLEFT(LJ1,2*NODNUM+NE) * = DMLEFT(LJ1,2*NODNUM+NE)+ VOLTMP*PMATRX(NE,J1) DMLEFT(LJ2,2*NODNUM+NE) * = DMLEFT(LJ2,2*NODNUM+NE)+ VOLTMP*PMATRX(NE,J2) 71 CONTINUE 1000 CONTINUE C C . 剛性マトリックスの内容の確認 IF(ITRSTP.EQ.1) THEN WRITE( 6,*) WRITE( 6,*) WRITE( 6,3000) WRITE(16,*) WRITE(16,*) WRITE(16,3000) DO 2000 I=1,2*NODNUM+NELNUM DO 2001 J=1,2*NODNUM+NELNUM IF(DMLEFT(I,J).EQ.0.0) THEN MAP(J)= '.' ELSE MAP(J)= '*' END IF 2001 CONTINUE WRITE( 6,3001) I,(MAP(J),J=1,2*NODNUM+NELNUM) WRITE(16,3001) I,(MAP(J),J=1,2*NODNUM+NELNUM) 2000 CONTINUE WRITE( 6,*) WRITE( 6,*) WRITE(16,*) WRITE(16,*) END IF C C . フォーマット文 3000 FORMAT( ' @@@ 剛性マトリックスの内容 @@@') 3001 FORMAT(I5,3X,66A1) RETURN END C C FUNC. YESTRS DOUBLE PRECISION FUNCTION YESTRS(EH,ER) IMPLICIT REAL*8(A-H,P-Z) C COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM C YESTRS= S0*(SA+SB*EH)**PH IF(PR.GT.0.0) THEN YESTRS= YESTRS*ER**PR ELSE IF(ER.LT.ERMINI) YESTRS= YESTRS*(ER/ERMINI)**ERRATE END IF RETURN END C C FUNC. YEDERV DOUBLE PRECISION FUNCTION YEDERV(EH,ER) IMPLICIT REAL*8(A-H,P-Z) C COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM C YEDERV= 0.0 IF(PR.GT.0.0) THEN V1= S0*(SA+SB*EH)**PH V2= PR*ER**(PR-1.0) YEDERV= V1*V2 END IF RETURN END C C SUBR. SOLVER SUBROUTINE SOLVER(NODNUM,NELNUM,REXFAC,REXNRM) IMPLICIT REAL*8(A-H,P-Z) CHARACTER*1 INCRED DIMENSION VELSOL(50),SCALEF(66),NNWORK(66) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLOPRS/ PRESUR(16) COMMON/BLKRED/ INCRED(66) COMMON/BLKPRS/ VDPRES COMMON/MATRX8/ DMRIGH(66) COMMON/MATRX9/ DMLEFT(66,66) C C . 初期化 NODMAX= 2*NODNUM+NELNUM C C . スケーリング RSCALE= 0.0 DO 11 I=1,NODMAX VMXTMP= 0.0 DO 12 J=1,NODMAX VCURNT= DABS(DMLEFT(I,J)) IF(VCURNT.GE.VMXTMP) VMXTMP= VCURNT 12 CONTINUE SCALEF(I)= VMXTMP IF(VMXTMP.EQ.0.0) GO TO 11 VCURNT= DABS(DMRIGH(I))/DSQRT(VMXTMP) IF(VCURNT.GE.RSCALE) RSCALE= VCURNT 11 CONTINUE L1= 0 DO 21 I=1,NODMAX IF(INCRED(I).NE.' ') GO TO 21 L1= L1+ 1 L2= 0 DO 22 J=1,NODMAX IF(INCRED(J).NE.' ') GO TO 22 L2= L2+ 1 DMLEFT(L1,L2)= DMLEFT(I,J)/DSQRT(SCALEF(I))/DSQRT(SCALEF(J)) 22 CONTINUE IF(RSCALE.EQ.0.0) GO TO 21 DMRIGH(L1 )= DMRIGH(I )/DSQRT(SCALEF(I))/RSCALE 21 CONTINUE C C . 連立一次方程式の求解 NDEGRE= L1 NARRAY= 2*MINDX1+MINDX2 EPS= 1.0D-14*DFLOAT(L1) CALL HGAUSS(DMLEFT,NDEGRE,NARRAY,DMRIGH,NNWORK,IERROR,EPS) C C . 解の代入 L1= 0 DO 31 I=1,2*NODNUM VELSOL(I)= 0.0 IF(INCRED(I).NE.' ') GO TO 31 L1= L1+ 1 VELSOL(I)= DMRIGH(L1)/DSQRT(SCALEF( I))*RSCALE 31 CONTINUE DO 32 I=1,NELNUM IF(INCRED(2*NODNUM+I).NE.' ') GO TO 32 L1= L1+ 1 PRESUR(I)= DMRIGH(L1)/DSQRT(SCALEF(2*NODNUM+I))*RSCALE 32 CONTINUE REXREF= 0.0 REXNRM= 0.0 DO 33 I=1,NODNUM VX(I)= VX(I)+REXFAC*VELSOL(2*I-1) VY(I)= VY(I)+REXFAC*VELSOL(2*I ) IF(INCRED(2*I-1).EQ.'C') VX(I)= 0.0 IF(INCRED(2*I ).EQ.'C') VY(I)= 0.0 IF(INCRED(2*I ).EQ.'P') VY(I)= VDPRES REXREF= REXREF+ DSQRT(VX ( I )**2+VY ( I )**2) 33 REXNRM= REXNRM+ DSQRT(VELSOL(2*I-1)**2+VELSOL(2*I )**2) REXNRM= REXNRM/ REXREF C RETURN END SUBROUTINE HGAUSS( A,N,ND,AN,NANS,NSTOP,EPS ) C C SUBPROGRAM FOR LINEAR EQUATION PROBLEM C (A) X = B C FOR REAL GENERAL MATRIX C C * USAGE C CALL HGAUSS( A,N,ND,AN,NANS,NSTOP,EPS ) C C INPUT C A R*8 2-DIM. ARRAY (N1,N) C CONTAINING REAL GENERAL MATRIX. C N I*4 ORDER OF MATRIX C ND I*4 ROW SIZE OF THE 2-DIM. ARRAY( A ). C AN R*8 1-DIM. ARRAY ( N ) C CONTAINING RIGHT HAND VECTOR. C EPS R*8 PARAMETER TO CHECK SINGULARITY OF THE MATRIX. C STANDARD VALUE = 1.0D-14 C C OUTPUT C AN R*8 1-DIM. ARRAY ( N ) C SOLUTION VECTOR OF LINEAR EQUATION. C NSTOP I*4 ERROR CODE C =1, FOR NORMAL EXECUTION C =2, FOR SINGULAR MATRIX C =3, FOR INVALID ARGUMENT C C WORKING AREA C NANS I*4 1-DIM. ARRAY ( N ) C C C * NOTE C THE MATRIX (A) AND VECTOR B ARE DESTROYED. C C * METHOD C GAUSS ELIMINATION C C C IMPLICIT REAL*8 ( A-H,O-Z ) DIMENSION A(ND,N),AN(N),NANS(N) DATA ZERO/0.0D0/ , ONE/1.0D0/ , CONST/0.01D0/ C C ARGUMENTS CHECK C IF ( N ) 600,600,100 100 IF ( N-ND ) 110,110,600 110 IF ( EPS ) 120,130,130 120 EPS = 1.0D-14 130 CONTINUE NSTOP = 1 NM = N-1 EPSS = EPS*CONST C C SEARCH FOR LARGEST ELEMENT C DO 350 K=1,N AMAX = ZERO DO 200 I=K,N AX = A(I,K) IF( DABS(AMAX) .GE. DABS(AX) ) GO TO 200 IROW = I AMAX = AX 200 CONTINUE IF( DABS(AMAX) .LE. EPS ) GO TO 500 IF( K .EQ. N ) GO TO 350 C C GAUSS ELIMINATION STARTS C NANS(K) = IROW IF( NANS( K ) .EQ. K ) GO TO 300 A(IROW,K) = A(K,K) A(K,K) = AMAX 300 CONTINUE PIVOT = -ONE/AMAX K1 = K+1 DO 320 J=K1,N T = A(IROW,J)*PIVOT A(IROW,J) = A(K,J) A(K,J) = -T DO 310 I=K1,N A(I,J) = A(I,J)+T*A(I,K) 310 CONTINUE 320 CONTINUE C T = AN(IROW)*PIVOT AN(IROW) = AN(K) AN(K) = -T IF( DABS(T)-EPSS ) 350,350,330 330 CONTINUE DO 340 I=K1,N AN(I) = AN(I)+T*A(I,K) 340 CONTINUE 350 CONTINUE NANS(N) = N C C BACKWARD TRANSFORMATION C AN(N) = AN(N)/A(N,N) IF( N .EQ. 1 ) GO TO 999 DO 410 J=1,NM JJ = N-J+1 JJJ = JJ-1 T = -AN(JJ) DO 400 I=1,JJJ AN(I) = AN(I)+T*A(I,JJ) 400 CONTINUE 410 CONTINUE GO TO 999 500 CONTINUE NSTOP = 2 WRITE(6,1001) AMAX , EPS , K GO TO 999 600 CONTINUE NSTOP = 3 WRITE(6,1002) N , ND GO TO 999 999 RETURN C C FORMAT STATEMENT C 1001 FORMAT(1H0,'(SUBR. HGAUSS) MATRIX IS NEARLY SINGULAR.' / 15X, 1'THE PIVOT ',D25.16,' IS LESS THAN OR EQUAL TO EPS ',D25.16,'.', 2 / 15X,'(',I4,'-TH ELIMINATING PROCESS)') 1002 FORMAT(1H0,'(SUBR. HGAUSS) N=',I4,', ND=',I4,', N MUST BE POSITIVE 1 AND LESS THAN OR EQUAL TO ND.'/15X,'RETURN WITH NO CALCULATION.') END C C SUBR. PRTOUT SUBROUTINE PRTOUT(NODNUM,NELNUM,ICAL,POWERI,POWERF,DTSTEP) IMPLICIT REAL*8(A-H,P-Z) CHARACTER*1 INCRED DIMENSION STRESS(3),DISPLA(8) C COMMON/BLOK00/ MINDX1,MINDX2 COMMON/BLKCOR/ X(25),Y(25) COMMON/BLKVEL/ VX(25),VY(25) COMMON/BLKINC/ INCDNS(16,4) COMMON/BLKRED/ INCRED(66) COMMON/BLKPRS/ VDPRES COMMON/BLKFRC/ FRCOEF COMMON/MAPROP/ S0,SA,SB,PH,PR,ERMINI,ERRATE,YNGMOD,POISNM COMMON/MATRX0/ BMATRX(16,4,8) COMMON/MATRX1/ CMATRX(16,8,8) COMMON/MATRX2/ PMATRX(16, 8) COMMON/MATRX3/ DETVOL(16) COMMON/MATRX4/ DM(4,4) COMMON/BLOPRS/ PRESUR(16) COMMON/SRBLK1/ EQRATE(16),STRATE(16,4) COMMON/SRBLK2/ EQHIST(16) C C . 節点座標，速度の印刷 WRITE( 6,100) WRITE(16,100) DO 101 I=1,NODNUM WRITE( 6,102) I,INCRED(2*I-1),INCRED(2*I ), * X(I),Y(I),VX(I),VY(I) 101 WRITE(16,102) I,INCRED(2*I-1),INCRED(2*I ), * X(I),Y(I),VX(I),VY(I) C C . 加工荷重，変形エネルギーの計算 IF(ICAL.LE.1) THEN TPOWER= POWERI+ POWERF FOLOAD= TPOWER/ VDPRES WRITE( 6,200) FOLOAD,TPOWER,POWERI,POWERF WRITE(16,200) FOLOAD,TPOWER,POWERI,POWERF ELSE WRITE( 6,500) WRITE(16,500) END IF C C . 応力，ひずみ，ひずみ速度の印刷 DO 201 I=1,NELNUM XC= 0.0 YC= 0.0 DO 251 J=1,4 XC= XC+ X(INCDNS(I,J))/4.0 251 YC= YC+ Y(INCDNS(I,J))/4.0 IF(ICAL.LE.1) THEN ER= EQRATE(I) EH= EQHIST(I) SM= YESTRS(EH,ER) CORFAC= 2.0*SM/3.0/ER STRESS(1)= CORFAC*STRATE(I,1)+PRESUR(I) STRESS(2)= CORFAC*STRATE(I,2)+PRESUR(I) STRESS(3)= CORFAC*STRATE(I,3) WRITE( 6,202) I,XC,YC,(STRESS(J),J=1,3),PRESUR(I),EH,ER WRITE(16,202) I,XC,YC,(STRESS(J),J=1,3),PRESUR(I),EH,ER ELSE EXX= 0.0 EYY= 0.0 EXY= 0.0 ESS= 0.0 DO 510 JNUM=1,4 KENODE= INCDNS(I,JNUM) DISPLA(2*JNUM-1)= VX(KENODE) 510 DISPLA(2*JNUM )= VY(KENODE) DO 520 LNUM=1,8 EXX= EXX+ BMATRX(I,1,LNUM)*DISPLA(LNUM) EYY= EYY+ BMATRX(I,2,LNUM)*DISPLA(LNUM) EXY= EXY+ BMATRX(I,3,LNUM)*DISPLA(LNUM) 520 ESS= ESS+ BMATRX(I,4,LNUM)*DISPLA(LNUM) SXX= DM(1,1)*EXX+DM(1,2)*EYY+DM(1,3)*EXY+DM(1,4)*ESS SYY= DM(2,1)*EXX+DM(2,2)*EYY+DM(2,3)*EXY+DM(2,4)*ESS SXY= DM(3,1)*EXX+DM(3,2)*EYY+DM(3,3)*EXY+DM(3,4)*ESS SSS= DM(4,1)*EXX+DM(4,2)*EYY+DM(4,3)*EXY+DM(4,4)*ESS IF(ICAL.EQ.2) SSS= POISNM*YNGMOD/ * (1.0+POISNM)/(1.0-2.0*POISNM)*(EXX+EYY) WRITE( 6,501) I,XC,YC,SXX,SYY,SXY,SSS,EXX,EYY,EXY,ESS WRITE(16,501) I,XC,YC,SXX,SYY,SXY,SSS,EXX,EYY,EXY,ESS END IF 201 CONTINUE C C . 節点座標の更新 DO 301 I=1,NODNUM X(I)= X(I)+VX(I)*DTSTEP 301 Y(I)= Y(I)+VY(I)*DTSTEP C C . 相当ひずみの積分 DO 401 I=1,NELNUM 401 EQHIST(I)= EQHIST(I)+EQRATE(I)*DTSTEP C C . フォーマット文 100 FORMAT(//' @@@ 節点座標，速度／変位 @@@', * //' Num bc X Y', * ' Vx Vy') 102 FORMAT(I10,3X,2A1,3X,2D12.4,3X,2D12.4) 200 FORMAT(//' @@@ 応力，ひずみ，ひずみ速度 @@@', * ' 加工荷重=',D12.4, * /' Tp=',D12.4,' Ip=',D12.4,' Fp=',D12.4, * //' Num X Y', * ' Sxx Syy Sxy', * ' Prs Eqh Eqr') 202 FORMAT(I10,3X,2D12.4,3X,4D12.4,3X,D12.4,3X,D12.4) 500 FORMAT(//' @@@ 応力，ひずみ @@@', * //' Num X Y', * ' Sxx Syy Sxy Szz', * ' Exx Eyy Exy Ezz') 501 FORMAT(I10,2D10.2,3X,4D12.4,3X,4D10.2) C RETURN END C C SUBR. BMACAL SUBROUTINE BMACAL(ICAL,XE,YE,XINTPN,YINTPN,GZAI,EATA, * SHPFN0,SHPFN1,SHPFN2,BMATMP,DETTMP) IMPLICIT REAL*8(A-H,P-Z) DIMENSION XE(4),YE(4),GZAI(4),EATA(4), * SHPFN0(4),SHPFN1(4),SHPFN2(4),BMATMP(4,8),A(2,2) C C . 形状関数の計算（４節点双１次要素） DO 10 JNUM=1,4 SHTMP1= (1.0+XINTPN*GZAI(JNUM)) SHTMP2= (1.0+YINTPN*EATA(JNUM)) SHPFN0(JNUM)= SHTMP1 *SHTMP2 /4.0 SHPFN1(JNUM)= GZAI(JNUM)*SHTMP2 /4.0 SHPFN2(JNUM)= SHTMP1 *EATA(JNUM)/4.0 10 CONTINUE C C . ［Ｂ］マトリックスの計算 SX0 = 0.0 SX1 = 0.0 SX2 = 0.0 SY1 = 0.0 SY2 = 0.0 DO 21 JNUM=1,4 XX= XE(JNUM) YY= YE(JNUM) SX0 = SX0 + SHPFN0(JNUM)*XX SX1 = SX1 + SHPFN1(JNUM)*XX SX2 = SX2 + SHPFN2(JNUM)*XX SY1 = SY1 + SHPFN1(JNUM)*YY SY2 = SY2 + SHPFN2(JNUM)*YY 21 CONTINUE A(1,1)= SX1 A(1,2)= SY1 A(2,1)= SX2 A(2,2)= SY2 CALL INVERS(A,DETTMP) DO 22 JNUM=1,4 DS1= A(1,1)*SHPFN1(JNUM)+ A(1,2)*SHPFN2(JNUM) DS2= A(2,1)*SHPFN1(JNUM)+ A(2,2)*SHPFN2(JNUM) BMATMP(1,2*JNUM-1)= DS1 BMATMP(1,2*JNUM )= 0.0 BMATMP(2,2*JNUM-1)= 0.0 BMATMP(2,2*JNUM )= DS2 BMATMP(3,2*JNUM-1)= DS2/ 2.0 BMATMP(3,2*JNUM )= DS1/ 2.0 BMATMP(4,2*JNUM-1)= 0.0 BMATMP(4,2*JNUM )= 0.0 22 IF(ICAL.EQ.1.OR.ICAL.EQ.3) BMATMP(4,2*JNUM-1)= SHPFN0(JNUM)/ SX0 C RETURN END C C SUBR. INVERS SUBROUTINE INVERS(A,DET) IMPLICIT REAL*8(A-H,P-Z) DIMENSION A(2,2),WK(2,2) C C . 逆行列の計算 DET= A(1,1)*A(2,2)-A(1,2)*A(2,1) DO 1 I=1,2 DO 1 J=1,2 1 WK(I,J)= A(I,J) A(1,1)= WK(2,2)/ DET A(1,2)=-WK(1,2)/ DET A(2,1)=-WK(2,1)/ DET A(2,2)= WK(1,1)/ DET RETURN END