# Optimization of Engineering Processes Including Heating in Time-Dependent Domains

Alfred Schmidt, Eberhard Bänsch, Mischa Jahn, Andreas Luttmann, Carsten

Niebuhr, Jost Vehmeyer

► To cite this version:

Alfred Schmidt, Eberhard Bänsch, Mischa Jahn, Andreas Luttmann, Carsten Niebuhr, et al. Optimization of Engineering Processes Including Heating in Time-Dependent Domains. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.452-461, 10.1007/978-3-319-55795-3_43. hal-01626905

HAL Id: hal-01626905

https://inria.hal.science/hal-01626905

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimization of Engineering Processes Including Heating in Time-Dependent Domains

Alfred Schmidt¹*, Eberhard Bänsch², Mischa Jahn¹, Andreas Luttmann¹, Carsten Niebuhr¹, and Jost Vehmeyer¹

¹ Zentrum für Technomathematik, Universität Bremen, Germany
http://www.math.uni-bremen.de/zetem

² Department Mathematik, Friedrich-Alexander-Universität Erlangen-Nürnberg, Germany http://www.mso.math.fau.de/applied-mathematics-3{schmidt,mischa,andreasl,niebuhr,vehmeyer}@math.uni-bremen.de,baensch@math.fau.de

**Abstract.** We present two models for engineering processes, where thermal effects and time-dependent domains play an important role. Typically, the parabolic heat equation is coupled with other equations. Challenges for the optimization of such systems are presented.

The first model describes a milling process, where material is removed and heat is produced by the cutting, leading to thermomechanical distortion. Goal is the minimization of these distortions.

The second model describes the melting and solidification of metal, where the geometry is a result of free-surface flow of the liquid and the microstructure of the re-solidified material is important for the quality of the produced preform.

**Keywords:** Optimization with PDEs, time-dependent domain, heat equation, thermoelasticity, free surface flow

## 1 Introduction

The optimization of industrial engineering processes often lead to the treatment of coupled, nonlinear systems of PDEs. Here we want to investigate applications where an important part of the nonlinearities is created by time-dependent domains, which are not prescribed but who are part of the solution itself. The treatment of such models in the context of PDE constraints in optimal control problems typically generates additional challenges in both, the solution of the forward problems and the treatment and storage of adjoint solutions.

In the following, we present two models for engineering processes, where heating of metal workpieces and time-dependent domains play an important role. Thus, a parabolic heat equation is coupled with other equations. Challenges for the optimization of such coupled systems are presented.

The first model describes a milling process, where material is removed and heat is produced by the cutting. This leads to a thermomechanical distortion of

<footnote>* corresponding author</footnote>
the workpiece during the cutting process and leads to an incorrect removal of material. An optimization of the cutting path and speed, varying chip thickness and thus heat production, etc., should give reduced distortion during the process and lead to a correct workpiece shape.

The second model describes the melting and solidification of metal heated by a laser beam. Due to free-surface flow, the shape of the liquid part depends on capillary boundary conditions, and heat transport on the flow field. The microstructure of the re-solidified material, which is important for subsequent process steps, depends on the temperature gradients near the moving liquid-solid interface and its velocity. Accelerating the process for mass production on the one hand and improving the microstructure on the other hand compete for an optimized process.

Both optimization problems can be formulated in an abstract setting as

$$ \min_{u \in U_{ad}} J(u, y) \quad \text{under the constraint} \quad y = S(u) \qquad (1) $$

where $u$ denotes the control, $U_{ad}$ the set of admissible controls, $y$ the state, $J$ the error functional to be minimized, and $S$ is the control-to-state operator, given by a nonlinear PDE. $J$ is typically given by a deviation of $y$ (or something derived from it) from a desired function $y_d$ plus some regularization by a norm of $u$, like

$$ J(u, y) = d(y, y_d) + \alpha \|u\|_p^p. \qquad (2) $$

For both applications, we will first state and describe the primal problem giving the solution operator $S$, and later cover some details of the associated optimization problem.

