Welding simulation problem using DFLUX subrutine, Fortran

[minamir]

Hi Guys,

I've got a question regarding welding simulation. I'm modeling a plate of size 180*120*17mm and a moving heat source using Abaqus and DFLUX subroutine.
Therefore, I wrote the input file, define the properties, loads, convection and radiation, etc. I also wrote a FORTRAN code -DFLUX- for moving heat source.

Heat transfer is problematic. The heat transfers all the way to the plate surface at the the same time so basically it doesn't simulate a moving heat source. Do you have any idea what could have caused this?

Here is the core of my Fortran code (also attached an image of this code for better formatting):

pointy = COORDS(2)
pointz = COORDS(3)
! try 1
ytop = wstart(2) + hcyl
ybottom = wstart(2)
DO i = 1,n
IF (t.GE.Tstart.AND.t.LE.Tend) then
IF(pointy.GE.ybottom.AND.pointy.LE.ytop) then
r = (((pointx - wx)**2 + (pointz - wz)**2)**0.5)
IF (r.LE.rcyl) then
Flux(1) = Power/Vol
WRITE (6,*) t, Flux(1), r, Vol
EXIT
END IF
END IF
Flux(1)=0
ELSE
Flux(1)=0
END IF
END DO
Flux(2)=0.


Thank you in advance

 

[Mustaine3]

I don't see how the time is changing the coordinates of the center.

 

[minamir]

I put the entire code:


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



C Abaqus Fortran Subroutine to model the moving cylindrical heat source
C Last edited: 2015-06-08
C - Fixed error in coordinates (Now y is top direction, as in the model)

c FLUX(1) will be passed into the routine as the magnitude of
c the flux specified as part of the element-based or surface-based flux definition.
c FLUX(2) is dq/dT, the rate of change of the flux with respect to the temperature
c at this point.
c The convergence rate during the solution of the nonlinear equations in an increment
c is improved by defining this value.
c SOL Estimated value of the solution variable (temperature in a heat transfer
c analysis) at this time at this point.
c KSTEP Step number.
c KINC Increment number.
c TIME(1) Current value of step time (defined only in transient analysis).
c TIME(2) Current value of total time (defined only in transient analysis).
c NOEL Element number.
c NPT Integration point number in the element or on the element's surface.
c The integration scheme depends on whether
c this is a surface or a body flux.
c COORDS An array containing the coordinates of this point.
c JLTYP Identifies the flux type for which this call to DFLUX is being made
c TEMP Current value of temperature at this integration point (defined only
c for a mass diffusion analysis).
c
c PRESS Current value of the equivalent pressure stress at this integration point
c (defined only for a mass diffusion analysis).
c SNAME Surface name for a surface-based flux definition (JLTYP=0). For a body flux
c or an element-based surface flux
c the surface name is passed in as blank.
C

REAL Qtot, rcyl, hcyl, tweld, Lweld, Vweld, Vol, Tstart, Tend,
1 Power
REAL wx, wy, wz, pointx, pointy, pointz, ztop, zbottom, r, t, pi
INTEGER n
DIMENSION Flux(2)
DIMENSION TIME(2)
DIMENSION Coords(3)
DIMENSION wstart(3)
DIMENSION wend(3)
pi = 4*atan(1.0)

C _______________________________________
C Heat source parameters
! Total welding heat input [micro joules]
Qtot = 38e+9*0.75

! radius of torch:
rcyl = 2.4

! torch depth:
hcyl = 0.7

! number of weld pass
n=1

tweld = 26.4
Lweld = 60

! Mean welding power input (micro watt )
Power = Qtot/ tWeld
! welding speed 60mm/26.4
Vweld = Lweld/tweld
C _______________________________________

! start point

wstart(1) = 0.0
wstart(2) = 0.0
wstart(3) = -30.0

! end point

wend(1) = 0.0
wend(2) = 0.0
wend(3) = 30.0

! volume of welding cylinder
Vol = (pi*rcyl**2)*hcyl
! weld start and end times

! current time
t = TIME(2)

Tstart = 0
Tend = Tstart + tweld
! torch location

wx = wstart(1)
wy = wstart(2)
wz = wstart(3) + Vweld*t

! Coordinates of the point under consideration
pointx = COORDS(1)
pointy = COORDS(2)
pointz = COORDS(3)
! try 1
ytop = wstart(2) + hcyl
ybottom = wstart(2)
DO i = 1,n
IF (t.GE.Tstart.AND.t.LE.Tend) then
IF(pointy.GE.ybottom.AND.pointy.LE.ytop) then
r = (((pointx - wx)**2 + (pointz - wz)**2)**0.5)
IF (r.LE.rcyl) then
Flux(1) = Power/Vol
WRITE (6,*) t, Flux(1), r, Vol
EXIT
END IF
END IF
Flux(1)=0
ELSE
Flux(1)=0
END IF
END DO
Flux(2)=0.

 

[Mustaine3]

I don't debug your subroutine.

Print out much more data and check through that what is wrong.

 

[minamir]

Hi Mustain3,

- I debugged this code quite extensively using Visual Studio debugger and the output seems legitimate to me. because time is increasing in my printout. that's why I got confused what's going on here.

 

[Mustaine3]

Then print wz and r and check what this is doing. And maybe add a printout within an IF to see if it's fulfilled when expected.

 

[minamir]

Yes Mustaine3 you're right. It can't print wz and r!

 

[minamir]

hmm no it prints