Welding Simulation Along an Elliptical Path in Abaqus Using a DFLUX Subroutine

Moving-heat-source welding simulation is one of the most powerful capabilities in Abaqus for predicting residual stress, distortion, and the thermal history of welded components. Most tutorials demonstrate a heat source traveling in a straight line or a circle, because those paths are easy to parameterize. Real fabrication, however, often follows more complex geometries- flanges, nozzle attachments, and sealing welds that trace an ellipse.

An elliptical weld path introduces a subtlety that trips up many engineers: the torch does not move at a constant angular rate. If you naively write angle = ω·t, the heat source will race through the sharply curved ends of the ellipse and crawl along the flatter sides, giving you the wrong heat input per unit length and an unphysical weld pool.

This article explains the correct formulation and provides a complete, tested DFLUX subroutine that keeps the travel speed constant along an ellipse.

Why an Ellipse Is Not Just a Stretched Circle

For a circle of radius $R$R, arc length and angle are proportional: $s = R\theta$. Move at angular speed $\omega = v/R$ω=v/R and the linear speed is automatically constant.

For an ellipse with semi-axes $a_e$​ and $b_e$​, the position is:

$$y_c(\theta) = a_e\cos\theta, \qquad z_c(\theta) = b_e\sin\theta$$

The speed at which the point travels as $\theta$ increases is the magnitude of the derivative:

$$\left|\frac{d\mathbf{r}}{d\theta}\right| = \sqrt{a_e^2\sin^2\theta + b_e^2\cos^2\theta}$$

This quantity changes with position — it equals $b_e$be​ at the ends of the major axis and $a_e$ae​ along the minor axis. Because there is no closed-form inverse of the elliptic arc-length integral, we must handle it numerically.

The physical requirement is a constant travel speed $v$, so the arc length covered must satisfy:

$$s(\theta) = \int_0^{\theta}\sqrt{a_e^2\sin^2\phi + b_e^2\cos^2\phi}\; d\phi = v\,t$$

The simulation therefore has to:

1. Compute the total perimeter of the ellipse.
2. Convert the current analysis time into an arc length, $s = v\,t$.
3. Invert $s(\theta)$ to recover the parametric angle $\theta$.
4. Build a local weld frame (travel direction + transverse + depth) at that point.
5. Evaluate the Goldak double-ellipsoid flux.

Step 1 – Compute the Perimeter Numerically

We integrate the arc-length integrand using Simpson’s rule, which converges very quickly for this smooth periodic function. With 720 panels, the perimeter is accurate to roughly $10^{-10}$.

Step 2 – Wrap the Arc Length

For repeated laps, the traveled arc length is wrapped back into one perimeter:

$$s_{\text{wrapped}} = s – P\left\lfloor \frac{s}{P}\right\rfloor$$

Step 3 – Invert Arc Length to Angle

We first march through the precomputed sub-arcs to bracket the correct interval, then refine with Newton–Raphson using the exact derivative $f'(\theta) = |\mathbf{r}'(\theta)|$f′(θ)=∣r′(θ)∣.

This converges to machine precision in a handful of iterations.

Step 4 – Build the Local Weld Frame

At the source center, the unit tangent (travel direction) and normal (transverse direction) are:

$$\mathbf{t} = \frac{(-a_e\sin\theta,\; b_e\cos\theta)}{|\mathbf{r}’|}, \qquad \mathbf{n} = \frac{(b_e\cos\theta,\; a_e\sin\theta)}{|\mathbf{r}’|}$$

Any material point is then projected onto this frame to get the Goldak local coordinates $\xi$ (along travel), $\eta$ (transverse), and $\zeta$ (depth).

Understanding the Goldak Parameters

The double-ellipsoid model splits the heat source into a front quadrant and a rear quadrant, each with its own length ($c_f$, $c_r$​) and fraction ($f_f$​, $f_r$​, with $f_f + f_r = 2$). This asymmetry captures the steep temperature gradient ahead of the arc and the longer trailing tail behind it.

The normalization constant $6\sqrt{3}/(abc\,\pi\sqrt{\pi})$ guarantees the volume integral of the flux equals the applied power, so energy is conserved regardless of path shape.

Setting Up the Analysis in Abaqus

1. Define a transient heat transfer step (*Heat Transfer, Transient) with a time increment small enough that the source moves a fraction of its own length per increment — typically $\Delta t \le c_f / v$Δt≤cf​/v.
2. Apply a body flux on the weld region: *Dflux, elset=WELD, BFNU. The NU (non-uniform) flag routes the load to your subroutine.
3. Link the subroutine at job submission with abaqus job=weld user=dflux.for.
4. Add radiation and convection on exposed surfaces to recover realistic cooling.
5. Set the total step time to cover the desired number of laps. For this example the perimeter is about $0.2549$0.2549 m, so one full lap at $v = 0.002$v=0.002 m/s takes roughly 127.4 s.

Adapting the Subroutine to Your Model

• Different plane orientation. The code assumes the ellipse lies in the 2–3 plane with COORDS(1) as depth. If your weld sits in the X–Y plane, use COORDS(1) and COORDS(2) for the path and COORDS(3) for depth.
• Different start point. Replace STG = V*T with STG = S0 + V*T, where S0 is the arc-length offset of your starting position.
• Off-center ellipse. Add the center coordinates when computing YC and ZC.
• Reverse direction. Negate the travel speed V or the tangent components.

Conclusion

Simulating a welding heat source on an elliptical path requires more than stretching a circular model – it demands correct handling of arc length so the torch maintains constant travel speed. By combining Simpson integration for the perimeter, Newton–Raphson inversion for the angle, and the Goldak double-ellipsoid flux in a clean local frame, this DFLUX subroutine delivers physically accurate heat input for any elliptical weld.

This approach generalizes directly to other parametric curves – splines, racetrack paths, or scan patterns in additive manufacturing ; by swapping the position functions and their derivatives while keeping the same arc-length machinery.