# 2 Thermomechanics of Milling Processes

## 2.1 Application

During a milling process, heat is produced by the cutting tool and transfered into the workpiece, and mechanical load is generated by cutting forces. Due to the resulting thermomechanical deformation of the workpiece, the final shape of the processed workpiece deviates from the desired shape, making a postprocessing finishing necessary. Deformations are relatively large especially when producing fine structures like thin walls for lightweight constructions. In order to reduce the shape deviation, an optimization of the tool path and other process parameters is desirable, taking into account the thermomechanical deformations.

## 2.2 Model

The mathematical model for the thermomechanics of the process includes thermoelasticity of the workpiece, energy and forces introduced by the process, and most importantly the cutting process itself, which leads to a time-dependent domain, whose shape influences the process and vice versa. As typical under
the assumption of small deformations, the model is formulated in a reference configuration.

Let $\Omega(t) \subset \mathbb{R}^3$ denote the time-dependent domain in the reference configuration, $Q_T := \{(x,t) : x \in \Omega(t), t \in (0,T)\}$ the space-time-domain, $\theta : Q_T \to \mathbb{R}$ the temperature, and $\mathbf{v} : Q_T \to \mathbb{R}^3$ the deformation.

The change of geometry by material removal as well as energy and forces introduced by the cutting process are provided by a process model [3, 5], taking into account the tool path and velocity, chip thickness, temperature and deformation, and other global and local parameters and properties. Here, we rely on a macroscopic model where microscopic processes like chip formation are not directly considered, but their effects considered via the process model. Let us denote by $\Gamma(t) \subset \partial\Omega(t)$ the contact zone of the cutting tool at time $t$, $\Gamma_T := \{(x,t) : x \in \Gamma(t), t \in (0,T)\}$, $q_T : \Gamma_T \to \mathbb{R}$ the normal heat flux, and $g_T : \Gamma_T \to \mathbb{R}^3$ the forces introduced at the cutting surface.

In the notion of optimal control problems, the state $y$ consists of the domain, temperature, and deformation, $y = (\Omega, \theta, \mathbf{v})$, while the control $u$ consists of the process parameters like tool path, feed rate, rotational velocity, etc. The material removal and thus the domain $\Omega(t)$ depend on the cutting process (control $u$) and the deformation $\mathbf{v}$.

The coupled model includes the parabolic heat equation and quasistatic, elliptic thermoelasticity

$$ \dot{\theta} - \nabla \cdot (\kappa \nabla \theta) = 0, \quad (3) $$

$$ -\nabla \cdot \sigma = f_v(\theta) \quad (4) $$

on $\Omega_T$ with stress tensor $\sigma = 2\mu D\mathbf{v} + \lambda \mathrm{tr}(D\mathbf{v})\mathbb{I}$ and $D\mathbf{v} = \frac{1}{2}(\nabla\mathbf{v} + \nabla\mathbf{v}^T)$ the symmetric gradient or strain tensor. On the contact zone $\Gamma_T$, we have boundary conditions for heat flux and mechanical forces given by the process model,

$$ \kappa \nabla \theta \cdot \mathbf{n} = q_T(u, \theta, \mathbf{v}), \quad (5) $$

$$ \sigma \mathbf{n} = g_T(u, \theta, \mathbf{v}). \quad (6) $$

The workpiece is clamped, which is reflected by Dirichlet conditions on a subset $\Gamma_D \subset \partial\Omega \setminus \Gamma_T$,

$$ \mathbf{v} = 0 \quad \text{on } \Gamma_D, \quad (7) $$

while cooling conditions and free deformation apply on the rest of the boundary,

$$ \kappa \nabla \theta \cdot \mathbf{n} = r(\theta_{\text{ext}} - \theta), \quad (8) $$

$$ \sigma \mathbf{n} = 0. \quad (9) $$

Initially, the temperature is typically constant at room temperature, thus $\theta = \theta_0$ on $\Omega(0)$.

