Abaqus热源子程序Dflux教程

有限元法模拟焊接过程，要根据热源分布数学模型来描述焊接过程的输入热源，常用的有两种：

1、Gauss模型（高斯热源）

熔池深度小，可近似认为热量从表面传入，面热源

假设热源沿x方向以速度v移动，热源表达式则为：

Q为热源有效功率

2、Goldak模型（双椭球热源）

熔池深度大，考虑热源在深度方向的变化，体热源

假设热源沿x方向以速度v移动，热源表达式则为：

其中，f1，f2表示热量在原点前后两部分间分配的比例，f1+f2=2

假设对所有step使用同一热源（可以针对不同step设置不同热源，可添加parameter）

      SUBROUTINE DFLUX(FLUX,SOL,JSTEP,JINC,TIME,NOEL,NPT,COORDS,
     &JLTYP,TEMP,PRESS,SNAME)

C     FLUX：FLUX(1)流经某点的热量,FLUX(2)热量对温度的变化率
C     SOL：所求解变量的估算值，如热分析中表示某点温度估算值
C     JSTEP：载荷步
C     JINC：增量步
C     TIME：TIME(1)当前载荷步时间，TIME(2)当前总时间
C     NOEL：单元号
C     NPT：单元积分点号
C     COORDS：坐标数组
C     TEMP：当前温度
C     PRESS：当前等效应力
C     SNAME：面热源时，surface的名称

      INCLUDE 'ABA_PARAM.INC'

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

      CHARACTER*80 SNAME

C	wu,焊接电压
C	wi,焊接电流
C	effi,焊接效率系数
C	q,电弧有效热功率W
	wu=16
	wi=80
	effi=0.5
    q=wu*wi*effi
      
C   v,热源移动速度 m/s
C   TIME(2),当前分析步总时间
	v=0.005
	d=v*TIME(2)

C   x,y,z，当前积分坐标
	x=COORDS(1)
	y=COORDS(2)
	z=COORDS(3)

C	定义起始坐标
	x0=-0.05
	y0=0
	z0=0.003

C	a1,a2,b,c为双椭球形状参数
	a1=0.0028
    a2=0.0056
	b=0.0019
	c=0.0032

C	f1,f2为热源分配系数
	f1=1
    f2=2.0-f1
	PI=3.1415926

C	焊接方程定义
	heat1=6.0*sqrt(3.0)*q/(a1*b*c*PI*sqrt(PI))*f1
    heat2=6.0*sqrt(3.0)*q/(a2*b*c*PI*sqrt(PI))*f2
      
	shape1=exp(-3.0*(x-x0-d)**2/a1**2-3.0*(y-y0)**2/b**2
   $-3.0*(z-z0)**2/c**2)
    shape2=exp(-3.0*(x-x0-d)**2/a2**2-3.0*(y-y0)**2/b**2
   $-3.0*(z-z0)**2/c**2)
      

C	JLTYP=0表示面热源
C   JLTYP=1表示体热源
	JLTYP=1

C	fortran中的关系运算符
C	.GT. 大于(>)
C	.LT. 小于(<)
C	.EQ. 等于(==)
C	.GE. 大于等于(>=)
C	.LE. 小于等于(<=)
C	.NE. 不等于(!=)
	
	if (x .GE. (x0+d)) then
		FLUX(1)=heat1*shape1
    else
        FLUX(1)=heat2*shape2
	endif
	return
	end