Rapid Simulation of Laser Processing of Discrete Particulate Materials

Abstract

The objective of this paper is to develop a computational model and corresponding solution algorithm to enable rapid simulation of laser processing and subsequent targeted zonal heating of materials composed of packed, discrete, particles. Because of the complex microstructure, containing gaps and interfaces, this type of system is extremely difficult to simulate using continuum-based methods, such as the Finite Difference Time Domain Method or the Finite Element Method. The computationally-amenable model that is developed captures the primary physical events, such as reflection and absorption of optical energy, conversion into heat, thermal conduction through the microstructure and possible phase transformations. Specifically, the features of the computational model are (1) a discretization of a concentrated laser beam into rays, (2) a discrete element representation of the particulate material microstructure and (3) a discrete element transient heat transfer model that accounts for optical (laser) energy propagation (reflection and absorption), its conversion into heat, the subsequent conduction of heat and phase transformations involving possible melting and vaporization. A discrete ray-tracking algorithm is developed, along with an embedded, staggered, iterative solution scheme, which is needed to calculate the optical-to-thermal conversion, particle-to-particle conduction and phase-transformations, implicitly. Numerical examples are given, focusing on concentrated laser beams and the effects of surrounding material conductivity, which draws heat away from the laser contact zone, thus affecting the targeted material state.

Appendix: Geometrical Ray Theory

Following a somewhat classical analysis found in, for example, Elmore and Heald [19], Cerveny et al. [8] and others, we consider the propagation of a general disturbance, ψ, governed by a generic wave equation:

$$ \nabla^2 \psi=\frac{1}{c^2(\boldsymbol{x})}\frac{\partial^2 \psi }{\partial t^2}. $$

(7.1)

Here c(x) is a spatially varying wave speed corresponding to a general inhomogeneous medium, where c(x)=c _o in a homogeneous reference medium and where the refractive index is defined as n=c _o /c(x). Consider a trial solution of the form

$$ \psi(\boldsymbol{x},t)=A(\boldsymbol {x})e^{j(k_oS(\boldsymbol{x})-\omega t)}, $$

(7.2)

where A(x) is the amplitude of the disturbance, and where k _o =ω/c _o =2π/λ is the wave number in the reference medium. The function S(x) (dimensions of length) is known as the “Eikonal”, which in Greek means “image”. One can interpret a set of waves as simply a family of surfaces for which the values of k _o S(x) differ in increments of 2π. Substituting the trial solution into the wave equation, one obtains

$$ \begin{aligned}[b] &k_o^2A \bigl(n^2-\nabla S \cdot\nabla S \bigr) \\ &\quad{}+ jk_o \bigl(2\nabla A \cdot\nabla S+A \nabla^2 S \bigr)+ \nabla^2 A=0. \end{aligned} $$

(7.3)

There are a variety of arguments to motivate so-called “Ray Theory”. Probably the simplest is simply to require that, as k _o →∞, each of the k _o -terms, the zeroth-order k _o -term, the first-order k _o -term and the second-order k _o -term, must vanish. Applying this requirement to the second-order k _o -term yields

$$ n^2=\nabla S \cdot\nabla S=\|\nabla S \|^2. $$

(7.4)

For a uniform medium, n=const, provided ∇² A=0 and an initial plane wave surface S=const, then (7.3) implies

$$ S(\boldsymbol{x})=n(\alpha x+\beta y+\phi z), $$

(7.5)

where α, β and ϕ are direction cosines. More generally, when n≠0, then (7.4) implies

$$ \nabla S(\boldsymbol{x})=n(\boldsymbol{x}) \hat {\boldsymbol{s}}(\boldsymbol{x}), $$

(7.6)

where $\hat{\boldsymbol{s}}(\boldsymbol{x})$ is a unit (direction) vector. From elementary calculus, recall that ∇S is perpendicular to S=const. This allows for the natural definition of continuous curves, called rays, that are everywhere parallel to the local direction $\hat{\boldsymbol{s}}(\boldsymbol{x})$. Rearranging first-order k _o -term of (7.3)