## 2.3 Numerical Discretization of the Forward Problem

The system (3-9) of PDEs is discretized using a finite element method on an adaptively locally refined tetrahedral mesh [10, 13], using piecewise polynomial
functions for the temperature and the components of deformation. Time-discre-
tization is based on a semi-implicit time stepping scheme.

The time-dependent domain Ω(t) is approximated by a subset Ω_h(t) of the triangulation, where the completely cut off elements are ignored. The cutting process is simulated by a dexel method [11], which is able to compute the inter-action of the tool with the (deformed) workpiece very efficiently, giving ˜Ω_h(t) and ˜Γ_h(t). At the same time, chip thickness and other cutting parameters are computed and the process model returns approximations of the heat flux q_Γ,h and forces g_Γ,h at the cutting surface ˜Γ_h(t). Based on that, the finite element method computes Ω_h(t) and Γ_h(t) and projects the boundary data onto Γ_h(t). Finite element approximations of temperature and deformation are computed on Ω_h(t). This is done in every time step of the finite element method. The overall method is described in [4, 5].

Figure 1 shows the mesh, temperature, and deformation from the simulation of a milling process. The mesh is adaptively refined in order to approximate the geometry Ω(t) well by Ω_h(t). The process removes layers of material to mill a pocket into a rectangular bar of metal. Especially the final thin backward wall is prone to deformations larger than the given tolerance, making it hard to produce the desired shape.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Simulation of a milling process: mesh and temperature (top) and deformation (bottom, amplified by a factor 100). The tool is at the moment cutting near the backward left lower corner.</figcaption></figure>
## 2.4 Optimization

The goal of a process optimization is a compensation and reduction of geometry errors while not slowing down the process too much. Control is given by the variation of process parameters like tool path and velocity, speed of tool rotation, etc. Such variations result in changes of the tool entry situation and thus, reflected by the process model, changes in heat source and cutting forces, and finally a change in the thermomechanic deformation. Admissible controls in $U_{ad}$ are defined by restrictions on the machining process.

Given a prescribed final geometry $\Omega_d$ of the workpiece, the optimization functional should include the deviation of the process geometry from the desired one, as well as the process duration. Considering the geometry error, different criteria are possible, especially comparing geometries during the whole cutting process or only in the end. For the latter, this would nevertheless include geometry error terms over time, as material which was removed before cannot be added later on again.

Another approach to an error functional includes geometry deviations near the cutting zone $\Gamma_T$ during the whole process:

$$J(u, \Omega, \theta, \mathbf{v}) = \int_{0}^{T} \|\Omega(t) - \Omega_{d}(t)\|_{\Gamma(t)}^{2} + \lambda\|u\|^{2}. \quad (10)$$

We show the effect of such an error functional in a model situation, where a L-shaped geometry is produced from a rectangular plate. This can be seen as a slice through the original workpiece, see the left of Figure 2. The non-optimized control does not consider the thermal extension and leads to increased material removal resulting in a recessed surface after the workpiece has cooled down. Figure 2 shows on the right the geometry error in the contact zone over the process time. The general optimal control problem requires to find the spatial tool positions and the cutting parameters which minimize (10). In a first investigation of the model problem, simple raster milling is performed and control is given by traditional setting parameters, i.e. cutting depth, radial and tangential feed and cutting velocity. Figure 3 shows the resulting surface for the non-optimized process and for the process with optimal parameters.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Model geometry (left), deformation $\mathbf{v}$ in the cutting zone at $t = 10.5s$ (middle), and maximal deformation in the cutting zone $\Gamma(t)$ over time (right).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** L-shaped workpiece near the end of the milling process, with initial control (left) and optimized parameters (right). Colors depict the modulus of deformation.</figcaption></figure>

