サブルーチンデバッグ
前に公開したもののバグフィックス版
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=0.0D0 INTEGER,SAVE:: LINC=0 INTEGER,SAVE:: LSTEP=0 C SCL :積分点における解変数 C (伝熱解析の場合は温度, 質量拡散解析の場合は濃度) の現在の推定値. C KSTEP :ステップ番号. C KINC :インクリメント番号. C TIME(1):現在のステップ時間 (非定常解析でのみ定義される). C TIME(2):現在の全時間 (非定常解析でのみ定義される). C NOEL :要素番号. C NPT :要素内部または要素表面上の積分点番号. 積分点の数と位置は, C 表面流束の場合と物体流束の場合で異なる. C COORDS:この積分点の現在の座標の配列. C 幾何学非線形性を考慮しているステップでは現在の座標. C それ以外は変形前の座標である. C JLTYP :DFLUX が呼び出されている現在の流束タイプ識別番号. C 流束のタイプは, 物体流束, 表面ベース流束, 要素ベース流束 C とすることができる. 要素ベースの表面流束の場合, C この変数は DFLUXが呼び出されている現在の要素面を識別する. C 同一要素に複数の非一様分布流束が指定されている場合の識別に C 用いることができる. 識別番号は以下のとおりである. C JLTYP 流束タイプ C 0 表面ベースの流束 C 1 BFNU C TEMP :伝熱解析における温度は,変数 SOL で引き渡されている. C PRESS :この積分点における現在の相当静水圧応力 (質量拡散解析でのみ定義される). C SNAME :表面ベースの流束 (JLTYP=0) を定義する場合の表面の名前. C 物体流束と要素ベースの表面流束の場合, C 表面の名前には空白が引き渡されてくる. C Initilize finished C FLUX(1)とFLUX(2)を定義する。 C FLUX(1): C この積分点においてモデル内部へ向かう流束の値. C 伝熱解析の場合の単位は, 表面熱流束ではJT-1 L-2, C 物体熱流束では JT-1 L-3 である. C 時間変化を指定して作用する流束を変化させる非定常伝熱解析では, C 流束の値として, 時間インクリメント終了時の値ではなく, C 時間インクリメントを通じての時間平均流束を定義しなければならない. C 質量拡散解析の場合の単位は, 表面流束ではPLT-1, C 物体熱流束では PT-1である. C FLUX(1)には, 要素ベースあるいは表面ベースの非一様分布流束の C 定義の一部として指定されている流束の値が引き渡されてくる. C 値が指定されていない場合は, 0 が引き渡されてくる. C この流束は, 出力オプションを用いて出力することはできない. C FLUX(2): C 伝熱解析ではこの積分点における熱流束の変化速度,dq/dθ. C 単位は表面熱流束ではJT-1 L-2 θ-1, C 物体熱流束では JT-1 L-3 θ-1. C 質量拡散解析ではこの積分点における質量濃度流束の変化速度,dq/dc. C 表面流束では LT-1, 物体流束ではT-1. C インクリメントにおける非線形方程式の収束の速度は, C この値を定義することにより向上する. C 熱流束が温度に強く依存する場合, C または流束が質量濃度に強く依存する場合に定義することが推奨される. C ガウス熱源用の変数の定義(単位系は他の物性値等とあわせる) C Q0 :熱量[J/sec)][JT-1] C FF :比。実際の溶け込みに合わせて調整. C FR :比。実際の溶け込みに合わせて調整. C A0 :溶接線方向の前方への熱源の長さ[mm][L] C AR :溶接線方向の熱源の長さの比(大体1.5-2) C BL :溶接線直角方向の熱源の長さ[mm][L] C CL :深さ方向の熱源の長さ[mm][L] C VW :溶接速度[mm/sec][LT-1] C X :溶接線直角方向の座標 C Y :溶接線方向の座標 C Z :高さ方向の座標 C PI :円周率 C T :溶接開始からの経過時間 Q0=1.5D5 FF=1. FR=2. A0=1.5 AR=2. BL=1.5 CL=6. VW=7.5 PI=4.0D0*ATAN(1.0D0) T0=0. C ABAQUSの全体座標からこのルーチンのための座標系への変換を行う C ワーク用の変数の定義 C XA0,YA0,ZA0:全体座標での溶接開始点の座標 C XAV,YAV,ZAV:溶接方向のベクトル C XBV,YAV,ZAV:溶接線より直上へのベクトル C szM:溶接部のメッシュの大きさ C detT:インクリメントでの時間増分の大きさ C WLEN:溶接長さ XA0=-100. YA0=0. ZA0=10. XAV=-1.0 YAV=0. ZAV=0. XBV=0. YBV=0. ZBV=1. szM=5.0D0 WLEN=200.0D0 VTMP=VW*(TIME(1)-T0) IF(WLEN.LT.VTMP) THEN FLUX(1)=0.0 FLUX(2)=0.0 GOTO 111 END IF IF(KSTEP.NE.LSTEP) THEN TLAST=0.0D0 LSTEP=KSTEP END IF detT=TIME(1)-TLAST IF(KINC.NE.LINC) THEN TLAST=TIME(1) LINC=KINC END IF 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(1)-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*.5D0 gauss_const=6.D0*SQRT(3.D0)*FTMP*Q0/(ATMP*BL*CL*PI*SQRT(PI)) gauss_x=EXP(-3.D0*X*X/(BL*BL)) gauss_x1=EXP(-3.D0*(X-szM1)**2/(BL*BL)) gauss_x2=EXP(-3.D0*(X+szM1)**2/(BL*BL)) gauss_y00=EXP(-3.D0*((Y-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y01=EXP(-3.D0*((Y-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y02=EXP(-3.D0*((Y-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y0=(gauss_y01+4*gauss_y00+gauss_y02)/6.0D0 gauss_y10=EXP(-3.D0*((Y-szM1-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y11=EXP(-3.D0*((Y-szM1-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y12=EXP(-3.D0*((Y-szM1-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y1=(gauss_y11+4*gauss_y10+gauss_y12)/6.D0 gauss_y20=EXP(-3.D0*((Y+szM1-(VW*(TT1-0.5*detT)))**2)/(ATMP*ATMP)) gauss_y21=EXP(-3.D0*((Y+szM1-(VW*(TT1-detT)))**2)/(ATMP*ATMP)) gauss_y22=EXP(-3.D0*((Y+szM1-(VW*TT1))**2)/(ATMP*ATMP)) gauss_y2=(gauss_y21+4*gauss_y20+gauss_y22)/6.D0 gauss_z=EXP(-3.D0*Z*Z/(CL*CL)) gauss_z1=EXP(-3.D0*(Z-szM1)**2/(CL*CL)) gauss_z2=EXP(-3.D0*(Z+szM1)**2/(CL*CL)) gauss_x_all=(gauss_x1+4.*gauss_x+gauss_x2)/6.D0 gauss_y_all=(gauss_y1+4.*gauss_y0+gauss_y2)/6.D0 gauss_z_all=(gauss_z1+4.*gauss_z+gauss_z2)/6.D0 gauss_all=gauss_const*gauss_x_all*gauss_z_all*gauss_y_all FLUX(1)=gauss_all FLUX(2)=0.D0 111 continue RETURN end
相変わらず解析は流れない。
