求问我这个移动高斯柱体热源哪儿出问题了,请大家指正

     SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,

     

1 JLTYP,TEMP,PRESS,SNAME)

      

INCLUDE 'ABA_PARAM.INC'

      

parameter(one=1.d0)

      

DIMENSION COORDS(3),FLUX(2),TIME(2)

      

CHARACTER*80 SNAME

      

C     热源功率（热源瞬时给焊件的热能）Q,热效率n；热源半径Rh，单位m ;热源穿透深度80微米，0.00008m

      

Q=170

      

n=1

      

Rh=0.0001

      

PI=3.1415

      

h=0.00008

     

 

C     移动速度v=0.05m/s

      

v=0.05

      

x=COORDS(1)

     

y=COORDS(2)

      

z=COORDS(3)

 

C     定义移动热源起始点。

      

x0=-0.001

      

y0=0

      

z0=0

      

d=v*TIME(1)

C     JLTYP＝1，表示为体热源      

      

JLTYP=1 

      

C     计算热源热流密度，FLUX,单位J/(m^2*s)

      

R=sqrt((x-x0-d)**2+(z-z0)**2)

      

heat=3*Q*n/(h*PI*Rh**2)

      

shape1=exp(-3*R**2/Rh**2)

     

if (JSTEP.eq.one) then

      

FLUX(1)=heat*shape1

      

    

      

RETURN

      

END