**Adjoint problem.** The optimization shown above was done with the standard MATLAB optimization toolbox which just calls the finite element package to solve the forward problem. For a more involved algorithm, the computation and storage of the adjoint solution would be used. The corresponding system of adjoint problems consists of the adjoint (backward) heat equation on the time-dependent domain (now given from the forward problem), coupled to the quasistatic adjoint elasticity equation. Due to the rather long process time, the three-dimensional time-dependent domain, and adaptive meshes for approximation of domain and solution, the handling of such adjoint solutions in the optimization procedure is a challenge.

# 3 Material Accumulation by Laser Heating

## 3.1 Application

For the production of micro components (like micro-valves, etc., with diameters smaller than 1mm) by cold forming, a necessary pre-forming step is to accumulate enough material for a subsequent cold forming step. This can be done by partial melting and solidification of a half-finished product like a thin wire [14]. Due to the small scale, the dominant surface tension of the melted material leads to a nearly spherical form, leading to an accumulated solid sphere attached to the wire after solidification which is called preform. Due to the industrial background, very high process speeds are requested. However, for the subsequent forming step, the microstructure of the material is important. Thus, besides the speed of the process and an accurate size of the preform, its microstructure is part of the optimization goal. Formation of dendritic structures and their spacings, or other phases, during the solidification are strongly influenced by the liquid-solid interface velocity and local temperature gradients [9]. Thus, temperature, phase transitions, and the geometry are important aspects of the corresponding optimization problem.
## 3.2 Model

We consider the time-dependent domain $\Omega(t)$ consisting of the solid and partially melted parts of the material (metal). Melting and solidification are typically modeled by the Stefan problem, including temperature $\theta$ and energy density $e$ as variables in the space-time-domain $Q_T := \{(x, t) : x \in \Omega(t), t \in (0, T)\}$:

$$ \dot{e} + \mathbf{v} \cdot \nabla e - \nabla \cdot (\kappa \nabla \theta) = 0, \quad \theta = \beta(e), \qquad (11) $$

where $\beta(s) := c_1 \min(s, 0) + c_2 \max(0, s - L)$ and $L$ denotes the latent heat of solid-liquid phase transitions. $\beta$ is only a piecewise-smooth function with a constant part, making (11) a degenerate parabolic equation.

The liquid subdomain $\Omega_l(t)$ is given by all points where the temperature is above the melting temperature (which is assumed to be 0 after some scaling),

$$ \Omega_l(t) := \{x \in \Omega(t) : \theta(x, t) > 0\}. \qquad (12) $$

The shape of the melted (and later on re-solidified) material accumulation is mainly influenced by the surface tension of the liquid, together with gravitational forces etc., which means free-surface flow. Parabolic Navier-Stokes equations with capillary boundary condition is the main model component for this, with solenoidal velocity field $\mathbf{v}$ and pressure $p$ in $\Omega_l(t)$,

$$ \dot{\mathbf{v}} + \mathbf{v} \cdot \nabla \mathbf{v} - \nabla \cdot \sigma = f_v(\theta), \quad \nabla \cdot \mathbf{v} = 0, \qquad (13) $$

with stress tensor $\sigma = \frac{1}{Re}D\mathbf{v} - p\mathbb{I}$. Here, $f_v$ denotes the forces introduced by gravity due to a temperature-dependent density and the Boussinesq approximation. The shape of the liquid subdomain is given through the capillary boundary condition on the free surface $\Gamma(t)$ of the melted subdomain, where the surface tension (proportional to the mean curvature of the surface) balances the normal stress. In a differential geometric PDE formulation, the mean curvature vector $H_\Gamma$ is given by the Laplace-Beltrami-operator $-\Delta_\Gamma$ applied to the coordinates of the surface (represented by the embedding $id : \Gamma(t) \to \mathbb{R}^3$), giving another (nonlinear) elliptic equation in the coupled system. Additionally, the normal component of the fluid velocity should be equal to the normal velocity $V_\Gamma$ of the capillary surface. Both lead to the following equations on $\Gamma(t)$:

$$ \sigma \mathbf{n} = \frac{1}{We} H_\Gamma = -\frac{1}{We} \Delta_\Gamma id, \quad \mathbf{v} \cdot \mathbf{n} = V_\Gamma. \qquad (14) $$

