A greedy algorithm for optimal heating in powder-bed-based additive manufacturing

TL;DR

This paper presents a greedy optimization algorithm for controlling heat input in powder-bed additive manufacturing, improving temperature management and melt uniformity during the process.

Contribution

It introduces a novel greedy control scheme that divides the complex heat optimization into smaller problems, enhancing melt pool control in additive manufacturing.

Findings

Successfully manages heat accumulation at turning points

Achieves more uniform melt depths

Resolves heat conduction issues

Abstract

Powder-bed-based additive manufacturing involves melting of a powder bed using a moving laser or electron beam as a heat source. In this paper, we formulate an optimization scheme that aims to control this type of melting. The goal consists of tracking maximum temperatures on lines that run along the beam path. Time-dependent beam parameters (more specifically, beam power, spot size, and speed) act as control functions. The scheme is greedy in the sense that it exploits local properties of the melt pool in order to divide a large optimization problem into several small ones. As illustrated by numerical examples, the scheme can resolve heat conduction issues such as concentrated heat accumulation at turning points and non-uniform melt depths.

Tables

Table 1. Table 1 . Parameter values in Example 6.1 .

PDE specific input
Thermal conductivity	$\lambda$	$20$	$\left(W/\mathrm{mK}\right)$
Thermal diffusivity	$\kappa$	$8.45e\text{-}6$	$\left(m^2/s\right)$
Initial temperature	$u_{\mathrm{init}}$	$1000$	$\left(K\right)$
Absorbed beam power	$P$	$100$	$\left(W\right)$
Greedy algorithm specific input
Reference surface temperature	$u_{\mathrm{melt}}$	$1800$	$\left(K\right)$
Reference melt temperature	$u_{\mathrm{surf}}$	$2800$	$\left(K\right)$
Secondary beam position, width	$w$	$100$	$\left(\mu m\right)$
Secondary beam position, depth	$d$	$50$	$\left(\mu m\right)$
Weight 1	$W_1$	$0.7$
Weight 2	$W_2$	$0.3$
Window size	$q-p+1$	$5$
Segments frozen in each iteration	$r$	$1$
Initial spot size on segment $k$	${\sigma }_k$	$0.2\forall k$	$\left(\mathrm{mm}\right)$
Initial speed on segment $k$	$v_k$	$0.5\forall k$	$\left(m/s\right)$
Bounds, spot size	$({\sigma }_{\min },{\sigma }_{\max })$	$(1e\text{-}2,1e0)$	$\left(\mathrm{mm}\right)$
Bounds, speed	$(v_{\min },v_{\max })$	$(1e1,1e4)$	$\left(\mathrm{mm}/s\right)$

Equations

$${\mathcal{C}}^{\mathrm{s}}=\{{\mathbf{x}}^{\mathrm{s}}(t):t\in \mathcal{T}\},$$

$$\Phi=\Phi(x,y,t)=\frac{P(t)}{2\pi\sigma(t)^{2}}\,\mathrm{exp}\Big{(}-\tfrac{\left(x-x^{\mathrm{s}}\left(t\right)\right)^{2}+\left(y-y^{\mathrm{s}}\left(t\right)\right)^{2}}{2\sigma(t)^{2}}\Big{)}.$$

$$(t_{n-1},t_{n}]=\big{(}t_{n}^{\mathrm{i}},t_{n}^{\mathrm{f}}\big{]},\,\,\,\,n=1,2,\ldots,N,$$

$$\left(P(t),\sigma(t),v(t)\right)=(P_{n},\sigma_{n},v_{n})\,\,\,\mathrm{if}\,\,t\in\big{(}t_{n}^{\mathrm{i}},t_{n}^{\mathrm{f}}\big{]},$$

$$\rho c_p\frac{\partial u}{\partial t}-\nabla \cdot (\lambda \nabla u)$$

$$(\lambda \nabla u)\cdot \hat{z}$$

$$u(\cdot ,0)$$

$$u(\mathbf{x},t)$$

$$\sigma (t)$$

$$v(t)$$

$$\gamma ({\mathbf{x}}^{\mathrm{s}})=\sum _{k=1}^{n-1}|{\mathbf{x}}_k^{\mathrm{f}}-{\mathbf{x}}_k^{\mathrm{i}}|+|{\mathbf{x}}^{\mathrm{s}}-{\mathbf{x}}_n^{\mathrm{i}}|,\mathrm{if}{\mathbf{x}}^{\mathrm{s}}\in {\ell }_n^s,n=1,2,\dots ,N.$$

$$\begin{array}{rl}\sigma (\gamma ({\mathbf{x}}^{\mathrm{s}}))& =\sum _{k=1}^N{\sigma }_k{\chi }_{(\gamma ({\mathbf{x}}_k^{\mathrm{i}}),\gamma ({\mathbf{x}}_k^{\mathrm{f}})]},\\ v(\gamma ({\mathbf{x}}^{\mathrm{s}}))& =\sum _{k=1}^N{v}_k{\chi }_{(\gamma ({\mathbf{x}}_k^{\mathrm{i}}),\gamma ({\mathbf{x}}_k^{\mathrm{f}})]}.\end{array}$$

$$\omega (t)=\{\mathbf{x}\in \Omega :u(\mathbf{x},t)\ge u_{\mathrm{melt}}\}.$$

$${\mathcal{C}}^{\mathrm{wd}}=\{\mathfrak{F}({\mathbf{x}}^{\mathrm{s}}(t)):t\in \mathcal{T}\}$$

$${\mathcal{C}}^{\mathrm{wd}}=\{(x^{\mathrm{s}}(t)+w\sin {\theta }_1,y^{\mathrm{s}}(t)-w\cos {\theta }_1,-d):t\in \mathcal{T}\}.$$

