TOPOLOGICAL DERIVATIVES AND A LEVEL SET APPROACH FOR AN INVERSE ELECTROMAGNETIC CASTING PROBLEM

In this paper we describe a new method for the topology design of the inductors used in the electromagnetic casting technique of the metallurgical industry. The method is based on a recently proposed topology optimization formulation of the inverse electromagnetic casting problem, and uses level-sets together with first and second order topological derivatives of the objective functional, which are herein given, to efficiently find the optimal solution. Results for two numerical examples show that the technique described can be efficiently used in the design of suitable inductors.


INTRODUCTION
The electromagnetic casting is an industrial technique that allows for contactless heating, shaping and controlling of chemical aggressive, hot melts.It makes use of the repulsive forces that an electromagnetic field produces on the surface of a mass of liquid metal.In the presence of an electromagnetic field, the liquid metal changes its shape until an equilibrium relation between the electromagnetic pressure and the surface tension is satisfied.The direct problem in electromagnetic casting consists in determining the equilibrium shape of the liquid metal.In general, this problem can be solved either directly studying the equilibrium equation defined on the surface of the liquid metal, or minimizing an appropriate energy functional.The main advantage of this last method is that the resulting shapes are mechanically stable [1][2][3].
The inverse problem consists in determining the electric currents and the induced exterior field, for which the liquid metal takes on a given desired shape.In [4] and [5] we studied the inverse electromagnetic casting problem considering inductors made of a single solid-core wire with a negligible area of the cross-section, and the more realistic case where each inductor is a set of bundled insulated strands.In both cases the number of inductors was fixed in advance.In a recent paper we overcome this constraint, and looked for configurations of inductors considering different topologies with the purpose of obtaining better results [6].
Following the general treatment of the inverse problem considered in [6], we formulate the inverse problem as an optimization problem.The topological asymptotic expansion of the shape functional considered is given, and an efficient topology optimization algorithm based on the level-set technique proposed by Amstutz [7] is used to solve the inverse problem.
Two examples are presented in Section 6 to show the that the optimization method described can efficiently find the solution to a given design problem.Finally, the conclusions are presented in Section 7.