On the solid-liquid interface and in the solid subdomain, the velocity vanishes,

$$ \mathbf{v} = 0 \quad \text{in } \Omega(t) \setminus \Omega_l(t). \qquad (15) $$

The heating is done through a laser pointing at a spot on the boundary, modeled by a time- and space-dependent energy density $q_L$, and cooling conditions apply on the whole boundary,

$$ \kappa \nabla \theta \cdot \mathbf{n} = q_L + r(\theta_{ext} - \theta) \quad \text{on } \partial\Omega(t). \qquad (16) $$

A more detailed description of the model can be found in [8].
## 3.3 Numerical Discretization of the Forward Problem

The forward problem with given boundary conditions is discretized by a finite element method which combines a Stefan problem solver with a free-surface Navier-Stokes solver. The latter is based on the Navier code [1], the combined approach is described in more detail in [2, 7]. Locally refined (triangular or tetrahedral) meshes are needed in order to approximate the large variations in temperature near the heating zone and the surrounding of the solid-liquid interface sufficiently well, while keeping the overall numerical costs acceptable. Due to the big changes in geometry, starting from a thin wire and ending in a relatively large spherical accumulation, several remeshings are necessary during the simulation in order to avoid a degeneration of mesh elements.

Figure 4 shows a typical mesh, temperature field, and velocity field during the melting, with the laser pointing to the center bottom of the material. Due to rotational symmetry, a 2D FEM with triangular meshes could be used. As the wire is melted from below, the growing sphere is moving upwards and thus the velocity vectors are pointing upwards, too. In Figure 5, we show several stages

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Melting the end of a thin wire: liquid material accumulation with mesh, solid-liquid interface, and liquid flow field.</figcaption></figure>

during the solidification after the heating is switched off (from a simulation with shorter heating period). In contrast to our results for energy dissipation during the heating process [6], the cooling after switching off the laser is mainly done by conducting heat upwards into the wire, which results in a downward movement of the interface. The additional cooling by the boundary conditions leads later on to an additional solidification at the boundary of the liquid region. During solidification, the shape does not change much anymore, so the velocities are typically quite small and not shown here.

## 3.4 Optimization

The goals for an optimization of the process are on the one hand a high speed in order to be able to produce thousands of micro preform parts in a short time,
<figure><img src="image_6.png" /><figcaption>Fig. 5. Material accumulation at the end of a thin wire: Solid-liquid interface and isothermal lines at 5 different times during solidification.</figcaption></figure>

while generating material accumulations of the desired size, and on the other hand generating a microstructure which is well suited for the subsequent cold forming step and usage properties of the work piece. Models for microstructure evolution during solidification show a dependency of the microstructure quality on the speed $V_{\Gamma_I}$ of the liquid-solid interface $\Gamma_I(t)$ and the temperature gradients near the interface [9]. A higher speed and larger gradient typically result in a microstructure which gives better forming characteristics and useful properties. Thus, the error functional has to include parts for geometry approximation, for speed, and for the microstructure generation during solidification,

$$ J(u, \Omega, \theta, \mathbf{v}) = \int_{0}^{T} \|\Omega(t) - \Omega_{d}(t)\|^{2} \\ + \lambda_{1} \int_{0}^{T} \left( \|V_{\Gamma_I} - V_{d}\|_{\Gamma_I(t)}^{2} + \|\nabla\theta - G_{d}\|_{\Gamma_I(t)}^{2} \right) + \lambda_{2}\|u\|^{2}. \quad (17) $$

The control parameters $u$ are given by the time-dependent intensity and location of the laser spot, thus they enter the system of equations via the heating energy density $q_L(u)$. An additional control variable might be the intensity $r$ of outer cooling, for example by adjusting the flow velocity of a cooling gas. $U_{ad}$ is given by control restrictions which follow from limits of the process, like an upper bound of the laser power for preventing an evaporation of the material.