$$\mathcal{M}(\mathbf{x};\mathbf{d})=\underset{t\in \mathcal{T}}{\max }\{u\left(\mathbf{x},t;\mathbf{d}\right)\}.$$

$$f(\bm{d})=\big{\{}f_{1}(\bm{d}),f_{2}(\bm{d})\big{\}},$$

$$\begin{split}f_{1}(\bm{d})&=\int_{\mathcal{C}^{\mathrm{wd}}}\alpha(\bm{x}^{\mathrm{wd}})\cdot\big{(}\mathcal{M}(\bm{x}^{\mathrm{wd}};\bm{d})-u_{\mathrm{melt}}\big{)}^{2}\,\mathrm{d}x,\\ f_{2}(\bm{d})&=\int_{\mathcal{C}^{\mathrm{s}}}\alpha(\bm{x}^{\mathrm{s}})\cdot\big{(}\mathcal{M}(\bm{x}^{\mathrm{s}};\bm{d})-u_{\mathrm{surf}}\big{)}^{2}\,\mathrm{d}x.\end{split}$$

$$\alpha (\mathbf{x})=\{\begin{array}{l}0\text{ if }\mathbf{x}\text{ is near the start or end of a hatch line,}\\ 1\text{ otherwise}\end{array}$$

$$\begin{array}{r}\begin{array}{rl}& \mathrm{minimize}J(\mathbf{d})=W_1{f}_1(\mathbf{d})+W_2{f}_2(\mathbf{d})\\ & \mathrm{subject}\mathrm{to}\end{array}\\ \begin{array}{rlr}& \mathrm{state}eq.\text{eqref}eq:pde,& (\mathrm{PDE}\mathrm{constraint})\\ & {\mathbf{d}}_{\min }\le \mathbf{d}\le {\mathbf{d}}_{\max }.& (\mathrm{Parameter}\mathrm{constraint})\end{array}\end{array}$$

$${\mathcal{C}}_{p,q}^{\mathrm{s}}$$

$${\mathcal{C}}_{p,q}^{\mathrm{wd}}$$

$${\mathcal{M}}_{p,q}(\mathbf{x};{\mathbf{d}}_{p,q})=\underset{t\in {\mathcal{T}}_{p,q}}{\max }\{u\left(\mathbf{x},t;{\mathbf{d}}_{p,q}\right)\}.$$

$$g_1({\mathbf{d}}_{p,q})$$

$$g_2({\mathbf{d}}_{p,q})$$

$${\beta }_{p,q}({\mathbf{x}}^{\mathrm{s}})={\chi }_{(\gamma ({\mathbf{x}}_p^{\mathrm{i}}),\gamma ({\mathbf{x}}_{p+r}^{\mathrm{i}})]}+\sum _{k=p+r}^{q-1}{\left(\frac{\gamma ({\mathbf{x}}_q^{\mathrm{f}})-\gamma ({\mathbf{x}}_k^{\mathrm{i}})}{\gamma ({\mathbf{x}}_q^{\mathrm{f}})-\gamma ({\mathbf{x}}_p^{\mathrm{k}})}\right)}^2{\chi }_{(\gamma ({\mathbf{x}}_k^{\mathrm{i}}),\gamma ({\mathbf{x}}_k^{\mathrm{f}})]}$$

$$\begin{array}{r}\begin{array}{rl}& \mathrm{minimize}{J}_{p,q}({\mathbf{d}}_{p,q})=W_1{g}_1({\mathbf{d}}_{p,q})+W_2{g}_2({\mathbf{d}}_{p,q})\\ & \mathrm{subject}\mathrm{to}\end{array}\\ \begin{array}{rlr}& \mathrm{state}eq.\text{eqref}eq:pde,& (\mathrm{PDE}\mathrm{constraint})\\ & {\mathbf{d}}_{{p,q}_{\min }}\le {\mathbf{d}}_{p,q}\le {\mathbf{d}}_{{p,q}_{\max }}.& (\mathrm{Parameter}\mathrm{constraint})\end{array}\end{array}$$

$$\begin{array}{r}\sigma ({\mathbf{x}}^{\mathrm{s}})=\sum _{k=1}^N{F}^{\sigma }({\mathbf{x}}_k^{\mathrm{i}};{\mathbf{\Lambda }}^{\sigma }){\chi }_{(\gamma ({\mathbf{x}}_k^{\mathrm{i}}),\gamma ({\mathbf{x}}_k^{\mathrm{f}})]},\\ v({\mathbf{x}}^{\mathrm{s}})=\sum _{k=1}^N{F}^v({\mathbf{x}}_k^{\mathrm{i}};{\mathbf{\Lambda }}^v){\chi }_{(\gamma ({\mathbf{x}}_k^{\mathrm{i}}),\gamma ({\mathbf{x}}_k^{\mathrm{f}})]},\end{array}$$

$${\sigma }_k$$

$$v_k$$

$$\begin{rcases*}F^{\sigma}(\bm{x}^{\mathrm{s}};\bm{\Lambda}^{\sigma})=F_{l}^{\sigma}(\bm{x}^{\mathrm{s}};\bm{\Lambda}_{l}^{\sigma})\\ F^{v}(\bm{x}^{\mathrm{s}};\bm{\Lambda}^{v})=F_{l}^{v}(\bm{x}^{\mathrm{s}};\bm{\Lambda}_{l}^{v})\\ \end{rcases*}\,\,\mathrm{if}\,\,\bm{x}^{\mathrm{s}}\,\,\mathrm{on\,\,hatch\,\,line}\,\,l$$

$${\mathcal{M}}_{p,q}(\mathbf{x};{\mathbf{d}}_{p,q})\approx \frac{1}{K}\log \left({\int }_{{\mathcal{T}}_{p,q}}\exp (K\cdot u(\mathbf{x},s;{\mathbf{d}}_{p,q}))\mathrm{d}s\right),$$

$${\mathcal{C}}^{\mathrm{wd}}=\{(x^{\mathrm{s}}(t),y^{\mathrm{s}}(t)+w,-d):t\in \mathcal{T}\}.$$

Full text

1. Introduction

Powder bed fusion (PBF) is a type of additive manufacturing (AM) where metal powder is melted by a laser or electron beam in a layer-wise fashion to enable the production of geometrically complex parts [7]. AM undergoes continuous progress towards a technology that is robust and efficient, but there are still issues when it comes to quality and repeatability.

It is important for completed parts to meet the mechanical requirements and quality standards specified by the applications for which they are manufactured. The qualities of a completed part, such as tensile strength and surface roughness, depend on the melting process, which in turn is governed by several dozen material and process parameters such as preheating temperature, powder packing ratio, and beam speed, among other. Correlations between process parameters, process signatures (such as melt pool size and temperature), and product qualities are presented in [17] and references therein. Due to these correlations, the design of process settings requires critical attention. However, this involves extensive and costly experimental work and the resulting melting schemes implemented in machines rely on an excessive amount of parameters and functions in order to account for the dynamics of the melting process. This approach makes it difficult to optimize PBF and limits both the number of applications and the number of materials available for manufacturing.

As of late, computation based techniques are used to improve PBF related process control, i.e., to tackle the question of how to select different process parameters in order to build parts with desired properties [17, 13]. A common approach is design of (computational) experiments (DoE) [18, 25], which is an exploratory tool used to identify parameters that influence the qualities of the completed part. In [15], a DoE with a finite element model determines that the beam power and beam speed are the two process controllable parameters that have the largest impact on peak temperature during a single track melt. In [10], a DoE based on a simple thermal model [3] is used to determine optimal process parameters for building high density parts. Also this study shows that the beam power and beam speed have the largest impact on melt pool width and depth. Empirical modelling techniques such as artificial neural networks are also used to determine the optimal selection of process parameters, see, e.g., [20] as well as [5] and references therein.

While DoE is useful for developing a deeper understanding of the melting process, it can be a tedious affair to use them for optimization. An optimal control approach might be better suited for that purpose, as it starts off with preset, desired melt characteristics and seeks corresponding optimal process parameters via a mathematical optimization problem.

Attempts of optimizing melt pool characteristics have been made based on the two-phase Stefan problem, where the free boundary between solid and liquid is understood via the Stefan condition. In the context of PBF, the free boundary characterizes the size and shape of the melt pool. In [8], the two-phase Stefan problem in a container is considered, and the temperature on the container boundary is optimized with respect to a desired transient evolution of the free boundary. Optimization problems based on quasi-steady state formulations of the two-phase Stefan problem are solved in [24, 2]. Here the desired free boundary between liquid and melt is prescribed and the goal involves tracking the melt temperature on the prescribed free boundary.

This paper is part of an effort that aims to reduce the number of parameters needed in the design of melt schemes. We present a simulation based framework that seeks to facilitate process optimization and material development, while keeping computational costs at a minimum. Thus, computational efficiency is prioritized as we trade a certain level of detail for a fast optimization scheme that can be applied to the melting of large domains. The scheme is efficient for three reasons:

• $\triangleright$

The continuum thermal model for describing heat conduction includes an analytic solution that allows for fast and pointwise computation of temperatures during melting [16, 4]. The model assumes that the beam parameters are piecewise constant in time. We remark that beam parameters in actual AM machines are often set to vary in such a discontinuous way. The model does not include phase change, and does in particular not capture the solid-liquid interface. Instead, we use the term melt pool to simply denote the region where the temperature is larger than the melt temperature.

• $\triangleright$

The formulation of the optimization problem involves a severe model order reduction. Rather than striving for some desired temperature distribution that is difficult to express, the goal consists of tracking preset reference maximum temperatures on lines that run along the beam path.

• $\triangleright$

The resulting optimization problem is solved by a greedy algorithm that divides it into several small problems that are easier to solve. These sub-problems are solved consecutively as the beam traverses the powder bed.

Together, this framework comprises an optimization scheme for process control as suggested in [22].

As noted earlier, DoE suggests that the beam parameters have the largest impact on the temperature distribution during melting and, ultimately, on the quality of the completed part. For this reason, we choose the beam power, spot size and speed as control variables. Our scheme allows for time-dependent beam parameters, which increases their ability to control the melting process. An ability to optimize these beam parameters is useful not only for validation, but also in order to speed up the development of process settings for new (and old) materials. Since optimization is an iterative process, the solver of the forward problem needs to be highly efficient. The analytic solution provides such efficiency.

The remainder of this paper is organized as follows. Section 2 describes the thermal model and its corresponding analytic solution. The optimization problem is formulated in Section 3. Section 4 and Section 5 propose two greedy algorithms for solving said problem. In Section 6 we apply our optimization scheme in numerical examples. Here we also detail how the combination of the two greedy algorithms can aid in the development of so called beam parameter functions. Additional comments and concluding remarks are given in Section 7.

2. Thermal model

Consider the heat equation on the lower half space $\Omega=\mathbb{R}^{2}\times\mathbb{R}_{-}$ during a time span $\mathcal{T}=(0,T]$. Let $\Gamma$ denote the surface boundary $z=0$ and let $u_{\mathrm{init}}$ denote the constant initial temperature. The beam travels on the surface $\Gamma$ along a preset, piecewise linear path

[TABLE]

where $\bm{x}^{\mathrm{s}}(t)=(x^{\mathrm{s}}(t),y^{\mathrm{s}}(t),0)$ is the position of the center of the beam at time $t$. The heat flux $\Phi$ due to a scanning electron beam is modeled as a Gaussian function

[TABLE]

The three beam parameters are the absorbed beam power $P(t)$, the beam spot size $\sigma(t)$, and the beam speed $v(t)=|\bm{v}(t)|$. Here $\bm{v}(t)=(v_{x}(t),v_{y}(t),0)=(v(t)\cos{\theta(t)},v(t)\sin{\theta(t)},0)$, where $\theta(t)$ is the angle between the positive $x$-axis and the direction of the path. This angle is known for any $t$ since the beam path $\mathcal{C}^{\mathrm{s}}$ is preset and it follows that $v$ uniquely defines the vector $(v_{x},v_{y},0)$. The position of the beam $\bm{x}^{\mathrm{s}}(t)$ depends on the speed with which the beam has traveled the path $\mathcal{C}^{\mathrm{s}}$ up to time $t$. The beam parameters are often set to vary in a piecewise constant fashion in AM machines. The following definition makes the concept of piecewise constant beam parameters more precise.

Definition 1**.**

Given times $0=t_{0}<t_{1}<\ldots<t_{N}=T$,

[TABLE]

is a partition of $\mathcal{T}$ consisting of $N$ segments. Index $n$ indicates the $n^{th}$ segment in the partition, and a segment in turn is a collection of the following data:

• –

$t_{n}^{\mathrm{i}},\,\,t_{n}^{\mathrm{f}}$; an initial time and final time, respectively,

• –

$(P_{n},\sigma_{n},v_{n})$; a triplet of power, spot size, and speed such that

[TABLE]

• –

$\ell_{n}^{\mathrm{s}}$; a line traversed by the beam between times $t_{n}^{\mathrm{i}}$ and $t_{n}^{\mathrm{f}}$,

• –

$\bm{x}_{n}^{\mathrm{i}},\,\,\bm{x}_{n}^{\mathrm{f}}$; the initial position and final position of $\ell_{n}^{\mathrm{s}}$, respectively.

• –

$\theta_{n}=\mathrm{tan}^{-1}\left(\frac{y_{n}^{f}-y_{n}^{i}}{x_{n}^{f}-x_{n}^{i}}\right)$; the angle between the positive $x$-axis and the direction of the path.

With $\hat{z}$ the outward unit normal of $\Omega$, the heat transfer problem can be written as

[TABLE]

where $\rho$, $c_{p}$, $\lambda$ denote density, heat capacity, and thermal conductivity, respectively. These material parameters are assumed to be constant. This gives us the thermal diffusivity $\kappa=\lambda/\rho c_{p}$. Problem (1) has an analytic solution [16, 4]. We refer to these sources for a detailed derivation of this solution and merely outline it here.

Proposition 1**.**

Given a partition as in Definition 1, the solution of problem (1) can be written as

[TABLE]

where $u_{n,k}^{\mathrm{I}}$ is the temperature due to the earlier scanning of segment $k<n$ and $u_{n}^{\Phi}$ is the temperature due to the current scanning of segment $n$.

The analytic expressions of $u_{n,k}^{\mathrm{I}}$ and $u_{n}^{\Phi}$ are derived in [4], wherein it also described how the solution can be efficiently computed.

The power is largely determined by the beam current, and this current can not be adjusted at a fast rate. Therefore, $P$ is set to be constant for all $t\in\mathcal{T}$ for the remainder of this paper. It should be noted, however, that if one would be interested in optimizing all three beam parameters, the following extends to the case of non-constant power as well.

A remark on our thermal model is in order, as it does not include cooling, the latent heat of fusion nor a description of the solid-liquid intersection between powder and melt pool, which makes the notion of a melt pool quite fuzzy. Here we use the term melt pool to simply denote the volume where the temperature is larger than the melt temperature. Hence we use the isothermal $\{\bm{x}(t)\mathrel{\mathop{\ordinarycolon}}u(\bm{x},t)=u_{\mathrm{melt}}\}$, where $u_{\mathrm{melt}}$ is the melt temperature of the powder, to represent the solid-liquid interface. Furthermore, since the model is a continuum model, it breaks down on the mesoscale where we see phenomena such as balling, inter-capillary effects, Plateau–Rayleigh instabilities, thermal expansion, among other [11, 12]. Despite these restrictions, it is anticipated that effective parameters, tuned via comparisons with experiments, can be used to make the thermal model reliable enough for control and optimization. The optimization problem described below also aligns with the overarching aim to reduce the number of parameters needed for process control.

3. Formulation of the optimization problem

The goal is to optimize the melting process with respect to the beam spot size and beam speed. Since these beam parameters are defined in a piecewise constant fashion according to Definition 1, we can write

[TABLE]

where $\chi$ is the indicator function. Since we aim to optimize the speed, either the times $t_{n}^{\mathrm{i}}$, $t_{n}^{\mathrm{f}}$ or the positions $\bm{x}_{n}^{\mathrm{i}}$, $\bm{x}_{n}^{\mathrm{f}}$, $n=1,2,\ldots,N$, in Definition 1 will have to vary during optimization. Due to the beam path being preset and reasons that will become clear in the next section, it is better to fix the positions. Therefore it is more appropriate to express $\sigma$ and $v$ as space dependent functions instead. To this end, introduce the scanning distance

[TABLE]

Then we have the following beam parameter functions:

[TABLE]

From (2), a decision vector can immediately be extracted as $\bm{d}=(\bm{\sigma},\bm{v})=\{(\sigma_{k},v_{k})\}_{k=1}^{N}=(\sigma_{1},v_{1},\sigma_{2},v_{2},\ldots,\sigma_{N},v_{N})$. The variables in the decision vector are bounded due to practical limitations, $\bm{d}_{\mathrm{min}}\leq\bm{d}\leq\bm{d}_{\mathrm{max}}$.

The control of the melting process is a multiobjective optimization problem due to the many correlations between process parameters, process signatures, and product qualities. The qualities of the final part are strongly dependent on the temperatures obtained during the melting process [19]. The characteristics of the melt pool are important process signatures. If the melt pool is too small relative to the line offset (i.e., the distance between two hatch lines) and layer depth, powder might be left unmelted between hatch lines or between layers, causing discontinuities and porosity. Furthermore, high surface temperatures might result in too much evaporation and subsequent recoil pressure, which can result in undesired material transport such as ejection of molten materials that later cause defects [6] or formations of small ridges that prohibit the deposition of new powder layers and thus cause the manufacturing process to terminate [23]. Qualitatively, therefore, the choice of cost functional can be motivated by the desire to

1. (1)

maintain a uniform and appropriately sized melt pool during melting, and 2. (2)

avoid too high surface temperatures.

3.1. A reductive approach

Consider a beam scanning along a path $\mathcal{C}^{\mathrm{s}}$. Denote by $\omega$ a supposed desired melt pool

[TABLE]

It is difficult to express $\omega(t)$ explicitly. Instead, we consider the final solidified volume. With our purely thermal model, the powder-solid interface of this volume is dependent on the maximum temperature and would be easier to explicitly define than the melt pool $\omega(t)$. However, even further reductions can be made by isolating particular curves on this powder-solid interface. More precisely, we introduce a secondary path $\mathcal{C}^{\mathrm{wd}}$ chosen such that it lies on a desired powder-solid interface. The secondary path is related to the beam path by some function $\mathfrak{F}\mathrel{\mathop{\ordinarycolon}}\mathcal{C}^{\mathrm{s}}\rightarrow\mathcal{C}^{\mathrm{wd}}$ and we write

[TABLE]

and let $\bm{x}^{\mathrm{wd}}=\mathfrak{F}(\bm{x}^{\mathrm{s}})$. For example, if the beam path consists of one segment, then a simple example of a secondary path is

[TABLE]

The idea is that the secondary path relates to $\mathcal{C}^{\mathrm{s}}$ via a width $w$ and a depth $d$. With this, the description of the optimal melting reduces to two paths; $\mathcal{C}^{\mathrm{s}}$ (preset) and $\mathcal{C}^{\mathrm{wd}}$ (chosen with respect to $\mathcal{C}^{\mathrm{s}}$). This reductive approach is illustrated in Figure 1.

As we shall see in the following section, the steps taken above allow us to formulate a simple optimization problem that is efficient in the sense that we, instead of tracking some desired transient melt pool $\omega(t)$ in a volume, only track two scalar values $u_{\mathrm{melt}}$ and $u_{\mathrm{surf}}$ on paths.

3.2. Mathematical formulation

Following the reduction in Section 3.1, we are now interested in maximum temperatures on paths running along the beam path since these temperatures determine the size of the subsequent solidified volume. Before we formulate the problem, we need the following.

Definition 2**.**

$\mathrm{(Hatch\,\,line).}$* Given a partition as in Definition 1, two segments $k$ and $k+1$, $1\leq k<N-1$, are connected if $\bm{x}_{k}^{\mathrm{f}}=\bm{x}_{k+1}^{\mathrm{i}}$. A sequence $\{i\}_{i=k}^{K}$ of connected segments form a hatch line if $\theta_{n}=\theta_{m}$ $\forall n,m\in[k,K]$ and $\theta_{k-1}\neq\theta_{k}$ and $\theta_{K}\neq\theta_{K+1}$.*

The total number of hatch lines $M$ satisfies $1\leq M\leq N$.

Define also the maximum temperature field

[TABLE]

We want $\mathcal{M}(\bm{x};\bm{d})=u_{\mathrm{melt}}$ for all $\bm{x}$ on $\mathcal{C}_{i}^{\mathrm{wd}}$ in order to ensure a uniform and thorough melting. Similarly, $\mathcal{M}(\bm{x};\bm{d})$ should not exceed some maximum allowed temperature $u_{\mathrm{surf}}$ on $\mathcal{C}^{\mathrm{s}}$.

The resulting objective vector becomes

[TABLE]

where

[TABLE]

Here

[TABLE]

is a weight that excludes intervals from the cost functional if they are very close to the start point or end point of a hatch line. This type of weight is inserted because the total heat supplied to a region near a start point, for instance, is comparatively small since the beam only moves away from it rather than passing it. As a consequence, it can be difficult to reach the reference temperatures in these regions and if included, they deteriorate the overall performance of the optimizer. Therefore, it is better to ignore these intervals in the goal functional and instead let them be covered by the contouring stage, in which the beam scans along the boundary of the shape being melted. The contouring stage also improves the surface finish of the part [21]. Effectively, this choice of $\alpha$ simply means that the domains of integration in (4) become slightly smaller.

The functionals (4) are of tracking type where we use the $L_{2}$ norm to minimize the distances between the actual maximum temperature and the desired maximum temperatures. The reason for tracking the surface temperature rather than penalizing only too high temperatures is that it helps restricting the shape of the melt pool. Without this restriction, one could potentially end up with an extremely wide and shallow type of melting, for instance.

The analytic solution presented in Section 2 allows for pointwise computations of temperatures, and that is why we can easily compute temperatures on lines, which saves a large amount of computation time compared to doing so on surfaces or in volumes. This remark highlights a big motivation behind the reduction carried out in Section 3.1.

The objective vector (3) needs to be translated into a scalar valued cost functional in order to use standard nonlinear programming solvers. We use a scalarization known as the weighting method. In this method the weighted sum of the objectives is minimized. We introduce weights $W_{i}\geq 0,\,\,i=1,2$. The scalarized optimization problem becomes

[TABLE]

The choice of weights should represent the relative importance of the objectives; important objectives are weighted more heavily.

4. A first greedy algorithm for solving the scalarized optimization problem

In order to speed up the optimization, we propose a method that makes use of the fact that the melt pool, and hence maximum temperatures, are localized to the beam. The proceeding involves a division of $\mathcal{T}$ into subintervals on the form $\bigcup_{i=p}^{q}(t_{i}^{\mathrm{i}},t_{i}^{\mathrm{f}}]$. Local optimization problems are solved on these subintervals and optimal parameter pairs $(\sigma_{p},v_{p})$ are frozen sequentially. When given parameter pair(s) has been frozen, the local problem is translated in time (and space) and the initial condition is updated. As such, this greedy type of algorithm divides the optimization problem (6) into several smaller optimization problems that are faster to solve.

The goal functional in (6) involves maximum temperatures over time near the beam path, which is a property that is local to the beam itself. Given a point $P$ on, say, $\mathcal{C}^{\mathrm{s}}$, it is known that $P$ will obtain its largest temperature during a time window when the beam, and the melt pool it generates, passes $P$. Therefore, it is the values of the beam parameters during this particular time window that has the highest influence on maximum temperature at $P$. The reasoning is similar for a point on $\mathcal{C}^{\mathrm{wd}}$, the only difference being that it takes slightly longer to reach the maximum temperature on $\mathcal{C}^{\mathrm{wd}}$ since heat diffusion is not an instantaneous process. Therefore, we decide to split the optimization problem (6) into multiple subproblems that are solved sequentially in time while parameter pairs are frozen as we go along.

In order to formalize the method, we make the following definition.

Definition 3**.**

$\mathrm{(Window).}$* Given a partition as in Definition 1, a time window $\mathcal{T}_{p,q}=\cup_{i=p}^{q}(t_{i}^{\mathrm{i}},t_{i}^{\mathrm{f}}]$ is defined as a set of adjacent segments in $\mathcal{T}$. The size of $\mathcal{T}_{p,q}$ is the number of segments that constitutes it. The local beam path and local secondary path corresponding to $\mathcal{T}_{p,q}$ are given by*

[TABLE]

Hence the beam traverses the path $\mathcal{C}_{p,q}^{\mathrm{s}}$ during time $\mathcal{T}_{p,q}$.

Define a local decision vector $\bm{d}_{p,q}=\{(\sigma_{k},v_{k})\}_{k=p}^{q}$ and a local maximum temperature field

[TABLE]

Similarly, we define the local objectives

[TABLE]

Here $\beta_{p,q}$ is used to prioritize minimization of the errors over the earlier segments in the window. This weight accentuates the error on the segment(s) that is about to become frozen and it plays a crucial role; since the greedy algorithm never returns to a segment once it has been frozen it is important that the solver prioritizes this segment. Here we let $\beta_{p,q}$ be piecewise constant over the segments and determined by a function that decreases quadratically along $\mathcal{C}_{p,q}^{\mathrm{s}}$. See Figure 2. Formally, we have

[TABLE]

on $\mathcal{C}_{p,q}^{\mathrm{s}}$. We define $\beta_{p,q}(\bm{x}^{\mathrm{wd}})$ in a similar fashion. The choice of a quadratic underlying function is based on tests that investigate how the weight affects the optimization results.

Now, by employing the same scalarization as for the global problem (6), the resulting scalarized subproblem becomes

[TABLE]

We can now formulate the greedy algorithm. This is done in Algorithm 1, and some complementary comments are given below. An illustration of the main idea is given in Figure 3.

1. $\mathscr{L}4$:

Future parameter pairs are updated as well, because the values found in the current window are likely a better guess than the initial one. 2. $\mathscr{L}5$:

The segments corresponding to the frozen parameter pairs are removed from the window. In the implementation, certain checks can be made to determine whether freezing should take place or not (as in, the amount of pairs to freeze). We leave out the details. 3. $\mathscr{L}7$-$9$:

The window size is updated in preparation for the next iteration. The min operator is used to handle the ending when $q=N$. Note also that the parameters $q$ and $r$ may depend on $p$ (i.e., on the location of the window). For instance, in a region where the beam path is complex or where the lengths of the segments are small, we might require a large window size $q-p+1$, and hence a large $q$. Furthermore, the size of the melt pool should be taken into account when choosing the value of $q$.

The presented greedy algorithm can be utilized as a standalone tool for process optimization. If the beam path consists of $N$ segments, the total number of parameters to optimize becomes $2N$. Now, depending on the design of the layer being melted, the value of $N$ may be very high. In the next section we present a second version of the greedy algorithm that has the ability to significantly lower the size of the decision vector.

5. A second greedy algorithm based on fitting beam parameter functions

The greedy algorithm of this section expands on an example in [4]. As we shall see, it is similar to Algorithm 1 in most regards, but it is based on educated guesses of how the beam parameters should behave.

Recall from (2) the expression for piecewise constant beam parameters. The subsequent optimization makes no assumption on the behavior of the beam parameters along the beam path. An alternative approach is to do curve fitting of prespecified beam parameter functions. To this end, we write

[TABLE]

where $\bm{\Lambda}=(\bm{\Lambda}^{\sigma},\bm{\Lambda}^{v})$, the coefficients in our beam parameter functions, become our new decision vector that we want to optimize. Given a decision vector, the beam parameters are then evaluated as (see (2))

[TABLE]

These parameter functions are defined with respect to the hatch lines as

[TABLE]

The reason for splitting $F^{\sigma}$ and $F^{v}$ between hatch lines is that a new hatch line often requires a rapid jump in beam parameter values.

It follows from Definition 2 that hatch lines can simply be seen as an intermediate level between the segments and the beam path. Let $\mathpzc{S}\mathrel{\mathop{\ordinarycolon}}\{1,\ldots,M\}\rightarrow\{1,\ldots,N\}$ be an injective function that, given a hatch line $m$, returns the first segment in $m$. Let $\mathpzc{L}\mathrel{\mathop{\ordinarycolon}}\{1,\ldots,N\}\rightarrow\{1,\ldots,M\}$ be a surjective function that, given segment $n$, returns the hatch line that contains it. With these two functions it is possible to seamlessly work with both hatch lines and segments. For instance, the local beam path that consists of hatch lines $1$ to $3$ is $\mathcal{C}_{\mathpzc{S}(1),\,\mathpzc{S}(4)-1}^{\mathrm{s}}$.

In terms of implementation, the second greedy algorithm is in many ways similar to the first greedy algorithm from the previous Section 4. They both rely on Definition 1 and piecewise constant beam parameters and they both solve subproblems of the form (7). What separates them is the content of the decision vector $\bm{d}$ as illustrated in Figure 4. In essence, the first algorithm optimizes the beam parameters segment wise while the second algorithm optimizes the beam parameters line wise. The second greedy algorithm is detailed in Algorithm 2.

6. Numerical examples and discussion

We apply the greedy algorithm on a couple of single layer problems. The scalarized subproblems (7) are solved with the L-BFGS-B optimization algorithm [26, 1] provided by the optimization package of SciPy [9]. Given some iterate $\hat{\bm{d}}_{p,q}$, scipy.optimize approximates the gradient of $J_{p,q}(\hat{\bm{d}}_{p,1})$ using a 2-point finite difference estimation. Then L-BFGS-B, which is a quasi-Newton method, approximates the Hessian that enters in the local quadratic approximation of $J_{p,q}(\hat{\bm{d}}_{p,q})$. The solver options are selected to fit the scale of the problems considered here.

It is worth noting that the objective $J_{p,q}$ is not differentiable. At an initial stage, not only L-BFGS-B but also some gradient free methods were tested, and L-BFGS-B performed the best out of all solvers in those trials. It is not remarkable that L-BFGS-B performs well on nonsmooth problems as well (although we can not expect the same convergence as for a smooth optimization problem) [14]. On a related note, in the implementation we relax the scalarized subproblem (7) somewhat by replacing the (local) maximum temperature field $\mathcal{M}_{p,q}(\bm{x};\bm{d}_{p,q})$ with an approximation,

[TABLE]

for an appropriate scalar $K$, which improves performance slightly. Finally, while the choice of starting point/initial decision vector can have a large impact on the performance of the optimizer, efforts related to this choice are not the main focus here and so disregarded.

It is important to emphasize that for practical use of the optimization scheme, the material data that enter the thermal model need to be fit with respect to experiments or a more detailed model. The determination of effective parameters is crucial since the model is simple and based on several assumptions. In the following examples we use material parameters that represent Ti-$6$Al-$4$V.

In the current implementation of the beam scanning, there are no pauses between any segments during melting. For instance, a jump from one hatch line to the next is instantaneous. However, adding delay time for jumps is straightforward.

6.1. Example 1: segment wise optimization on snake pattern

We illustrate how the optimization scheme resolves certain heating related issues. The beam path is shown in Figure 5. It consists of $50$ segments, each with length $0.5\,\mathrm{mm}$. The line offset during hatching is $l_{\mathrm{off}}=200\,\mu\mathrm{m}$. The hatching is performed in a snake-like manner in the upward $y$-direction. The weight $\alpha(\bm{x})$ (see (5)) is zero on the first $0.4\,\mathrm{mm}$ and last $0.4\,\mathrm{mm}$ of a hatch line. While this beam path amounts to a simple rectangular shape, it still allows for several types of investigations.

The secondary path is set to

[TABLE]

To avoid porosity, $w$ and $d$ need be chosen such that unmelted gaps between lines are avoided. Here we set $w=l_{\mathrm{off}}/2=100\,\mu\mathrm{m}$ and $d=50\,\mu\mathrm{m}$. A complete list of parameter values is given in Table 1.

Figure 6 shows the solution of (6) as obtained by the greedy Algorithm 1. We see rapid variations in the optimized beam parameters, in particular close to the turning points where they seek to neutralize the concentrated influx of heat that occur in those regions. Problem (1) is then solved for the optimized beam parameters and the resulting maximum temperature in various slices of the domain is shown in Figures 7, 8, and 9. These figures include the initial maximum temperature for comparison.

The results indicate that despite the reductions leading up to its formulation, the greedy algorithm is able to control the heat generated during melting to a rather large degree.

One concern with the greedy algorithm is that it carries with it several uncertainties. Many trials are required to find proper values for the parameters that make up the scheme, such as $\alpha$, the window size and segment lengths, since they depend on the beam path and thermal diffusivity.

6.2. Example 2: segment wise optimization on nonparallel pattern

We melt the first quadrant of an annulus. The annulus has an inner radius $r_{i}=1\,\mathrm{mm}$ and outer radius $r_{o}=5\,\mathrm{mm}$. The beam path is shown in Figure 10.

It consists of 19 lines of length $4\,\mathrm{mm}$. Two adjacent lines differ by an angle of $5{}^{\circ}$. Consequently, the distance between them at $r_{i}$ and $r_{o}$ are about $87\,\mu\mathrm{m}$ and $436\,\mu\mathrm{m}$, respectively. Each line is divided into $8$ segments of equal length $0.5\,\mathrm{mm}$. The hatching is performed in the counter-clockwise direction. The thin gray lines indicate the secondary path $\mathcal{C}^{\mathrm{wd}}$. Once again, the weight $\alpha(\bm{x})$ is zero on the first $0.4\,\mathrm{mm}$ and last $0.4\,\mathrm{mm}$ of a hatch line. The remaining parameters are identical to the ones used in the previous example and are listed in Table 1.

Only the blue part of the beam path is considered during optimization. More precisely, we apply the greedy algorithm on the first $5$ lines only. The results from this optimization is then extended by letting the beam parameters on lines $6$-$19$ equal the optimal beam parameters on line $5$.

Figure 11 shows the solution of the optimization problem. The speed increases along each hatch line since the width between a hatch line and the corresponding secondary path decreases as the hatch line approaches $r=1\,\mathrm{mm}$ (see Figure 10).

After having extended this solution to the entire beam path, problem (1) is solved for the optimized beam parameters and the resulting maximum temperature is shown in Figures 12 and 13. These figures include the initial maximum temperatures for comparison.

The plots in Figure 13 confirm that the extension of the solution onto remaining lines $6$-$19$ works well. While the temperature on $\mathcal{C}^{\mathrm{wd}}$ increases near the inner radius of the annulus as the scanning progresses, this increase is small and does not justify applying the greedy algorithm on the entire beam path, $19$ lines, rather than just $5$ lines. This is just a small example of how the results from the greedy algorithm on a very small section can be utilized on larger sections of the build area. In general, this procedure offers an efficient method for improving process control: first examine and optimize typical problematic melting scenarios, then combine the results and extend them to the remainder of the layer.

6.3. Example 3: line wise optimization on snake pattern

We now use the second greedy algorithm to solve the problem introduced in the first example, Section 6.1, and compare the results with the results obtained by the first greedy algorithm. Motivated by the solution obtained by the first greedy algorithm, see Figure 6, we make the following ansatz. Let

[TABLE]

and

[TABLE]

for $l=1,\ldots,M$. Hence we associate $8$ coefficients with each hatch line. From the beam path in Figure 5, we have $M=5$. We solve optimization problem (6) according to Algorithm 2. The window always consist of $1$ hatch line, i.e., $10$ segments. The results are shown in Figure 14. The solution obtained by the first greedy algorithm is included for comparison. The results are similar. The corresponding optimal objectives $J^{\mathrm{opt}}$ are similar as well, as $J^{\mathrm{opt}}=124.18$ with the second greedy algorithm and $J^{\mathrm{opt}}=123.03$ with the first greedy algorithm ($J^{\mathrm{init}}=1655.21$).

We end this section with a comparison between the first greedy algorithm and second greedy algorithm. The first algorithm from Section 4 optimizes the beam parameters segment wise. It can be applied to general beam paths and there are no restrictions on the window. In particular it and can be used to look at specific problematic areas of the layer being melted in order to get an understanding of how the beam parameter functions should behave in those areas.

The second greedy algorithm from Section 5 optimizes the beam parameters line wise. This makes the second algorithm preferable in practical problems because it significantly decreases the dimension of the decision space; the number of lines $M$ is much lower than the number of segments $N$. In the above example, the hatch lines are fairly short, but we still get $40=8M<2N=100$ when comparing dimensions of the two decision spaces. For more realistic problems the number of segments could be orders of magnitude larger, making the second greedy algorithm more attractive.

The potential drawback of the second algorithm is that it might be difficult to make the initial ansatz for the beam parameter functions since they depend on the beam path (and secondary path). However, this is where the first algorithm can be utilized to give an initial estimate that shows the behavior of the desired beam parameters. This approach is what enabled us to choose the beam parameter functions in (9). For the future, we imagine the development of a database of parameter functions that have been generated for different melting scenarios and that can be shared and used in other, more detailed models.

7. Conclusions

We have formulated an optimization scheme for controlling the heat conduction during the melting process in powder-bed-based additive manufacturing. The scheme is efficient because it exploits that the melt pool is local to the beam and shows good capabilities despite several simplifications. The current choice of objectives prioritizes speed since it only requires temperature evaluations on lines rather than in entire volumes. The scheme should be useful for studying problematic areas of the melting process where particular care needs to be put into the choice of beam parameters.

The optimization scheme relies on a greedy algorithm. Two versions of a greedy algorithm have been presented and applied in this paper. The first one carries out optimization segment wise, which makes it applicable to many types of beam paths. The second one carries out optimization line wise, which can significantly reduce the dimension of the decision space. While the two algorithms are similar, they serve different purposes and we have detailed how they can be combined to improve process control.

The examples considered in this paper are fairly small. For more realistic problems the number of segments could be orders of magnitude larger. By design, the greedy algorithm becomes more attractive as the total amount of segments in the beam path increases; the division of the global problem (6) into subproblems (7) becomes more beneficial, relatively speaking, as $N$ (and $M$) increases. Furthermore, the examples are purely numerical. It is currently difficult to compare the results in Section 6 with experiments because existing machines lack the functionality required to match our numerical results. Because of this, a crucial next step is to implement the necessary code in the machine such that experimental validation becomes possible. Experiments are also needed for the generation of effective parameters; the analytic model is very simple and since it contains few parameters, they need to be carefully fit against experiments.

The optimization method also needs to be complemented with different types of testing. It requires effective reference temperatures $u_{\mathrm{melt}}$, $u_{\mathrm{surf}}$. Furthermore, the weights and the secondary path $\mathcal{C}^{\mathrm{wd}}$ need to be carefully chosen. Work related to this kind of testing has not been detailed here.

Acknowledgements

This work was supported by the Swedish Foundation for Strategic Research under the contract “Industrial PhD 2015 – ID15-0058”.