やっぱり…
今日も論文書かなきゃいけないのにサブルーチンの機能だけ増えた。
増えたのは熱源サブルーチンのインクリメントにおける時間増分の自動認識(実際はこんな大層なものではない)と、ひずみ応力変換サブルーチンのすっかり忘れていたせん断ひずみへの対応。
SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS, 1JLTYP,TEMP,PRESS,SNAME) C INCLUDE 'ABA_PARAM.INC' C DIMENSION FLUX(2), TIME(2), COORDS(3) CHARACTER*80 SNAME REAL*8,SAVE:: TLAST C Initilize finished C Q0 :熱量[JL-3T-1] C FF :比。実際の溶け込みに合わせて調整 C FR :比。実際の溶け込みに合わせて調整 C A0 :溶接線方向の前方への熱源の長さ[L] C AR :溶接線方向の熱源の長さの比 C BL :溶接線直角方向の熱源の長さ[L] C CL :深さ方向の熱源の長さ[L] C VW :溶接速度[LT-1] C X :溶接線直角方向の座標 C Y :溶接線方向の座標 C Z :高さ方向の座標 C PI :円周率 C T0 :溶接開始の時間 C szM:溶接部のメッシュの大きさ C detT:インクリメントでの時間増分の大きさ Q0=1500. FF=1. FR=4. A0=0.15 AR=2. BL=0.3 CL=0.7 VW=0.75 PI=4.0E0*ATAN(1.0D0) T0=0. szM=.2 detT=TIME(2)-TLAST TLAST=TIME(2) C ABAQUSの全体座標からこのルーチンのための座標系への変換を行う C XA0,YA0,ZA0:全体座標での溶接開始点の座標 C XAV,YAV,ZAV:溶接方向のベクトル C XBV,YAV,ZAV:溶接線より直上へのベクトル XA0=0 YA0=0.424264 ZA0=1.32426 XAV=1. YAV=0 ZAV=0 XBV=0 YBV=.424264 ZBV=.424264 C 座標計算 C 与えられたベクトルの大きさを1にそろえる LTEMP=SQRT(XAV*XAV+YAV*YAV+ZAV*ZAV) XAV=XAV/LTEMP YAV=YAV/LTEMP ZAV=ZAV/LTEMP LTEMP=SQRT(XBV*XBV+YBV*YBV+ZBV*ZBV) XBV=XBV/LTEMP YBV=YBV/LTEMP ZBV=ZBV/LTEMP C 与えられたベクトルより残りの座標軸となるベクトルを計算 XCV=YAV*ZBV-ZAV*YBV YCV=ZAV*XBV-XAV*ZBV ZCV=XAV*YBV-YAV*XBV C 座標軸を平行移動 XTMP=COORDS(1)-XA0 YTMP=COORDS(2)-YA0 ZTMP=COORDS(3)-ZA0 C 与えられたベクトルが基底になるように座標軸を回転 X=XCV*XTMP+XAV*YTMP+XBV*ZTMP Y=YCV*XTMP+YAV*YTMP+YBV*ZTMP Z=ZCV*XTMP+ZAV*YTMP+ZBV*ZTMP C 座標計算はここまで C 熱源の前方か否か TT1=TIME(2)-T0 TT=Y-(VW*TT1) IF (TT.LT.0) THEN ATMP=A0 FTMP=FF ELSE ATMP=A0*AR FTMP=FR END IF C 最終出力 szM1=szM*.5 gauss_const=6.E0*SQRT(3.)*FTMP*Q0/(ATMP*BL*CL*PI*SQRT(PI)) gauss_x=EXP(-3.E0*X*X/(BL*BL)) gauss_x1=EXP(-3.E0*(X-szM1)**2/(BL*BL)) gauss_x2=EXP(-3.E0*(X+szM1)**2/(BL*BL)) gauss_y00=EXP(-3.E0*((Y-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y01=EXP(-3.E0*((Y-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y02=EXP(-3.E0*((Y-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y0=(gauss_y01+4*gauss_y00+gauss_y02)/6 gauss_y10=EXP(-3.E0*((Y-szM1-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y11=EXP(-3.E0*((Y-szM1-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y12=EXP(-3.E0*((Y-szM1-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y1=(gauss_y11+4*gauss_y10+gauss_y12)/6 gauss_y20=EXP(-3.E0*((Y+szM1-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y21=EXP(-3.E0*((Y+szM1-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y22=EXP(-3.E0*((Y+szM1-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y2=(gauss_y21+4*gauss_y20+gauss_y22)/6 gauss_z=EXP(-3.E0*Z*Z/(CL*CL)) gauss_z1=EXP(-3.E0*(Z-szM1)**2/(CL*CL)) gauss_z2=EXP(-3.E0*(Z+szM1)**2/(CL*CL)) gauss_x_all=(gauss_x1+4.*gauss_x+gauss_x2)/6. gauss_y_all=(gauss_y1+4.*gauss_y0+gauss_y2)/6. gauss_z_all=(gauss_z1+4.*gauss_z+gauss_z2)/6. gauss_all=gauss_const*gauss_x_all*gauss_z_all*gauss_y_all FLUX(1)=gauss_all FLUX(2)=0.D0 RETURN end
*** READ FROM .fil *** C SUBROUTINE HKSMAIN C INCLUDE 'aba_param.inc' CHARACTER FNAME*80,FNAMEO*80 C 最大レコード数 PARAMETER(NMAX=1000000) C 最大エレメント数 PARAMETER(MAXE=1000000) C キーのためのオフセット PARAMETER(KOFFSET=2) C ヤング率[M L-1 T-2] PARAMETER(YOUNGR=210.0D+7) C ポアソン比 PARAMETER(POISON=0.3D+0) DIMENSION ARRAY(513),JRRAY(NPRECD,513),LRUNIT(2,1) DIMENSION PEX(MAXE),PEY(MAXE),PEZ(MAXE),PEA(MAXE),PES(MAXE) DIMENSION PEXY(MAXE),PEXZ(MAXE),PEYZ(MAXE) DIMENSION NELN(MAXE) EQUIVALENCE(ARRAY(1),JRRAY(1,1)) C C 入力ファイル名と出力ファイル名 C 入力ファイル名では拡張子は省略 C WRITE(6,*)'Enter the name of the input file (w/o extension):' READ(5,'(A)') FNAME WRITE(6,*)'Enter the name of the output file:' READ(5,'(A)') FNAME0 C DO 10 I=1,MAXE PEX(I) = 0.D+00 PEY(I) = 0.D+00 PEZ(I) = 0.D+00 PEXY(I) = 0.D+00 PEXZ(I) = 0.D+00 PEYZ(I) = 0.D+00 PES(I) = 0.D+00 PEA(I) = 0.D+00 NELN(I) = 0 10 CONTINUE NTHRES=0 C *** READ FROM .fil *** C 初期化 C 読み込む結果ファイルの数 NRU = 1 C デバイス番号 LRUNIT(1,1) = 8 C 結果ファイルがバイナリの場合2,アスキーの場合1 44 CONTINUE WRITE(6,*)'Enter the format of input file (1-ASCII, 2-binary):' READ(5,*)LRUNIT(2,1) IF (LRUNIT(2,1).NE. 1 .AND. LRUNIT(2,1) .NE. 2) THEN WRITE(6,*) 'ERROR! This number must be 1 or 2' GOTO 44 END IF C 書き出し形式、読み込みの場合は何の意味も持たない LOUTF = 0 CALL INITPF(FNAME,NRU,LRUNIT,LOUTF) C デバイス番号 JUNIT = 8 CALL DBRNU(JUNIT) C DO 100 J=1,NMAX C レコードの順次読み出し CALL DBFILE(0,ARRAY,JRCD) C EOFの検出 IF (JRCD .NE. 0) GOTO 110 KEY = JRRAY(1,KOFFSET) C エレメント番号 IF (KEY .EQ. 1) THEN NELNO = JRRAY(1,KOFFSET+1) NELN(NELNO) = NELNO END IF C C 塑性ひずみPE IF (KEY .EQ. 22) THEN PEX(NELNO) = ARRAY(KOFFSET+1) PEY(NELNO) = ARRAY(KOFFSET+2) PEZ(NELNO) = ARRAY(KOFFSET+3) PEXY(NELNO) = ARRAY(KOFFSET+4) PEXZ(NELNO) = ARRAY(KOFFSET+5) PEYZ(NELNO) = ARRAY(KOFFSET+6) PES(NELNO) = ARRAY(KOFFSET+7) PEA(NELNO) = ARRAY(KOFFSET+9) END IF C 100 CONTINUE 110 CONTINUE KEND = J-1 WRITE(6,*) 'KEND=',KEND WRITE(6,*) ' ' C *** WRITE TO FILE *** C OPEN(51,FILE=FNAMEO) C C WRITE(51,*)'**Initial conditions,TYPE=STRESS' DO K=1,MAXE IF (NELN(K) .NE. 0) THEN PEAL=PEX(K)+PEY(K)+PEZ(K) GR=YOUNGR/(2*(1+POISON)) SX=YOUNGR/(1+POISON)*(PEX(K)- 1 (POISON/(1.0-2.0*POISON))*(PEAL))) SY=YOUNGR/(1+POISON)*(PEY(K)- 1 (POISON/(1.0-2.0*POISON))*(PEAL))) SZ=YOUNGR/(1+POISON)*(PEZ(K)- 1 (POISON/(1.0-2.0*POISON))*(PEAL))) SXY=GR*PEXY(K) SXZ=GR*PEXA(K) SYZ=GR*PEYZ(K) WRITE(51,*)NELN(K),',',SX,',',SY,',',SZ,',',SXY,',',SXZ,',', 1 SZY END IF END DO C WRITE(6,*) ' ' WRITE(51,*) ' ' STOP END
