またもやサブルーチンが増えた
線熱源ということでノードに入熱していたわけだが、それだとあまりにも熱勾配が大きすぎるということでまたしても入熱用のサブルーチンを書いた。これで応力解析が流れてくれるとうれしいのだが。
C 線熱源用 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 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/(mm3*sec)][JL-3T-1] C A0 :熱源の大きさ[L] C PI :円周率 Q0=5.5D5 A0=3 PI=4.0D0*ATAN(1.0D0) C 時間変化 IF (TIME(1).GT.2.0D0) THEN GOTO 1234 ELSE IF (TIME(1).LT.5.0D-1) THEN QTMP=Q0*TIME(1)*2.0D0 ELSE IF (TIME(1).LT.1.5D0) THEN QTMP=Q0 ELSE QTMP=Q0*(2.0D0-TIME(1))*2.0D0 END IF C 溶接線の定義 C X0,Y0,Z0:溶接線の開始点 C WLEN:溶接線の長さ C XV,YV,ZV:溶接線の向き X0= Y0= Z0= WLEN=200.0D0 XV= YV= ZV= C 与えられたベクトルの大きさを1にそろえる LTEMP=SQRT(XV*XV+YV*YV+ZV*ZV) XV=XV/LTEMP YV=YV/LTEMP ZV=ZV/LTEMP C溶接線との距離を計算 PX=COORDS(1) PY=COORDS(2) PZ=COORDS(3) X2=X0+XV*WLEN Y2=Y0+YV*WLEN Z2=Z0+ZV*WLEN X1=(X0+X2)*0.5D0 Y1=(Y0+Y2)*0.5D0 Z1=(Z0+Z2)*0.5D0 EPS=(WLEN*5.0D-3)**2 TLEN=(X0-X2)**2+(Y0-Y2)**2+(Z0-Z2)**2 DO WHILE (TLEN.GT.EPS) TLEN0=(PX-X0)**2+(PY-Y0)**2+(PZ-Z0)**2 TLEN2=(PX-X2)**2+(PY-Y2)**2+(PZ-Z2)**2 IF (TLEN0.GT.TLEN2) THEN X0=X1 Y0=Y1 Z0=Z1 ELSE IF (TLEN2.GT.TLEN0) THEN X2=X1 Y2=Y1 Z2=Z1 ELSE GOTO 123 END IF TLEN=(X0-X2)**2+(Y0-Y2)**2+(Z0-Z2)**2 X1=(X0+X2)*0.5D0 Y1=(Y0+Y2)*0.5D0 Z1=(Z0+Z2)*0.5D0 END DO 123 CONTINUE TLEN=(X1-PX)**2+(Y1-PY)**2+(Z1-PZ)**2 C 最終出力 gauss_const=6.D0*SQRT(3.D0)*QTMP/(A0*PI*SQRT(PI)) gauss_len=EXP(-3.D0*TLEN/(A0*A0)) gauss_all=gauss_const*gauss_len FLUX(1)=gauss_all FLUX(2)=0.D0 1234 CONTINUE RETURN end