**Adjoint Problem.** For an efficient implementation of the optimization procedure, the solution of the adjoint problem will be used. Here, the system of adjoint equations includes the adjoint (linearized) Navier-Stokes system on the prescribed time dependent domain (from the forward solution), together with the (linear) Laplace-Beltrami equation on the prescribed capillary boundary. For the formulation of the adjoint Stefan problem, a regularization could be used which was used in [12] for the derivation of an aposteriori error estimate.

As in our first application, also here the efficient handling of the adjoint problem including adaptively refined meshes with remeshings for the time dependent domain and the corresponding system of solutions in the domain and on the capillary surface poses a challenge for the overall numerical optimization method.
**Acknowledgments.** The authors gratefully acknowledge the financial support by the DFG (German Research Foundation) for the subproject A3 within the Collaborative Research Center SFB 747 "Micro cold forming" and the project MA1657/21-3 within the Priority Program 1480 "Modeling, Simulation and Compensation of thermal effects for complex machining processes".

Furthermore, we thank the Bremen Institute for Applied Beam Technology (BIAS) and the Institute of Production Engineering and Machine Tools Hanover (IFW) for cooperation.

References

1. Bänsch, E.: Finite element discretization of the Navier-Stokes equations with a free capillary surface, Numer. Math. 88, 203-235 (2001)

2. Bänsch, E., Paul, J., Schmidt, A.: An ALE finite element method for a coupled Stefan problem and Navier-Stokes equations with free capillary surface, Int. J. Numer. Meth. Fluids 71, 1282-1296 (2013)

3. Denkena, B., Maaß, P., Niederwestberg, D., Vehmeyer, J.: Identification of the specific cutting force for complex cutting tools under varying cutting conditions, Int. J. of Maschine Tools and Manufacture 82-83, 42-49 (2014)

4. Denkena, B., Schmidt, A., Henjes, J., Niederwestberg, D., Niebuhr, C.: Modeling a Thermomechanical NC-Simulation, Procedia CIRP 8, 69-74 (2013)

5. Denkena, B., Schmidt, A., Maaß, P., Niederwestberg, D., Niebuhr, C., Vehmeyer, J.: Prediction of Temperature Induced Shape Deviations in dry Milling, Procedia CIRP 31, 340-345 (2015)

6. Jahn, M., Brüning, H., Schmidt, A., Vollertsen, F.: Energy dissipation in laser-based free form heading: a numerical approach. Prod. Eng. Res. Devel. 8, 51-61 (2014)

7. Jahn, M., Luttmann, A., Schmidt, A., Paul, J.: Finite element methods for problems with solid-liquid-solid phase transitions and free melt surface. PAMM 12, 403-404 (2012)

8. Jahn, M., Schmidt, A.: Finite element simulation of a material accumulation process including phase transitions and a capillary surface. Tech. Rep. 12-03, ZeTeM, Bremen (2012)

9. Kurz, W., Fisher, J.D.: Fundamentals of Solidification. Trans Tech Publications, Switzerland-Germany-UK-USA (1986)

10. Niebuhr, C., Niederwestberg, D., Schmidt, A.: Finite Element Simulation of macroscopic Machining Processes - Implementation of time-dependent Domain and Boundary Conditions, 14 p., Berichte aus der Technomathematik 14-01, University of Bremen (2014)

11. Niederwestberg, D., Denkena, B., Böß, V., Ammermann, C.: Contact Zone Analysis Based on Multidexel Workpiece Model and Detailed Tool Geometry, Procedia CIRP 4, 41-45 (2012)

12. Nochetto, R.H., Schmidt, A., Verdi, C.: A posteriori error estimation and adaptivity for degenerate parabolic problems. Math. Comput. 69, 1-24 (2000)

13. Schmidt, A., Siebert, K.G.: Design of adaptive finite element software: The finite element toolbox ALBERTA, Springer LNCSE Series 42 (2005)

14. Vollertsen, F., Walther, R.: Energy balance in laser based free form heading, CIRP Annals 57, 291-294 (2008)