$$ \frac{1}{A}\nabla A \cdot\nabla S =-\frac{1}{2} \nabla^2 S =-\frac{1}{2}\nabla\cdot(n\hat{\boldsymbol{s}}). $$

(7.7)

Recalling the directional derivative, $\frac{d (\cdot)}{ds}\stackrel{\mathrm{def}}{=}\hat{\boldsymbol {s}}\cdot\nabla (\cdot)$, we have

$$ \biggl(\frac{\nabla S}{\|\nabla S\|} \biggr)\cdot \nabla A= \biggl( \frac{\nabla S}{n} \biggr)\cdot\nabla A= \frac{dA}{ds}, $$

(7.8)

where s is the arc-length coordinate along the ray. With this definition, once S(x) is known, the component of ∇A in the $\hat {\boldsymbol{s}}(\boldsymbol{x})$ can be found from (7.7) and (7.8):

$$ \frac{1}{A}\frac{dA}{ds}=-\frac{1}{2n}\nabla\cdot(n\hat {\boldsymbol{s}}). $$

(7.9)

Thus, we are able to determine how the amplitude of the trial solution changes along a ray, but not perpendicular to the trajectory.

Geometrical “ray-tracing”, deals directly with the ray trajectories, rather than finding them as a by-product of the solution of the wave equation for the Eikonal function S and the resulting wave front. To eliminate S we need to look at the rate of change of the quantity $n\hat{\boldsymbol{s}}$ along the ray. Making repeated use of (7.6), we have

$$\begin{aligned} \frac{d(n \hat{\boldsymbol{s}})}{ds} =& \hat{\boldsymbol{s}}\cdot \nabla(\nabla S)= \frac{\nabla S}{n}\cdot\nabla(\nabla S) \\ =& \frac{1}{2n} \nabla(\nabla S \cdot \nabla S)= \frac{1}{2n} \nabla n^2=\nabla n, \end{aligned}$$

(7.10)

where $\frac{d (\cdot)}{ds}\stackrel{\mathrm{def}}{=}\hat {\boldsymbol{s}}\cdot \nabla(\cdot)$. The previous equation allows us to find the trajectories of a ray ($\hat{\mbox {\boldmath$s$}}$), given only the refractive index n(x) and the initial direction $\hat{\boldsymbol{s}}_{i}$ of the desired ray.

Remark

A more general derivation of the Eikonal equation can be found in a variety of textbooks, for example, Cerveny et al. [8], and starts by assuming a trial solution of the form

$$ \psi(\boldsymbol{x},t)=A(\boldsymbol{x})\varPhi\bigl(t-\varLambda (\boldsymbol{x}) \bigr) $$

(7.11)

where Λ is an Eikonal function, and the waveform function α is assumed to be of high frequency.¹⁶ This function is then substituted into the wave equation to yield

$$ \nabla^2A \varPhi+2\nabla A\cdot\nabla\varPhi+A \nabla^2\varPhi=\frac {1}{c^2}A\frac{\partial^2\varPhi}{\partial t^2}. $$

(7.12)

After using the chain rule of differentiation, this can be written as

$$\begin{aligned} &{\frac{\partial^2 \varPhi}{\partial\varLambda^2} A \biggl(\nabla\varLambda\cdot\nabla \varLambda- \frac{1}{c^2} \biggr)} \\ &{\quad{}+ \frac{\partial\varPhi}{\partial\varLambda} \bigl( 2\nabla A\cdot \nabla\varLambda+A \nabla^2 \varLambda\bigr) +\varPhi \nabla^2 A=0.} \end{aligned}$$

(7.13)

Analogous to the special case considered before, to motivate so-called “Ray Theory” one requires that the coefficients of $\frac{\partial^{2} \varPhi}{\partial\varLambda^{2}}$, $\frac{\partial\varPhi}{\partial\varLambda}$ and Φ are satisfied separately, in other words, the following hold

$$ \nabla\varLambda\cdot\nabla\varLambda-\frac{1}{c^2}=0 $$

(7.14)

and

$$ 2\nabla A\cdot\nabla\varLambda+A \nabla^2 \varLambda=0 $$

(7.15)

and

$$ \nabla^2 A=0. $$

(7.16)