THE MATHEMATICAL MODEL
The simplified model of the electromagnetic casting problem studied here concerns the case of a vertical column of liquid metal falling down into an electromagnetic field created by vertical inductors.We consider the quasi-static model and assume that the frequency of the imposed current is very high so that the magnetic field does not penetrate into the metal.We assume also that a stationary horizontal section is reached so that the two-dimensional model is valid.The equilibrium of the system is ensured by the static balance on the surface of the metal between the surface tension and the electromagnetic forces.This problem and other similar ones have been considered by several authors, we refer the reader to papers [1,2,[8][9][10][11] for the physical analysis of the simplifying assumptions of the model.
We denote by Ω ⊂ R 2 the exterior of the domain ω occupied by the cross-section of the metal column, which is assumed closed, simply connected, and with a non-void interior.The exterior boundary value problem regarding the inverse electromagnetic casting problem in terms of the flux function ϕ : Ω → R is: where µ 0 is the vacuum permeability, j 0 is the component of the current density vector normal to the cross-section, that we assumed compactly supported in Ω and such that the total current is zero: Γ is the boundary of Ω, • denotes the Euclidean norm, and the constant c is the value at infinity of the solution ϕ, which is also an unknown of (1) [4][5][6].
Problem (1) has unique solutions ϕ ∈ W 1 0 (Ω) and c ∈ R [12,13], where W 1 0 (Ω) is defined as: The cross-section area of the liquid metal column is known and equal to S 0 : The equilibrium of the liquid metal boundary between the magnetic pressure and the surface tension is characterized by the following equation [11,[14][15][16][17]: where C is the curvature of Γ seen from the metal, σ is the surface tension of the liquid and the constant p 0 is an unknown of the problem.Physically, p 0 represents the difference between the internal and external pressures.
In the direct problem the electric current density j 0 and the cross-section area S 0 are given, and one needs to find the shape of ω satisfying (4) such that the flux function ϕ solution to (1) satisfies also the equilibrium equation ( 5) for a real constant p 0 .

The Inverse Problem
In the inverse problem we have to determine the current density j 0 satisfying (2) such that the solution ϕ of (1) satisfies also the equilibrium equation (5).If the inverse problem has an exact solution, we say that the target shape is shapable, if it does not, we say that it is not shapable.Even considering a shapable shape, the inverse problem is inherently ill posed: small variation of the liquid boundary may cause dramatic variations in the solution j 0 of the inverse problem [9,18].In addition, the uniqueness of the solution in terms of j 0 cannot be ensured [6].Hence, we follow the approach proposed in [6], where the inverse problem is formulated as an optimization problem, in order to look for a solution (maybe just an approximate solution) minimizing an appropriate functional.
There are some known facts about the exact solutions of the inverse problem that are of main importance in what follows.First, at the solution the constant p 0 satisfies [6,9]: Hence, for a given target shape, the value of p 0 can be calculated using (6).Second, calling p = 2µ 0 (p 0 − σC), the equilibrium equation in terms of the flux function reads where κ = ±1, with the sign changes located at points where the curvature of Γ is a global maximum.The two possible ways of defining κ lead to the same solution j 0 but with the opposite sign [6].Third, a necessary condition for the existence of a solution to the inverse problem is the following [6]:

TOPOLOGICAL DERIVATIVE CONCEPT
The topological derivative measures the sensitivity of a given shape functional with respect to an infinitesimal singular domain perturbation, such as the insertion of holes, inclusions, source-terms or even cracks.The topological derivative was introduced by Eschenauer et al. [19] and rigorously investigated by Sokołowski and Żochowski [20] and by Céa et al. [21].Since then, this concept has proved extremely useful in the treatment of a wide range of problems, such as topology optimization [7,[22][23][24][25][26], and inverse analysis [27][28][29][30][31][32].
Let us consider that the domain Ω is subject to a non-smooth perturbation confined in a small ball B ε ( x) of radius ε and center x ∈ Ω.Then, if a given shape functional ψ(ε), associated to the topologically perturbed domain, admits the following topological asymptotic expansion [20] where ψ(0) is the shape functional associated to the original (unperturbed) domain and f i (ε), 1 ≤ i ≤ 2, are positive functions such that f i (ε) → 0, and f 2 (ε)/f 1 (ε) → 0, when ε → 0, we say that the functions x → D i T ψ( x), 1 ≤ i ≤ 2, are the topological derivatives of ψ at x (in this paper we call first and second order topological derivatives to T ψ can be seen as a second order correction of ψ(0) to approximate ψ(ε).In fact, the topological derivatives are scalar functions defined over the original domain that indicate, at each point, the sensitivity of the shape functional when a singular perturbation of size ε is introduced at that point.
In this paper we describe a new method for inverse electromagnetic casting problems based on the topological expansion (9), which is presented in details in the next section.

PROBLEM FORMULATION
The considerations made in Section 2.1 allow us to formulate the inverse problem as follows: for a given domain Ω limited by a smooth boundary Γ, determine an electric current density j 0 satisfying (2) with compact support spt(j 0 ) ⊂ Θ, where Θ is a given compact in Ω [9], and a real constant c such that the system has a solution ϕ ∈ W 1 0 (Ω).We assume that condition (8) is satisfied.Note that the above problem is ill-posed, since the boundary conditions on Γ are overdetermined.Therefore, the idea is to rewrite the inverse problem (10) as an optimization problem.
Consider the following shape functional based on the Kohn-Vogelius criterion [6], where the auxiliary function φ depends implicitly on j 0 and c by solving the following boundary-value problem where, denoting |Γ| = Γ ds, d(j 0 ) is defined as Note that φ can be computed for a current density j 0 not necessarily satisfying condition (2).The approach described here to deal with (10) is the following: determine an electric current density j 0 satisfying (2) and the constant c such that the solution φ to (12) minimizes the shape functional (11).We note that the minimum of the shape functional ( 11) is attained when φ ≡ 0 on Γ.This means that in this situation, from the well-posedness of problems ( 1) and ( 12), we have φ ≡ ϕ in Ω.
In a first step we can eliminate the variable c of the optimization problem defining it as the global minimum c * (j 0 ) of ( 11) for any fixed j 0 , i.e., we take c = c * (j 0 ) = arg min c J(φ(j 0 , c)).It is easy to see that the global minimum with respect to c is such that Γ φ ds = 0 [6].Hence, we can formulate an equivalent optimization problem as follows: minimize the shape functional (11), where φ ∈ W 1 0 (Ω) depends implicitly on j 0 only, by solving the following problem:

The topological derivative calculation
Associated to φ we define the function φ ε solution to the perturbed problem.In this context, the perturbation is characterized by changing the electric current distribution j 0 by a new one j ε which is identical to j 0 everywhere in Ω except in a small region B ε (x) ⊂ Ω, where B ε (x) denotes a ball of radius ε and center x.More precisely, j ε is given by where I is a given current density value and α = ±1 is the sign of the current density in B ε (x).In this way, the shape functional associated to the perturbed problem reads: where φ ε is unique solution in W 1 0 (Ω) to the following problem: Theorem 1.The functional ψ admits the topological expansion with the functions f 1 (ε) = πε 2 , f 2 (ε) = π 2 ε 4 and the topological derivatives where the function f ∈ W 1 0 (Ω) satisfies the following problem: Alternatively, the first order topological derivative can be computed efficiently using the adjoint state v, see [6], where v is the unique solution in W 1 0 (Ω) to the following problem:

A TOPOLOGICAL DERIVATIVE-BASED LEVEL-SET ALGORITHM
In this section we state a level-set topology design algorithm based on the topological derivatives obtained in the previous section.We consider domains Ω which possess a central symmetry.In this case we can devise a level set algorithm that generate a sequence of current density functions j 0 satisfying (2) at each iteration.In this case the current density j 0 is sought as the solution of the general optimization problem stated as It should be stressed that the design variable in problem (24) is the topology of the support of j 0 , and its sign.Hence, the use of the exact topological sensitivity information provided by the topological derivatives (19) and (20) emerges as a natural alternative in the development of a numerical optimization algorithm to tackle the problem.
We decompose the working compact set Θ ⊂ Ω into disjoint parts Θ + , Θ − , and Θ 0 , representing the regions of Θ where the current density j 0 is, respectively, positive, negative or zero.For clarity of exposition, let us assume that Θ + , Θ − , and Θ 0 are open domains (perturbations with center x on the domain boundaries will not be considered).Introduction of a circular region of current density αI and center x ∈ Θ 0 changes the value of the objective function of problem (24) according to the topological expansion (18).Note also that the topological expansion corresponding to introduction of a circular region of zero current and center x ∈ Θ + is also given by ( 18) with α = −1.For a circular region of zero current and center x ∈ Θ − we can use (18) with α = +1.
Consider now an optimal configuration of inductors with respect to the class of perturbations we are analyzing.Introduction of a small circular region of current density αI and center x ∈ Ω 0 = Ω − (Θ + ∪ Θ − ) should not increase the objective function.Hence, according to the expression of the first order topological derivative, an optimal configuration should satisfy the following necessary condition: Since α could be positive or negative, for the optimal configuration we have Since v is harmonic, it must be zero on the entire domain Ω, and by ( 19) the first order topological derivative D 1 T ψ vanishes in Ω for any optimal configuration, i.e., Having made the above considerations, we are ready to devise e topological derivativebased optimization algorithm based on the ideas presented in [7] to solve problem (24).The procedure relies on a level-set domain representation [33] and approximation of the topological optimality conditions by a fixed point iteration.In particular, the algorithm presented in [7] displays a marked ability to produce general topological domain changes uncommon to other methodologies based on a level-set representation and has been successfully applied in [7] to topology optimization in the context of two-dimensional elasticity and flow through porous media.
The main difficulty of applying the ideas in [7] to the present case is that the first order topological derivative vanishes at the solution according to (27).Thus, the first order topological derivative alone cannot be used to define the level sets of the fixed point algorithm stated in [7].Fortunately, the expected variation of the objective functional given by the topological expansion (18) can be used instead of the first order topological derivative to overcome this difficulty.
Following [7], to devise a level-set-based algorithm whose aim is to produce a topology that satisfies ( 31)-( 33), we choose a value for ε (in the numerical approach we define a mesh of cells in the domain Θ and take a value for ε related to the size of the cells) and define the functions With the above definitions and ( 28)-( 29), it can be easily established that the sufficient conditions (31)- (33) are satisfied if the following equivalence relations between the functions g + and g − and the level-set functions ψ + and ψ − hold where h : R → R must be an odd and strictly increasing function, e.g., In fact, if x ∈ Θ + , then ψ + (x) < 0. By (36) and since h preserves the sign we have g + (x) < 0, and by (34) we have EV (x, ε, −1) > 0, proving (31).The proofs of ( 32) and ( 33) are analogous considering x respectively in Θ − and Θ 0 .Let us see how to find the optimal topology of Θ + .The case of Θ − is completely analogous.Starting from an initial level-set function ψ + 0 ∈ L 2 (Θ) which defines the chosen initial guess for the optimum topology of Θ + , the algorithm produces a sequence (ψ + i ) i∈N of level-set functions that provides successive approximations to the sufficient condition for optimality.The sequence satisfies where co(ψ + n , h(g + n )) is the convex hull of {ψ + n , h(g + n )}.In the actual algorithm the initial guess ψ + 0 and subsequent iteration points can be normalized, see [7] for further details.The non-normalized sequences for both Θ + and Θ − are where t n ∈ [0, 1] is a step size determined by a line-search in order to decrease the value of the cost functional J(φ).The iterative process is stopped when for some iteration the obtained decrease in J(φ) is smaller than a given numerical tolerance.The angles θ + and θ − can be use as indicators of the accuracy of (36) and (37) at the final iteration and can be used to determine the necessity of a mesh refinement [7].
If there is no previous information about the optimum, a suitable initial guess that satisfies the inequality ψ + 0 + ψ − 0 ≥ 0 is ψ + 0 = ψ − 0 = 1, which corresponds to zero positive and negative currents.

NUMERICAL EXAMPLES
Two numerical examples are presented.The first example has a known solution, since the target shape considered is an equilibrium shape for given current density distributions.The next example consists of a more realistic design problem.The same boundary element method used in [6] was used here to solve the boundary value problems and to compute approximately the first and second order topological derivatives.
The target shape of the first example is the solution of a direct free-surface problem for a liquid metal column of cross-section area S 0 = π, considering six distributed currents of density I = 0.4 as shown in Fig. 1.The result obtained for a mesh of size 0.02, defined in the region shown in Fig. 1, is shown in Fig. 2.

CONCLUSIONS
We have described a new method for the topology design of inductors in electromagnetic casting.The method is based on the topology optimization formulation presented in [6] and uses level-sets together with first and second order topological derivatives to design suitable inductors.
The complete asymptotic expansion of the objective functional regarding the introduction of a small inductor was given in this paper, generalizing the results of [6].In the case of centrally symmetric geometries, the method described generates a sequence of solutions satisfying all the equality constraints.The set of examples considered show that the method described is effective and efficient, and therefore can be successfully used in the design of inductors in electromagnetic casting.

Figure 1 .
Figure 1.Example 1.(a) Initial configuration of the direct free-surface problem.(b) Target shape and exact solution.Black area: positive inductors, gray area: negative inductors, dashed line: target shape, thin solid line: boundary of the mesh of cells.

Figure 2 .Figure 3 .
Figure 2. Example 1. Solution for a mesh of cells of size 0.02 with β = 3. Black area: positive inductors, gray area: negative inductors, dashed line: target shape, thin solid line: equilibrium shape.

Figure 4 .
Figure 4. Example 2. Solution for a mesh of cells of size 0.02 and β = 3. Black area: positive inductors, gray area: negative inductors, dashed line: target shape, thin solid line: equilibrium shape.