Ngoc Anh Nguyen

Sorin Olaru

Pedro Rodriguez-Ayerbe

Inverse parametric

Inverse parametric linear/quadratic programming problem for continuous PWA functions defined on polyhedral partitions of polyhedra N. A. Nguyen 1 , S. Olaru 1 , P. Rodriguez-Ayerbe 1 Abstract-Constructive solution to inverse parametric linear/quadratic programming problems has recently been investigated and shown to be solvable via convex liftings [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF], [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF]. These results were stated and solved starting from polytopic partitions of a polytope in the parameter space. Therefore, the case of polyhedral partitions of unbounded polyhedra, was not handled by this method and deserves a complete characterization to address the general inverse optimality problem. This paper has as main objective to overcome the unboundedness limitation of the given polyhedral partition and to extend the constructive solution put forward in [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF] for this omitted case.

I. INTRODUCTION

Inverse parametric convex programming was previously investigated in [START_REF] Baes | Every continuous nonlinear control system can be obtained by parametric convex programming[END_REF] where it was proved that every continuous nonlinear function can be obtained via parametric convex programming without a guideline for the construction of such convex programming problems. The first constructive results have recently been presented via different approaches in [START_REF] Hempel | Every continuous piecewise affine function can be obtained by solving a parametric linear program[END_REF], [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF], [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF]. These developments focus on continuous piecewise affine functions, a particular but important class of nonlinear functions which receives great attention in control theory as attested by the references in different topics [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF], [START_REF] Johansen | On multi-parametric nonlinear programming and explicit nonlinear model predictive control[END_REF], [START_REF] Seron | Characterisation of receding horizon control for constrained linear systems[END_REF], [START_REF]Approximate explicit receding horizon control of constrained nonlinear systems[END_REF], [START_REF] Pistikopoulos | Multi-parametric programming[END_REF], [START_REF] Bemporad | An algorithm for approximate multiparametric convex programming[END_REF], [START_REF] Olaru | A parameterized polyhedra approach for explicit constrained predictive control[END_REF], [START_REF] Pistikopoulos | Perspectives in multiparametric programming and explicit model predictive control[END_REF], [START_REF] Grancharova | Explicit Nonlinear Model Predictive Control[END_REF].

It is well known that optimal solution to a parametric linear/quadratic programming problem is a piecewise affine (PWA) function defined over a polyhedral partition. Note that the continuity of optimal solution to parametric linear programming problems, is not guaranteed in many cases. However, it is proved in [START_REF]On the continuity and complexity of control laws based on multiparametric linear programs[END_REF] that a continuous optimal solution can always be selected. Conversely, inverse parametric linear/quadratic programming (IPL/QP) aims to recover an appropriate optimization problem, described by a constraint set and a cost function, having as the optimal solution a given continuous PWA function defined over a polyhedral partition.

This paper completes and generalizes the result in [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF], [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF] wherein IPL/QP problems are solved via convex liftings. In the respective works, this solution covers polytopic partitions (i.e. bounded polyhedral partitions) of the parameter space, since the convex lifting construction is based on continuity and convexity conditions on vertices of the polytopes composing the given polytopic partition. Therefore, if the given polyhedral partition contains unbounded polyhedra, N. A. Nguyen, S. Olaru, P. Rodriguez-Ayerbe are with the Laboratory of Signals and Systems (L2S, UMR CNRS 8506), CentraleSupélec-CNRS-Université Paris Sud, 3 rue Joliot Curie, Plateau de Moulon, 91192, Gif-sur-Yvette, France. Ngocanh.Nguyen, Sorin.Olaru, Pedro.Rodriguez-Ayerbe @supelec.fr this algorithm cannot be directly employed. This paper is motivated by this limitation to solve IPL/QP problems for more general class of continuous PWA functions, namely continuous PWA functions defined over polyhedral partitions of polyhedra. This result is useful to reduce implementation complexity of PWA control laws as shown in [START_REF]On the complexity of the convex liftings-based solution to inverse parametric convex programming problems[END_REF], [START_REF] Gulan | Implications of inverse parametric optimization in model predictive control[END_REF].

II. PROBLEM FORMULATION A. Generalities and basic notions

Throughout the paper, R, R + , R >0 , N, N >0 denote the field of real numbers, the field of non-negative real numbers, the set of strictly positive real numbers, the set of nonnegative integers and the positive integer set, respectively. The following index set is also defined for ease of presentation with respect to a given N ∈ N >0 :

I N = {i ∈ N >0 | i ≤ N }.
Given a d ∈ N >0 , by 1 d , we denote a vector of dimension d whose elements are equal to 1; i.e.

1 d ∈ R d .
Given a finite set S := {s 1 , . . . , s L }, then Card(S) denotes the cardinal number of the set S, i.e. Card(S) = L. conv(S) denotes the convex hull of S. Moreover, if R is a set of vectors, by cone(R), we denote the cone which has elements of R as its extreme rays.

A polyhedron is defined as the intersection of a finite number of halfspaces. A polytope is defined as a bounded polyhedron. Given a full dimensional polyhedron S ⊂ R d , then V(S) denotes the set of its vertices, and R(S) denotes the set of its extreme rays. We use k-face to denote a face of dimension k of S. A 0-face is also called a vertex, an 1-face is alternatively called an edge, a (d -1)-face amounts to a facet. Furthermore, by int(S), we denote the interior of S. If S is a vector space, and S is an arbitrary set, then Proj S S represents the orthogonal projection of set S onto space S. Also, we use dim(S) to denote the dimension of the affine hull of S.

Finally, given two sets P 1 , P 2 ⊂ R d , by P 1 ⊕P 2 , we denote the Minkowski addition defined as follows:

P 1 ⊕ P 2 = y ∈ R d | ∃x 1 ∈ P 1 , x 2 ∈ P 2 s.t. y = x 1 + x 2 .
The following definitions are extended in comparison with their counterparts presented and employed in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF], [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF], [START_REF]On the complexity of the convex liftings-based solution to inverse parametric convex programming problems[END_REF].

Definition 1: A collection of N ∈ N >0 full-dimensional polyhedra X i ⊂ R d is called a polyhedral partition of a polyhedron X in R d if: 1) X = i∈I N X i , 2) int(X i ) int(X j ) = ∅ with i = j, (i, j) ∈ I 2 N , Also, (X i , X j ) are called neighbours if (i, j) ∈ I 2 N , i = j and dim(X i ∩ X j ) = d -1.
If set X is compact, then this partition is called polytopic partition.

The definition of a cell complex was presented by Grünbaum in [START_REF] Grünbaum | Convex polytopes[END_REF]. For simplicity, a cell complex should be understood as a polyhedral partition whose facet-to-facet property is fulfilled; i.e. any pair of neighboring regions share a common facet.

In the remainder of this paper, a slight abuse of notation is used: X implies simultaneously a polyhedral partition/cell complex of a polyhedron and this polyhedron itself. Its meaning will be clarified from the context.

Definition 2: For a given polyhedral partition X = i∈I N X i ⊆ R d , a picewise affine lifting is described by a real-valued function z : X → R with:

z(x) = a T i x + b i for any x ∈ X i , (1)
and

a i ∈ R d , b i ∈ R, ∀i ∈ I N . Definition 3: Given a polyhedral partition X = i∈I N X i ⊆ R d , a piecewise affine lifting z(x) = a T i x + b i ∀x ∈ X i , is called convex piecewise affine lifting if the following conditions hold true: • z(x) is continuous over X , • for each i ∈ I N , z(x) > a T j x + b j for all x ∈ X i \X j and all j = i, j ∈ I N .
For ease of presentation, a convex piecewise affine lifting is called, throughout the rest of this paper, convex lifting. It is shown in [START_REF] Rybnikov | Polyhedral partitions and stresses[END_REF], [START_REF]Recognition of additively weighted voronoi diagrams and weighted delaunay decompositions[END_REF] that a polyhedral partition which admits a convex lifting, is a cell complex. Therefore, a convex lifting is always defined over a cell complex instead of a polyhedral partition.

For a given polyhedral partition, some conditions have to be fulfilled for the existence of a convex lifting. A summary of these conditions can be found in [START_REF] Rybnikov | Polyhedral partitions and stresses[END_REF]. We extract here only the necessary notation in connection with the present developments. Interested readers can find details in [START_REF] Rybnikov | Polyhedral partitions and stresses[END_REF], [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF].

Given a polyhedron C and its facet F , we use n(F, C) to denote the inward unit normal vector at F to C, meaning the unit vector normal to F at F and inward to C. The definition of stress is recalled here for completeness.

Definition 4:

(Rybnikov [20]) A real-valued function s(•) defined on the (d -1)-faces of a cell complex X ⊆ R d is called a d-stress if at each internal (d -2)-face F of X : F ⊂C s(C)n(F, C) = 0, (2)
where this sum ranges over all (d-1)-faces in the star of F (the (d -1)-faces such that F is their common facet). The quantities s(C) are the coefficients of d-stresses, and are called tension if the sign is strictly positive, and compression if the sign is strictly negative. Such a (d -2)-face F of cell complex X satisfying ( 2) is also called in equilibrium. Note also that a necessary and sufficient condition for such a cell complex X to be convexly liftable is the existence of a strictly positive d-stress defined over all (d -1)-faces of cell complex X .

B. Problem formulation

The solution to inverse parametric linear/quadratic programming problems is shown in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF], [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF] via convex liftings. This solution restricts the parameter space to a polytope and lets aside the important case of unbounded parameter space. This omission is tackled in this paper and the problem can be formally stated as follows: Problem statement: Given a polyhedral partition of a polyhedron X = i∈I N X i ⊆ R dx , we use f pwa (•) to denote a continuous PWA function: f pwa (x) : X → R du , defined over partition X . The objective is to determine

J(x, z, u), H x , H u , H z , K such that:    f pwa (x) = Proj R du arg min [z u T ] T J(x, z, u), s.t: H x x + H z z + H u u ≤ K. (3)
Some assumptions need to be stated for guaranteeing the validity of the proposal. Assumption 1: Polyhedral partition X is a cell complex. Assumption 2: Polyhedral partition X is convexly liftable. As discussed before, a polyhedral partition, admitting a convex lifting, is a cell complex, therefore, if a given partition satisfies Assumption 2, then it also fulfills Assumption 1.

Note that Assumption 2 is not restrictive, since, as shown in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF], any polyhedral partition can be subdivided into a convexly liftable cell complex such that the internal boundaries are maintained. This inclusion can be easily extended to unbounded polyhedral partitions by following the same methodology of the proof (Theorem IV.2 in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF]).

III. CONSTRUCTIVE SOLUTION

It is shown in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF], [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF] that constructive solution to such an inverse parametric linear/quadratic programming problem can be built upon convex lifting principle. However, the algorithm to construct convex liftings of the given cell complex put forward therein, is only applicable for cell complexes of polytopes. Therefore, for the cell complexes of polyhedra, this algorithm cannot be directly applicable. To adapt to unbounded parameter space, two intuitive approaches can be proposed:

• either to modify this algorithm subject to cell complexes of polyhedra, • or to adjust the given cell complex of polyhedra in such a manner that this algorithm can be directly of use. The latter approach will be of interest via a simply intermediate step in the sequel. More precisely, we restrict the given parameter space partition to an appropriate bounded region such that the convex liftability of this new partition and the initial one is treated in a similar manner. Accordingly, we propose here the use of a hyperbox large enough to contain in its interior all vertices of the initial partition. The existence of such a hyperbox is guaranteed by the fact that a polyhedral partition is the collection of a finite number of polyhedra, as a consequence the number of vertices is finite. By B dx ( ), we denote a hyperbox for a given scalar ∈ R >0 , defined as

B dx ( ) := x ∈ R dx | x ∞ ≤ .
The hyperbox used for our purpose needs to meet the following assumption: Assumption 3: Given a polyhedral partition X = i∈I N X i ⊆ R dx , the scalar is restricted to the values that ensure x ∈ int(B dx ( )) for every x ∈ i∈I N V(X i ).

The following observation guarantees the correspondence of their convex liftability property.

Proposition 1: If a given partition X of a polyhedron satisfies Assumption 1, then partition X is convexly liftable if and only if there exists a hyperbox B dx ( ), satisfying Assumption 3 such that the partition X ∩B dx ( ) = i∈I N (X i ∩ B dx ( )) is convexly liftable.

Proof: -→ It is already known ( [START_REF] Rybnikov | Polyhedral partitions and stresses[END_REF], [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF]) that the convex liftability of partition X leads to the existence of a strictly positive function s(•) defined over the (d x -1)-faces of partition X such that every internal (d x -2)-face is in equilibrium. Consider an internal (d x -2)-face of partition X denoted by F , we write F dx-1 (F ) to denote the set of (d x -1)-faces of partition X , sharing a common facet F , then the following inclusion holds true:

Ci∈F dx -1 (F ) s(C i ) n(F, C i ) = 0.
By the construction of B dx ( ), if F is an internal (d x -2)-face of partition X , then F ∩ B dx ( ) is also an internal (d x -2)-face of partition X ∩ B dx ( ). Furthermore, F ∩ B dx ( ) also satisfies the following property:

Ci∈F dx -1 (F ) s(C i )n(F ∩ B dx ( ), C i ∩ B dx ( )) = 0, due to the fact that n(F ∩ B dx ( ), C i ∩ B dx ( )) = n(F, C i ).
Therefore any internal face F ∩B dx ( ) is in equilibrium with strictly positive d x -stress: s(C i ∩ B dx ( )) = s(C i ) > 0, meaning the cell complex X ∩ B dx ( ) is convexly liftable. ←-The sufficient condition can be similarly proved. The goal is to present a construction of convex liftings for the given cell complex X of a polyhedron, based on an intermediate polytopic partition X ∩ B dx ( ). The following observation shows this close relation.

Proposition 2: Given a cell complex X of a polyhedron satisfying Assumption 2 and a hyperbox

B dx ( ) satisfying Assumption 3, if f : X ∩ B dx ( ) → R i.e. f (x) = a i x + b i for x ∈ X i ∩ B dx ( ),
is a convex lifting of the cell complex X ∩ B dx ( ), then the function g : X → R, defined as follows:

g(x) = a i x + b i for x ∈ X i ,
is also a convex lifting of partition X .

Proof: First, we observe that due to Assumption 3, the constraints characterized by B dx ( ) do not contribute to the division of the partition X ∩ B dx ( ).

From the hypothesis, f (x) is a convex lifting of partition X ∩ B dx ( ), then it can be deduced from two neighboring regions (X i ∩ B dx ( ), X j ∩ B dx ( )) that:

a i x + b i = a j x + b j ∀x ∈ X i ∩ X j ∩ B dx ( ), a i x + b i > a j x + b j ∀x ∈ (X i ∩ B dx ( ))\(X j ∩ B dx ( )).
Note also that the constraint a i x + b i = a j x + b j describes the set of points which lie on the hyperplane separating X i ∩B dx ( ) and X j ∩B dx ( ), then it also separates X i and X j . Therefore, the following inclusions hold true:

a i x + b i = a j x + b j ∀x ∈ X i ∩ X j , a i x + b i > a j x + b j ∀x ∈ X i \X j .
The same argument holds true for all pairs of neighboring regions in partition X ∩ B dx ( ), leading to:

g(x) = a i x + b i > a j x + b j , ∀x ∈ X i \X j , ∀j = i.
In other words, as per Definition 3, g(x) is a convex lifting of the cell complex X . The above results are meaningful to construct a convex lifting for the cell complex X of a polyhedron, based on a bounded one X ∩ B dx ( ). The structure of the latter partition is compatible with the algorithms presented in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF].

As previously shown, the construction of a convex lifting for partition X of a polyhedron can be carried out via the restriction of X into a large enough hyperbox in the same dimensional space. A natural question arises how this hyperbox can be chosen. This choice must clearly satisfy Assumption 3, and can be found through Algorithm 1.

Algorithm 1 Computation of a hyperbox B dx ( ) Input: A cell complex X of a polyhedron in R dx and a scalar a > 0. Output: .

1: V = i∈I N V(X i ). 2: Solve the following problem min s.t. ≥ 0, -( -a)1 dx ≤ x ≤ ( -a)1 dx , ∀x ∈ V. (4)
Note that such a constant a > 0 needs to be inserted in the set of constraints for guaranteeing Assumption 3 i.e. all vertices v ∈ i∈I N V(X i ) lie in the interior of B dx ( ). This constant can be freely chosen as long as it is strictly positive.

Based on these results, to state a solution to inverse parametric linear/quadratic programming problem via convex liftings for unbounded partition, we need to define some supplementary notation.

First, given a cell complex of a polyhedron X = i∈I N X i ⊆ R dx satisfying Assumption 2, a continuous PWA function f pwa (•) : X → R du is defined as

f pwa (x) = f i x + g i for x ∈ X i .
(5)

By (x), we denote a convex lifting of the cell complex X i.e. (x) = a i x+b i for x ∈ X i . Based on the above elements, the following sets are also defined:

V x = i∈I N V(X i ), R x = i∈I N R(X i ), V [x T z u T ] T = x T (x) f T pwa (x) T | x ∈ V x , R [x T z u T ] T =        r (r) f (r)   r ∈ R x , (r) = a i r f (r) = f i r if r ∈ R(X i )      , Π v = conv(V [x T z u T ] T ), Π r = cone(R [x T z u T ] T ), Π = Π v ⊕ Π r . (6)
With these notation, the following theorem presents a solution to IPL/QP problems via convex liftings for polyhedral partitions of a polyhedron. Theorem 3.1: Given any continuous PWA function f pwa (x) (5) defined over a cell complex of a polyhedron satisfying Assumption 2 and the sets defined in (6), then f pwa (x) is the image via the orthogonal projection of the optimal solution to the following optimization problem:

min [z u T ] T z s.t. x T z u T T ∈ Π. (7)
Proof: Consider x ∈ X i , for a region X i of the given partition X , according to the Minkowski-Weyl theorem for polyhedra (Corollary 7.1b in [START_REF] Schrijver | Theory of linear and integer programming[END_REF]), x can be described as follows:

x = v∈V(Xi) α(v)v + r∈R(Xi) β(r)r,
where α(v), β(r) ∈ R + and v∈V(Xi) α(v) = 1. As a consequence, the convex lifting at x, i.e. (x) can be described by:

(x) = a i x + b i = a i ( v∈V(Xi) α(v)v + r∈R(Xi) β(r)r) + b i , = v∈V(Xi) α(v)(a i v + b i ) + r∈R(Xi) β(r)(a i r).
Similarly,

f pwa (x) = f i x + g i = f i ( v∈V(Xi) α(v)v + r∈R(Xi) β(r)r) + g i , = v∈V(Xi) α(v)(f i v + g i ) + r∈R(Xi) β(r)(f i r).
It can be observed that if r is an extreme ray of X i , then r T a i r T is also an extreme ray of the affinely equivalent polyhedron Π [x T z] T of cell complex X defined as follows:

Π [x T z] T = conv(V [x T z] T ) ⊕ cone(R [x T z] T ),
where

V [x T z] T = x (x) | x ∈ V x , R [x T z] T = r (r) | r ∈ R x , (r) = a i r if r ∈ R(X i ) .
According to the definition of an affinely equivalent polyhedron, the optimal solution to the following problem:

min z z s.t. x T z T ∈ Π [x T z] T
falls in the lower facets of Π [x T z] T described by the convex lifting (x) of partition X . Now, consider a region X i , by

F (i) [x T z] T we denote the lower facet of Π [x T z] T such that Proj R dx F (i) [x T z] T = X i .
According to Proposition 5.1 in [START_REF] Nguyen | Inverse parametric convex programming problems via convex liftings[END_REF], every augmented point in V [x T z u T ] T corresponds to the vertices of Π v after lifting onto R dx+du+1 . Therefore, by preserving lifting onto R dx+du+1 , we guarantee the existence of a d x -face of Π, denoted by

F (i) [x T z u T ] T such that: Proj R dx+1 F (i) [x T z u T ] T = F (i) [x T z] T .
Note that such a d x -face is defined as follows:

F (i) [x T z u T ] T = F (i) 1 ⊕ F (i) 2 F (i) 1 = conv v T (v) f T pwa (v) T | v ∈ V(X i ) F (i) 2 = cone      r (r) f (r)   r ∈ R(X i ), (r) = a i r f (r) = f i r    .
It can be observed that at a point x ∈ X i , the minimal value of z can be obtained when x T z u T T lies in

F (i) [x T z u T ] T since Proj R dx +1 F (i) [x T z u T ] T = F (i) [x T z]
T and the optimal solution to:

min z z s.t. x T z T ∈ Π [x T z] T , falls at the point x T (x) T ∈ F (i) [x T z] T .
Therefore, optimal solution to (7) at x can be described by:

  x z * (x) u * (x)   = v∈V(Xi) α(v)   v a i v + b i f i v + g i   + r∈R(Xi) β(r)   r a i r f i r   , where α(v), β(r) ∈ R + and v∈V(Xi) α(v) = 1. It is clear that z * (x) u * (x) = a i x + b i f i x + g i = (x) f pwa (x)
, for x ∈ X i .

To complete the proof, we now need to show that the optimal solution to (7) is unique. In fact, at a point x ∈ X i , suppose there exist two different d x -faces denoted by F 1 , F 2 to which optimal solutions to (7) belong. Let

x T z 1 (x) u T 1 (x) T ∈ F 1 and x T z 2 (x) u T 2 (x) T ∈ F 2 , then it can be observed that z 1 (x) = z 2 (x), leading to Proj R dx +1 F 1 = Proj R dx +1 F 2 = F (i) [x T z] T . If u 1 (x) = u 2 (x), then F 1 , F 2 lie in a hyperplane of dimension d x + 1 which is orthogonal to the space of x T z
T . In this case, either f pwa (v) or f (r) or both are not uniquely defined for some vertices v ∈ V(X i ), some extreme rays r ∈ R(X i ). This contradicts the construction in [START_REF] Gulan | Implications of inverse parametric optimization in model predictive control[END_REF]. Therefore, F 1 = F 2 leading to the uniqueness of the optimal solution to [START_REF] Hempel | Every continuous piecewise affine function can be obtained by solving a parametric linear program[END_REF].

Based on the above result, Algorithm 2 summarizes a procedure to solve the inverse optimality problem to recover a continuous PWA function defined over a convexly liftable cell complex of a polyhedron.

Algorithm 2 Solution to IPL/QP problem via convex lifting. Input: A continuous PWA function f pwa (x) [START_REF] Grünbaum | Convex polytopes[END_REF]. Output: H x , H u , H z , K, J(x, z, u).

1: Call Algorithm 1 to find an appropriate B dx ( ). 2: Construct a convex lifting denoted as (x) of the cell complex X ∩ B dx ( ). 1 3: Compute (x) for the cell complex X . 4: Construct constraint set Π as in (6). 5: Solve the following optimization problem:

z * u * = arg min [z u T ] T z s.t. x T z u T T ∈ Π.
6: Project the optimal solution onto R du i.e.

u * = Proj R du z * u * .
Note that this result can be trivially extended for removal of redundant constraints as shown in [START_REF]On the complexity of the convex liftings-based solution to inverse parametric convex programming problems[END_REF]. Similarly, extensions of inverse optimality with quadratic cost functions and discontinuous PWA functions are referred to [START_REF] Nguyen | Any discontinuous pwa function is optimal solution to a parametric linear programming problem[END_REF].

IV. NUMERICAL EXAMPLES This section considers some numerical examples to illustrate the above results. All numerical examples in this paper are carried out in the environment of MPT 3.0 [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF].

Example 1

This subsection considers a continuous PWA function defined over the whole space R, as shown in Fig. 1. This partition is the union of several segments:

X 1 = {x ∈ R | x ≤ -0.5} , X i = {x ∈ R | 0.1(i -7) ≤ x ≤ 0.1(i -6)} , i = 2...11, X 12 = {x ∈ R | x ≥ 0.5} .
A convex lifting of this partition is presented in Fig. 2. A set of constraints for the formulation of inverse optimality problem, is shown in Fig. 3 as the shaded pink polyhedron. It is shown therein that the line along the x-axis, represents the given partition X , while the green curve represents the given PWA function associated with this partition. It can be observed that the solid pink line depicts the optimal solution to (7) whose image via the orthogonal projection onto R du coincides with the given PWA function. 1 An algorithm to construct a convex lifting of a bounded cell complex is presented in [START_REF]On the lifting problems and their connections with piecewise affine control law design[END_REF].

Example 2

The present subsection aims to illustrate the constructive solution with another continuous PWA function defined over a convexly liftable partition covering the whole space R 2 as seen in Fig. 4. This continuous PWA function is depicted in Fig. 5. A convex lifting of the given partition can be found in Fig. 6. Finally, the recovered PWA function is shown in Fig. 7, and can be seen to be similar to the original one.

V. CONCLUSIONS

The present paper shows a procedure to recover a continuous PWA function defined over a polyhedral partition of unbounded polyhedra. This solution is a simple extension of the constructive solution based on convex liftings put forward in previous results. Two numerical examples are presented to illustrate the constructive procedure proposed above.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Given PWA function defined over a partition of R.

Fig. 4 .

 4 Fig. 4. A partition covering the whole space R 2 .

Fig. 6 .Fig. 7 .

 67 Fig. 6. A convex lifting of the partition in Fig.4.


Zuosheng Lei

Bin Li

Rujun Wei

Yueming Zhou

Yunbo Zhong

Zhongming Ren

Two indexes to evaluate the flow field in slab continuous casting mold with electromagnetic stirring

Keywords: Slab Continuous Casting, Flow Control, Electromagnetic Stirring

Electromagnetic technologies, such as electromagnetic brake and electromagnetic stirring, are two efficient technologies to control the flow field in slab continuous casting mould. However, up to now there is still an important open problem to be answered firstly. It is that what kind of flow structure is our best choice to get good casting products when we use electromagnetic technologies ? We establish an physical model and an numerical model to study the flow field in slab continuous casting mold with electromagnetic stirring. During our physical experiments, we use liquid mercury as the analogue of liquid steel. In this study, we study the flow field under different casting parameters and stirring currents, respectively. At last, in order to evaluate the flow field in the mold, we propose two parameters. They are Mould Flux Entrapment Index and Velocity Uniform Index. Two parameters, named , are proposed in order to . Based on these two indexes, some optimized stirring parameters are proposed under different casting situations. These kinds of parameters are some reasonable criterions in order to find the optimized casting parameters when we use Electromagnetic flow control technologies.

Introduction

Fluid flow in the mold is very important to product quality in the slab continuous casting of steel [1] . There are two kinds of technologies to control the flow field in the CC mold. One is the submerged entry nozzle technology in that the submerge depth or the mechanical structure of SEN can be changed in order to modify the jet flow throughout the SEN exit port. The other one is electromagnetic technology in that Lorentz force can be generated in the up zone of liquid steel to realize the aim of controlling the flow field in the mold. Electromagnetic brake and electromagnetic stirring [2] are two typical electromagnetic technologies widely applied in the industrial process. Due to the high cost of empirical investigation by test-error type study on the relationship between the flow field in the mold and the CC slab quality, recently, mathematical modeling [1,3] and physical modeling using low melt point metal [4,5] have been developed to understand the fluid flow phenomena with electromagnetic field and we have achieve many valuable quantitative results. But up to now, most of the work about the flow field in mold can only tell us what it is rather than what should be. There is still an important open problem to be answered firstly.It is what kind of flow structure is our best choice to get the best casting products when we use electromagnetic technologies? As the complication of the multi-physical transport phenomena in the CC mold, it is indeed very difficult to answer the question. How to evaluate the flow field with the addition of electromagnetic field in the mold is a key problem when the standing position is improving the initial solidifying slab quality in the up part of the mold. In this study, two parameters, named Mold Flux Entrapment Index and Velocity Uniform Index, are proposed in order to evaluate the flow field in the mold. Based on these two indexes,some optimized stirring parameters are proposed for different casting situations. These kinds of parameters are some reasonable criterion in order to find the optimized casting parameters when we use Electromagnetic flow control technologies.

Physical modeling and mathematical modeling

A physical model using mercury is set up [START_REF] Zhou | EPM2015[END_REF] , whose photo and dimension of main parts are showed in Fig. 1(a) and(b). In this system, mercury is circulated by a electromagnetic pump from the mold bottom to the tundish and then flowing back to the up zone of the mold via the SEN. The vertical and horizontal velocity in the mold is measured by means of the ultrasonic Doppler velocimetry (UDV) in the measurement area showed in Fig. 1(b). The 1/5 scale mold with section size of 430mm × 46mm and 1000mm height is designed. And it is based on the same governing nondimensional Reynolds number and magnetic interaction parameter. The EMS Stirrer is a 6 poles linear induction motors with iron core height 80mm and width 450mm, and the maximum magnetic flux density is about 120 mT when the current is 60A. According to the same dimension of the physical model, a mathematical model using commercial packages Ansys® is built(Fig. 1 (c)). Fig. 1 (d) and (e) show the electromagnetic field in the mold on the center level of EMS Stirrer iron core. The governing Maxwell equations, N-S equations and boundary conditions will not be given here because these are not the critical points of this paper. Compared with the physical measurement flow field, the mathematical model is considered reasonable. The relevant results with or without EMS are showed in Fig. 2.

Valuation of the flow field in the mold

Then the validated mathematical model is used to simulate the flow field in the industrial mold under different EMS parameters and continuous casting conditions. Fig. 3 show some calculated results when the casting speed, the section

Where c v can be obtained by eq.

    1 / 2 2 1 2 1 2 1 2 2 c v               (2)
Where ρ 1 and ρ 2 are the density of liquid steel and mold flux, and  is the surface tension between liquid steel and mold flux. Eq.( 2) is the Kelvin-Helmholtz instability criterion between two immiscible flow liquid. It suggests that once the velocity difference of liquid steel and mold flux is greater than the critical velocity, the possibility of mold flux entrapment will increase. Although mathematic modeling describe above the flow field of mold flux is not calculated, an reasonable hypothesis is that the critical velocity should be 0.8m/s (Which can be determined by experiments in facts). Fig. 4 shows the free surface flow field and contour map of mold flux entrapment index under different EMS currents and casting speed. From the map, the region right to the red line is relatively dangerous for the viewpoint of mold flux entrapment. So the maximum EMS current should be selected under different casting speed.

As we know that the growth speed of initial solidifying shell is decided by two aspects. One is the cooling density from the water-cooled copper mold, and the other is the impact of the hot liquid steel flow in front of the initial shell. In order to prevent the crack formation caused by thermal stress, it is very important to keep the shell growing uniformly. It is the flow velocity in front of the initial solidifying shell that determine the convection heat transfer between the interface of the shell and liquid steel, so the uniform of the velocity in front of this interface is applied to valuate the impact of flow field on initial shell growth.

Due to the help of EMS Lorentz force, a horizontal circulation will form in the stirring zone, and it will affect the flow structure of the low part in the mold. A index named, Velocity Uniform Index (VUI), is introduced followed ref. [START_REF] Okazawa | [END_REF] to give a quantitative valuation of the effect of flow field on solidification on every level, which is defined as:

      100% m a x V x d x V U I V x d x         ( 3
)
Where Γ is a closed curve 20 mm away from the inner mold wall, see in the Fig. 5, and

x V is the absolute velocity of every point on Γ . It is obvious that VUI has a value between 0 and 1. This value is greater, better the flow field.Fig. 5 shows the VUI of different levels in the mold when the EMS currents are different. It can be seen that the VUI increases with the EMS current, but decreases a little in the deep part of the mold. Through the comprehensive index of two effects, an optimized EMS current can be achieved under a given Continuous casting condition.

Conclusions

In this study, the flow field in slab continuous casting mold with EMS is measured by physical modeling and calculated by numerical simulation under different casting parameters and stirring currents. Two parameters, named Mold Flux Entrapment Index and Velocity Uniform Index, are proposed in order to evaluate the flow field in the mold. Based on these two indexes, some optimized stirring parameters are proposed under different casting situation.

Fig. 1 Fig. 2

 12 Fig.1 Physical modeling and mathematical modeling System (a) Experimental system. (b) Schematic of experimental system. (c) mathematical modeling. (d) Calculate result of electromagnetic field when phase angle=0. (e) Current density in the mold.

Fig. 3 Fig. 4

 34 Fig.3 Calculate flow field in the mold with different EMS current (a) 0A, (b) 300A, (c) 400A, (d) 500A, (e) 600A, (f) 700A

UniformFig. 5

 5 Fig.5 VUI under different EMS currents

Acknowledgment

This project financially supported by National Science Foundation of China (NO. 50874133 and NO. 51274137).


R Ebert

C Resagk

New developments on Lorentz force velocimetry for weakly conducting fluids

Keywords: Lorentz force velocimetry, weakly conducting fluids, force measurement, non-invasive

. At applications for weakly conducting fluids such as electrolytes with conductivities in the range of 1 -10 S/m the challenging bottleneck is the detection of the tiny Lorentz forces in the noisy environment of the test channel. For the force measurement a state-of-the-art electromagnetic force compensation balance is used. Due to this device the mass of the Lorentz force generating magnets is limited. For enabling larger magnet systems and for higher force signals we have developed and tested a buoyancy based weight force compensation method which will be presented here.

Introduction

As it is well known from liquid metal MHD, characterized by relatively small magnetic Reynolds numbers [START_REF] Moffatt | Reflections on Magnetohydrodynamics[END_REF], Ohm's law for a moving liquid can be written in the following form:

𝒋 = 𝝈 (-𝛁Ф + 𝐯 × 𝑩 𝟎 ) (1)
Here σ is the electrical conductivity, denotes electric current density, Ф is the electric scalar potential, v is the velocity of the melt, and B 0 denotes the magnetic field. Ohm's law expresses the fact that in electrically conducting liquid, eddy currents are induced by the movement of the liquid through a magnetic field and the action of a curl-free electrical field, defined by the gradient of its scalar potential. Due to the solenoidal nature of B 0 , charge conservation and the mutual interaction of the eddy currents and the applied magnetic field, Lorentz force is generated, which density f defined by the relation f = j × B 0 . We now insert the characteristic values v, B, and L for velocity, magnetic induction, and length, respectively, and, upon integrating over the volume affected by the magnetic field we find the following important scaling relation for the induced Lorentz force.

𝐅~𝛔𝐯𝐁 𝟐 𝐋 𝟑 (2)
According to formula (2), the Lorentz force linearly depends on the product of electrical conductivity and conductor's velocity, hence measuring of the force under predefined conditions allows an estimation of the velocity, as it is described by [START_REF] Shercliff | The theory of electromagnetic flow-measurement[END_REF].

Experimental setup

The investigations are made at a closed loop channel (see fig. 2) where tap water is used as a working fluid. Its electrical Fig. 1: Principle of Lorentz force velocimetry: The fluid flow gets interspersed by a static magnetic field and hence, eddy currents with a secondary magnetic field are generated. According to the laws of electrodynamics, the resulting force on the source of the primary magnetic field is an accelerating force in the direction of flow. This force is called Lorentz force.

conductivity can be set easily by adding sodium chloride. Hence, the conductivity can be adjusted from 0.03 S/m up to 20 S/m. A centrifugal pump is driving the flow and enables average fluid velocities up to 4 m/s in the test section. The nozzle and a set of mesh and honeycomb filters provide a piston shaped velocity profile at the nozzle's outlet. To decouple vibrations from the pump the force measurement setup is placed at a separate fundament. Further, a fluid temperature control system keeps the temperature of the fluid constant with a peak to peak deviation of ± 25 mK. Hence, the electrical conductivity that varies about 2 %/K at room temperature can be kept constant sufficiently. As a velocity reference a commercial magneto-inductive-flowmeter is used which enables the calculation of a calibration factor between Lorentz force und fluid velocity respectively volume flow.

A state-of-the-art electromagnetic force compensation (EFC) balance is used in a 90° rotated orientation for measuring the Lorentz forces of the horizontal flow.

On the one hand the EFC balance allows high force resolution but on the other hand the mass of the magnet system that is hanged on its load bar is limited to approximately 1 kg to avoid breaking of its thin and very fragile solid state joints.

Even with rare-earth high-field magnet systems and numerical optimizations of the geometry of the magnets there unveils an upper limit in available Lorentz force and thereby a lower limit for the electrical conductivity of the fluid and generally a limited measurement resolution due to this mass restriction. It has been decided to test the compensation by buoyancy to use the above mentioned advantages of this concept: float bodies at both sides of the magnet system are placed in two tanks (à 380 x 110 x 100 mm^3) which are filled with a fluid (water and oil). The float bodies are mounted by side arms to the magnet system. For testing the setup in operating conditions of larger magnet systems, the mass of the magnet system has been increased from 1.008 kg to 3.01 kg by gluing 4 cuboids from flint glass (à 500 g, 77 x 46 x 41 mm^3) symmetrically so that the center of mass of the magnet system stays unvaried. In terms of magnetic field those glass bodies are only dead mass and they do not contribute to the generation of Lorentz force. This important condition needs to be considered in the evaluation later. The choice of dead mass instead of additional magnets have been made to allow a direct comparison of the noise levels with and without the new parts because here the amplitude of the Lorentz force is not changed in general.

The lift force of the float bodies has to compensate the additional weight of the glass cuboids of about 20 N and of course the weight of themselves and of the side arms. Water and paraffin oil have been selected and compared as buoyancy media in the tanks because they have a low electrical conductivity to avoid interaction with the Lorentz force measurements.

Especially the last-mentioned substance has a relatively high viscosity and a moderate heat conduction coefficient and hence temperatures differences will lead only to low thermal convection.

Overcoming limitation for the size of the magnet systems: Buoyancy compensated magnet systems For applying the LFV for weakly conducting fluids such as electrolytes with an electrical conductivity of 10 S/m and less the force measurement is the main problem because the occurring forces are in the same range smaller as the conductivities are: six orders of magnitude. In a 50 mm x 50 mm duct and a fluid velocity from 0 to 3.5 m/s the force resolution has to be less than 1 µN in the rather noisy environment of the flow channel. Forces can only be resolved by an electromagnetic force compensation balance (EFC balance). This value equals 10^-7 of the weight force of the magnet system and hence the force measurement system is mounted on a separate footing. Due to this necessary type of balance the weight force of the magnet system hanging on it is limited to 10 N. Larger magnet systems cannot be used without a compensation of the vertical acting force. The idea of weight force decreasing by using buoyancy would enable larger magnet systems and higher resolution for force and hence higher velocity resolution. As a side-effect this force has only a component in the opposite direction of the gravitational force without any horizontal components that would disturb the force signal. Fig. 3 and fig. 4 show the setup that has been built for using this effect with large float bodies, which have been constructed for water and oil as buoyancy media and with respect to the stability criteria for swimming objects [START_REF] White | Fluid Mechanics[END_REF].

Measurement results

As a reference for the evaluation of the concept of compensating weight by buoyancy a "dry" measurement series with the conventional setup has been made without the glass dead mass and without the float bodies and side arms. The gained data can be compared with two measurement series for the compensation setup: one with water and one with paraffin oil as buoyancy media, each with the additional dead mass which gets compensated now. At all measurements the working fluid in the channel (saltwater) had a fixed electrical conductivity of 10 S/m and a temperature of 21 °C. At first, it can be seen obviously that the new compensation concept is suitable for detecting and resolving the occurring Lorentz forces for both buoyancy media. The concern of too much noise due to thermal convection in the float bodies or from other sources can be invalidated, as the data in (fig. 6-a) to (fig. 6-c) show clearly.

For the compensation with water the measurement data (fig. 6-b) show a strong drift of the ground state (zero fluid velocity) that could not be removed completely by detrending. By the analysis of the setup there has been a noticeable increase in the weight of the float bodies due to leakages coming from small cavities in the mesh structure of the 3D printing process. In comparison to that, the drift of paraffin oil (fig. 6-c) is more linear and can be removed as good as in the reference measurement (fig. 6-a). In general there is a drift superimposing the Lorentz force signal which comes from small tilting of the setup due to small temperature changes and other effects. 6-d). This means that the measurement series for water and paraffin oil have a fairly good accordance with the reference measurement. The noise that remains after detrending and averaging can be quantified by analyzing the standard deviation (STD) of the signals. By extending the setup as it is described here the noise level of the force signal stays unvaried after applying our standard data processing: Even if the noise level of the raw data increases from 1.5 % in the reference setup up to 7 % in the compensation setup for fluid velocities of 3 m/s the filtered data show a STD for all setups in the range below 1 %. Hence, no influence of using the principal of buoyancy compensated magnet systems to the noise level has been measured.

Conclusion and Outlook

We demonstrated that the concept of buoyancy compensation is a feasible method which enables large scale magnet systems for Lorentz force velocimetry at weakly conducting fluids. The initially unknown effects of the additional float bodies like damping of force amplitudes or higher noise level have been measured: There is no negative influence to the measurement resolution -on the contrary there is completely no amplitude modification and the additional noise can be eliminated by the yet used moving average filter. Thus, a comprehensively positive evaluation of the concept can be given. Replacing the dead mass that have been used for the simulation of a higher magnet weight force would result in an at least 3 times higher force signal and hence in a much better SNR. Therefore, the results shown here are only an underestimation of what would be possible with a real full size magnet system. Accordingly, this technique is able to improve the resolution of the Lorentz force Velocimetry at least by a factor of 3 too. In general, there is no limitation for using even larger magnet systems such as Halbach magnet arrays with 10 kg weight to get an increase of resolvable fluid velocity or a decrease of electrical conductivity of the fluid of a factor of 10. Here, 90 % of the magnet system's weight force needs to be compensated. Further, new technologies like superconducting magnet systems in cryostats which definitely exceed the mass restriction became possible applications for the Lorentz force velocimetry.

Fig. 2 :

 2 Fig. 2: Schematic illustration of the test channel with (a) pump , (b) nozzle, (c) test section / duct, (d) mounting for Lorentz force velocimeter, (e) mechanical decoupling and (f) relax vessel.

Fig. 3 :

 3 Fig. 3: Principle setup of buoyancy compensation. Fig. 4: Sketch of the casting flow problem with magnet pole faces and numerical grid for the simulations.

  For quantifying the results, the changes in Lorentz force have to be plotted over the related changes in fluid velocity by determining the heights of the flanks in (fig.6-a) to (fig.6-c). Results are shown in (fig.6-d). It can be seen clearly, that there is no damping or distortion of the force for both buoyancy media, because the amplitudes between the plateaus and the adjacent zero-flow phases are nearly the same for each velocity domain. The largest deviation occurs between the (a)-(c): Filtered (detrended and running average filter) Lorentz force signal with lifting body for different media in the tank. Velocities are measured with the magneto-inductive-flowmeter. (d): Lorentz force over related fluid velocity: Good accordance to the reference measurement qualifies the new principle for enhancing the measurement resolution of the Lorentz force velocimetry at least by a factor of 3. reference measurement and the measurement with water: 5.6 % (fig.

Acknowledgment

The authors are grateful to the German Science Foundation (Deutsche Forschungsgemeinschaft) for financial support of the presented work in the framework of the Research Training Group "Lorentz force velocimetry and Lorentz force eddy current testing" (GRK 1567/2) at Ilmenau University of Technology.


J Hodek

M Zemko

P Shykula

Finite Element Model of Gear Induction Hardening

Keywords: gear hardening, finite element, induction heating

This paper presents a finite element model of a gear induction hardening process. The gear was surface-heated by an induction coil and quickly cooled by spraying water. The finite element model was developed as a three-dimensional model. The electromagnetic field, temperature field, stress distribution and microstructure distribution were examined. Temperature and microstructural characteristics were measured and used. The gear material data was obtained in part by measurement and in part by calculation.

Introduction

Gears are machine parts that are widely used in industry. The end-use properties of selectively hardened parts of gears, in particular their hardness, wear resistance and contact fatigue resistance can be controlled to a great extent through the parameters of the heat treatment used. FEM simulations [1][2][START_REF] Šuchman | Optimization of Induction Hardening Parameters through FEM Simulation[END_REF][START_REF] Hodek | 23 th International Conference on Metallurgy and Materials, Fem Simulation Induction Heating Process[END_REF] and measurement of material properties, temperatures and martensite layer thickness are useful tools for describing the properties of gears upon heat treatment.

Process Description

The gear teeth are induction-heated (Fig. 1) to the required heat treatment temperature and then cooled rapidly. By varying the heat treatment temperature, the inductor power and the cooling method, one can obtain various mechanical properties upon heat treatment. The power and frequency of the inductor used in the present process was 100 kW and 200 kHz, respectively. The heating time was 0.8 seconds. The diameter of the gear was 100 mm and the module was 1.64 mm.

Fig. 1 Induction heating of the gear Material Data Acquisition

The L75PT dilatometer was used for the thermal expansion measurement. From the thermal expansion profile, the austenitizing temperature was determined. The LFA1000 laser flash instrument was used for the determination of heattransport properties: the specific heat and thermal diffusivity as a function of temperature. Thermal conductivity was calculated from thermal diffusivity and measured density values. The rest of the properties required were calculated using the JMatPro software. JMatPro is simulation software which calculates a wide range of material properties from the chemical composition. It focuses particularly on multicomponent alloys used in industrial practice. The following temperature-dependent data were computed: electrical conductance, Young's modulus, Poisson's ratio and the flow stress for all components of the material's mixture.

Measurement of Temperature

Temperature profiles during heating and cooling were measured for several combinations of power and time. Thermocouples were used for the measurement (see Fig. 2). At every instant, four temperatures were measured simultaneously: at two points at the top of the tooth and at two points at its bottom. Average values were then calculated from the temperatures of the tooth top and bottom. Those values should be understood to be approximations because it was not possible to weld the thermocouples onto exactly same locations. During the temperature measurement, the gear was not revolving.

Martensitic layer after heat treatment Martensitic Layer Thickness

The thickness of the martensitic layer was measured in the centre of the tooth. The shape of the martensitic layer is shown in Fig. 3.

Mathematical Model

The mathematical model is defined by the following equations: 1) The electromagnetic field is described by equations for the vector A and scalar V potentials [START_REF] Marc | [END_REF]. For the sake of simplicity, the problem is considered to be a harmonic one, where μ stands for the magnetic permeability, ω is angular frequency and γ denotes electrical conductance:

(1)

(2) 2) The transient temperature field is expressed by the Fourier equation [START_REF]System Documentation[END_REF]:

(3) Here, c denotes specific heat, ρ stands for density, λ is thermal conductivity and w represents the heat loss during induction heating.

3) Stress and transformation field

In the DEFORM 3D software, strain, heat transfer and transformation are strongly related parameters. Any change in one of them has impact on the other two. In the model, a material is defined as a mixture of phases where rules for phase transformations apply [START_REF]System Documentation[END_REF]. The stress strain relation is defined as an elastic-plastic relation for all phases:

(4) Here, σ denotes the flow stress, c represents a material constant, ε is the effective plastic strain, n stands for the strain exponent, m is the strain rate exponent and y represents the initial yield value. The initial phase is pearlite. The model contains two phase transformation rules: a. Pearliteaustenite: diffusional transformation described by the Avrami equation. b.

Austenitemartensite: martensitic transformation characterised by the martensite start temperature and the 50%-martensite temperature. Owing to high cooling rates used in the process, other phase transformations only take place on a very limited scale. Hence, they are omitted from this analysis.

FE Model of Heat Treatment Procedure

The block diagram of the calculation is shown in Fig. 4. The entire calculation procedure consisted of two parts. In the first one, the electromagnetic and temperature fields were studied using the MSC Marc software. The data for the temperature field just after the end of heating was then fed into the DEFORM 3D software by means of the Python script. The cooling process, the stress distribution and transformation of microstructure constituents were then analysed using DEFORM 3D.

Fig. 4 FEM model

As the problem is geometrically symmetrical, the computational model only comprised one quarter of a single tooth (see Fig. 5). The models of the electromagnetic and temperature fields in the MSC Marc software comprised the gear, the inductor and a sufficiently large portion of the environment. The calculation of stresses, temperatures and transformations in DEFORM 3D only involved the gear. Relevant boundary conditions were defined for individual mathematical expressions. The model was meshed to comprise at least three elements across the skin depth (see Fig. 6) Fig. 5 FE modelone quarter of the gear tooth Fig. 6 The meshed

Model Calibration

The model was calibrated using the measured data in order to obtain agreement with the real-world processes that take place during heat treatment. 1) Temperature calibration: the inductor current was varied until the measured and calculated temperatures were in agreement (see Fig. 7). During cooling by the water spray, heat was removed rapidly from the gear. The heat transfer coefficient value for this cooling sequence was determined by inverse analysis. 2) Calibration against the martensitic layer depth: the diffusional transformation model and the martensitic transformation model were varied until the measured and calculated martensitic layer data were in agreement (see Fig. 8).

Results and Discussion

Several heating and cooling possibilities were computed. The following computation results were obtained for 0.8second heating and 8-second water cooling. The temperature of the tooth upon heating (see Fig. 9) is between 1000 and 1200 °C. The region where the temperature is sufficiently high becomes austenitized (see Fig. 10).

Conclusion

The above-described FE model characterises the behaviour of gears during surface hardening. The model integrates calculations of electromagnetic, temperature and strain parameters and accounts for phase transformations. It was calibrated against measured temperatures and the thickness of the martensitic layer. Temperature-dependent material properties were in part measured and in part calculated. Using the FE model, one can estimate the behaviour of the gear during heat treatment, namely the time profiles of temperature, phase and stress distribution.

In their follow-up research, the authors will focus on the relationship between the residual stress distribution and the contact fatigue resistance of gears.

Fig. 2

 2 Fig. 2 Temperature measurement Fig. 3 Martensitic layer after heat treatment Martensitic Layer Thickness The thickness of the martensitic layer was measured in the centre of the tooth. The shape of the martensitic layer is shown in Fig. 3.

Fig. 7

 7 Fig. 7 Comparison between calculated and measured temperatures

Fig. 9

 9 Fig. 9 Temperature distribution after heating Fig. 10 Austenite distribution after heating

Acknowledgment

These results were created by project Development of West-Bohemian Centre of Materials and Metallurgy No.: LO1412, financed by the MEYS of the Czech Republic.


V Nemkov
email: vsnemkov@gmail.com

How Accurate is Computer Simulation of Induction Systems?

Keywords: induction heating, computer simulation, modeling, accuracy, precision

Computer simulation is a widely accepted tool for design of induction coils and other components of induction systems, for development of new processes and equipment, for troubleshooting and for teaching and learning. Multiple programs and packages are being used for these purposes, from "home-made" programs to multiphysics packages such as Comsol, Ansys or Flux. One of the first questions, asked by new users, is: "How accurate are the results?". The author has long experience in computer simulation and design of various induction systems and understands well that there is no simple answer to this question. This presentation is an attempt to discuss the sources of inaccuracies and provide information on what the user must pay attention to. It may be useful for the software users and developers as well as for all users of induction heating. The study is related to induction heating itself, i.e. to electromagnetic and thermal processes with much less attention to other processes. Simulation of structural transformations, stresses and deformations as well as electromagnetic forces and magnetohydrodynamic processes are not considered in this presentation.

Introduction

Induction heating is used in many technologies that usually include multiple physical processes and material transformation. The core processes are electromagnetic and thermal which may be followed by phase or structural transformation, additional mechanical processes (rolling, forging, pinning), etc. A scheme of processes in heat treating of machined parts is shown in Fig. 1. Majority of induction-oriented programs allow the user to simulate only few mutually coupled processes, e.g. electromagnetic (EM) and thermal (TH). The most effective even if more challenging, is computer simulation of the whole technological chain [START_REF] Muelbauer | History of induction heating and melting[END_REF][START_REF] Cai | Integration of induction heat treat simulation into manufacturing cycle[END_REF] using multiphysic programs or a set of programs. For example, program Flux 2D/3D allows simulating EM and TH processes with relatively simple supplying circuitry and possibility to move the coil position relative to the part (scanning with step-variable speed). It is obvious that multistage simulation may be less accurate if there is no intermediate control and corrections. Only the most traditional processes of induction heat treating are considered in this presentation, i.e. EM and thermal. A question about the simulation accuracy can't be answered in a simple way,because different criteria and norms are being used in different cases. There are specific difficulties in defining accuracy of simulation of the induction processes. One of them is a different type of compared objects. They may be numbers (such as power or heating time), continuous 1D, 2D and 3D functions, varying in time or space. The other problem is an uncertainty of experimental data and limited direct access to some values, such as internal temperatures. Other values can't be measured during processing, e.g. internal stresses. A problem of accuracy is closely linked to optimization of designs and processes, which requires a formulation of criteria and comparison of "nearness" of the results at different steps [START_REF] Ya | Optimalcontrol of induction heating processes[END_REF].

Analysis

Accuracy is the nearness of a calculation to the true value.Common definition associates accuracy with systematic errors and precision with random errors, typically distributed according to a normal distribution law. Foreach particular case and model, a simulation is precise in sense of stability of the results but not accurate. Experiments are not precise and not quite accurate. Comparison of calculated and measured data for a very simple case of the furnace impedance, presented in a form of a vector ReZ + jImZ, is given on Fig. 2. Black dots represent results of simulation using different programs or input data. It is clear that an error may be small or big, depending on the program used, calculation procedure and reliability of the experiment. For matching purposes it is sufficient to be sure that the calculated coil impedance is within arated load range of the power supply. For the user of induction machine, a criterion of the simulation accuracy may be a correspondence of predicted hardness pattern to the measured or specified pattern allowance.Inaccuracies of computer simulation can be caused by multiple factors and can appear at any step of simulation. Possible sources of inaccuracy for a typical chain of induction simulationare presented on Fig. 3.

Considered processesmay not cover all the processes occurring inside the induction system that can influence the results of simulation. For example EM +Thermal simulation of one-side heating of a flat body does not take into account the part deformation which can strongly influence the coil parameters and power distribution.Another importantfactor is a mode of energizing the induction system. The most common mode is an application of known current or currents to the system input (inputs). A program Elta [START_REF] Nemkov | Learning and teaching induction heating using programELTA[END_REF] in addition to traditional modes of current or voltage values on the coil or generator terminals has modes of known power or specific power on the part surface. The latter mode allows themodeler to design and optimize the process without detailed design of the coils. Elta has also a unique option of the heating process simulation with automatic adjusting of frequency, which is a common case in practice.

Fig. 3: Graph of possible inaccuracy sources

System Geometry is one of the most important sources of inaccuracy. The modellers always try to avoid 3D simulation, describing the system as 2D, in order to simplify their job and reduce processing and analysis time. However, it is necessary to state that almost all induction systems are 3D. Even simulation of traditional multi-turn solenoidal coils as a set of parallel plain turns can lead to significant errors when number of turns isn't big. This issue is a technical problem too. There are three types of multiturn coils: a) simple helical coils, Fig. 4; b) helical coils with "equilizers" (wedge-like pieces of copper brazed to extreme turns in order to make the ends parallel); c) coils with "plane-parallel" turnswith axial connections between them, often called "dog legs". It looks like the latter design is completely 2D, but dogleg connections increasean angular distance between the leads, which can significantly influence the coil parameters and the technological effects, especially inthe cold crucible melting furnaces. The "global" 3D effects in the cold crucible furnace can result in excessive skull thickness in the lead area and even cause unintended melt rotation. Helical and some other "supposed-to-be 2D" coils can induce longitudinal voltage in the part. If the part is longer than the coil and has connection to the "ground" (such as driving rolls, Fig. 5), it can result in significant induced current Ig thru the ground and in technological effects, which can't be simulated with 2D approach. These effects are usually negative because of undesirable heating, power loss, sparking and fast wearing of the roll bearings, but sometimes the "ground currents" can be used intentionally for additional selective heating. Phenomenon of "ground current" is known in practice but the author did not meet any publication aboutsimulation of this subject. Even single turn cylindrical coils have 3D zones in the lead area. Application of magnetic controllers (shunts, concentrators) can significantly reduce or eliminate 3D effects in single-and multiturn coils, Fig. 6. In general, neglecting of the lead zonein many cases can lead to big errors.

Material properties are another major source of inaccuracy especially in the case of magnetic materials. Temperature dependence of physical properties for majority of materials that are not intended to work at high temperatures, aren't easy available. For example,the author found that the best source of specific heat capacity of carbon steels in wide range of temperaturesis more than 80 years old but it is still much demanded [START_REF] Saburo' Umino | On the specific heat of carbon steels[END_REF]. For accurate simulation of steel heating we need to take into accountnot only chemical composition but also their structure (annealed, hardened, etc.). There is one more factor that typically isn't considered. It is well known that in forge heating the surface of steel oxidizes forming a scale layer. This reaction is very exothermic and ideally it must be taken into account.

Reliable and accurate description of magnetic properties is much more complicated. In some programs (Quick Field, Ansys, etc.) permeability is supposed to be constant. It can lead to an error around 30% in value of absorbed power for a given field strength on the part surface [START_REF] Slukhotskiy | Inductorsforinductionheating, Energia[END_REF][START_REF] Nemkov | Theory and calculation of induction heating devices[END_REF]. In majority of other programs used for induction heating simulation, "frequency domain" approach is being used [START_REF] Slukhotskiy | Inductorsforinductionheating, Energia[END_REF][START_REF] Lr | Skin effect in ferromagnetic bodies[END_REF]. In this case permeability is supposed to be constant during the field half period but its value is being calculated as a function of the field rms or amplitude. This approach can give "technically correct" power value but does not allow accounting for hysteresis or higher harmonics in "linear" parameters (field strength, flux density, current and voltage). Experimental tests and adjustment of simulation and design procedure continue until now [START_REF] Zedler | Investigation of relative magnetic permeability as input data for numerical simulation of induction surface hardening[END_REF]. "Time-domain" approach is based on simulation of EM "transient" processes during each period of the field wave. Though attempts to simulate eddy-currents problem in magnetic bodies have been made around 60 years ago [START_REF] Muelbauer | History of induction heating and melting[END_REF], the time-domain approach isn't used in common practice of induction heating. We can expect that with powerful modern computers we can return to this approach for more accurate simulation.

There is another "grey area" in simulation of induction heating of magnetics. It is a "stripe" or "zebra" effect. Under certain conditions we can observe periodical or quasi-periodical bright stripes on the part surface when its temperature passes thru the Curie point, Fig. 7. Thiseffect was described by Prof. G. Babat and other scientists in 1950-1960s [START_REF] Babat | Induction heating of metals and its industrial applications[END_REF].

Several explanations have been proposed for zebra effect with insufficient proofs. Only recently Dr.S. Dzliev published several papers on simulation of zebra effect [START_REF] Dzliev | Instability in induction heating of magnetic steels[END_REF] explaining it by "instability of permeability" close to Curie point, Tc. This study and simulation at Centre for Induction Technology using Flux 2D program (Fig. 8, right) showed that zebra effect can appears only when dependence of permeability upon temperature is close to rectangular with sharp drop near Tc. This dependence (Fig. 8,left) may be approximated by following expression:

Here μ rt is permeability at room temperature.Lower values of n correspond to higher permeability and therefore lower field strength.

𝜇 T = 1 + 𝜇 𝑟𝑡 -1 (1 -T T c 𝑛 ) (1)
Simulation experience showed that zebra effect takes place only when value of nin equation ( 1) is high, probably above 16. It explains why this effect appears more often in high frequency hardening when the field strength is much lower than at the low frequency hardening. To opinion of the author, zebra effect is caused by EM end effects on the plane of contact of magnetic and partially non-magnetic pieces of the part [START_REF] Nemkov | Theory and calculation of induction heating devices[END_REF]. More simulation tests with experimental validation are necessary for complete explanation and prediction of this effect. It is desirable to perform time-domain simulation.

Simulation type and method. It is difficult to compare accuracy of different simulation types and methods. Some considerations about Time Domain and Frequency Domain simulation of magnetic material processing had been presented above. By means of Time Domain simulation of heating magnetic materials, it is possible not only to obtain more detailed and accurate information about the electromagnetic processes but also validate and possibly improve the Frequency Domain method.

Algorithmof calculationsand post-processingcan significantly influence the program accuracy. Special attention must be paid to the algorithm and processing setting in simulation of strongly non-linear tasks such as impulse hardening orintensive quenching, otherwise the calculation can non converge or be completely incorrect. Post-processing is very important for evaluation of the results using selected criteria of accuracy.

Conclusions

-Computer simulation is a powerful tool for study and design of induction heating systems, but it isn't quite accurate.

-Incomplete or incorrect description of the system geometry and of a set of processes can result in big errors.

-There are still problems in simulation accuracy due to a scarce database of the temperature and field dependent properties of steels and due to insufficient study of the magnetic material characterization (account for hysteresis losses, "equivalent" permeability calculation in harmonic approach, etc.). -Accuracy depends on many factors including the problem formulation, correctness of the input data, used program and experience ofthe modelers. -It is very important that the modeler knew physics and practice of induction heating as well as assumptions and limits incorporated into the used software. -The modelershould make "credibility evaluation" at each step of simulation (model construction, accuracy of input data, distribution of calculated EM field, power and temperature) in order to improve accuracy and avoid mistakes. -Alternative simulation and experimental data can provide final or additional credibility to the obtained results.

-For each particular area of application it is useful to create a data or even knowledge base that can help to implement the best simulation practice. -Properly made simulation of induction heating may be effectively incorporated into the whole simulation and optimization chain of the product manufacturing [START_REF] Cai | Integration of induction heat treat simulation into manufacturing cycle[END_REF][START_REF] Ferguson | Modeling stress and distortions of full-float truck axle during induction hardening[END_REF].

Fig. 1 :

 1 Fig. 1: Mutually coupled processes of induction heating Fig. 2: Comparison of simulation and experimental data

Fig. 4 :

 4 Fig. 4: Typical helical induction coil Fig. 5: "Ground" eddy current in induction systemwith helical coil

Fig. 7 :Fig. 8 :

 78 Fig.7: Stripe effect: 477 kHz ( left) and 242 kHz (right). Photos: M. Divilkovsky[START_REF] Babat | Induction heating of metals and its industrial applications[END_REF]

Acknowledgment

Support from Fluxtrol, Inc., USA, and kind advices of colleagues from different countries are gratefully acknowledged.


J R Alves

J Barlier

S Marie

C Beraudo

F Bay

Modelling of large displacements and large strains in coupled electromagnetic/solid mechanics/heat transfer problems

Keywords: Electromagnetic Forming, Induction Heating, Finite Elements, Mesh adaptation. Influence zone Influence zone

Modelling of large displacements and large strains in coupled electromagnetic/solid mechanics/heat transfer problems

Introduction

Several EPM processes involve large displacements or large strains of the solid bodies. We can mention magnetic forming applications [START_REF] Psyk | Electromagnetic forming -A review[END_REF] where the part to be formed undergoes large strains and strain rates,induction heat treatment applications where the inductors may have a large motion with respect to the part to be heated/treated. This problem is one of the core subjects shown in this work.

We shall present here the work carried out on the coupled simulation package MATELEC/FORGE® [START_REF] Fourment | State-of-the-Art 3D simulation Tool as a Strategic Investment[END_REF], where MATELEC is the FORGE® electromagnetic module developed in our laboratory. This tool enables solving the electromagnetic equations by means of the finite elements method. The sequentially-coupled approach between MATELEC and the thermo-mechanical solver enables an optimized treatment of the physics, preserving reasonable computational cost for complex cases. The main industrial applications treated here are electromagnetic forming and induction heating.

A full finite elements approach has been retained for solving the electromagnetic equations. This means that besides the solid parts, a surrounding medium (air) has to be meshed in order to allow electromagnetic coupling between the conductors. The movement of the solids within the surrounding domain is included in the solver by means of automatic r-adaptation or remeshing techniques.

Simulation Strategy

One of the main challenges for solving multi-physical problems is how to set-up the resolution scheme in terms of the coupling between the physics. As shown in Fig. 1, all the physical aspects are highly dependent one on each other. The two main options that can be proposed are: (a) To include all the physics in a single solver for a simultaneous resolution (strongly-coupled). Or (b) solve each physics in a separate fashion, communicate and correct iteratively. We have chosen to follow this sequentially-coupled approach because solving the physics separately has the advantage of allowing independent management of the meshesmeaning that each problem will have a finite element space adapted to its own requirements. ) and the Joule heating term ( Q  ) at each element of the solid body from the EM module. The temperature evolution ( T ) and the velocity field ( v 

) or the new position ( x 

) of the nodes of the solid body are transmitted from the TM module. The main difference between the electromagnetic forming application and the induction heating process is that the communication is done after each time-step increment in the first case or after a full electric period for the latter. The electromagnetic wave propagation phenomenon is described by the Maxwell's equations. The latter is commonly expressed in potential form as shown in set (1) for many numerical applications.

               A t A r   1 1 0 (1.a)   0       (1.b)
Where, A  is the magnetic potential,  is the electric potential,  is the material electric conductivity,  the magnetic permeability,  the electrical permittivity, subscript " 0 " represents the value for void and " "r the relative value with respect to the void. The thermo-mechanical problem is described by the conservation laws. Equation ( 2) represents the conservation of energy with heat sources (induction heating) and the effect of self-heating by mechanical deformation. p C denotes the heat capacity, the main unknown T the temperature, and k for the heat conduction within the body. b dt

       :       vol p q T k dt dT C
v d         (3.a)  P v       (3.b)
Where  is the Cauchy stress tensor, b  represent any volumetric force such as gravity (which is not taken into account in our analysis) or the Lorentz body forces, in which case

B J F b lorentz        (4)
The inertia is described by the material density  and the material time derivative of the velocity v d  . The hydrostatic pressure is represented by P and the bulk modulus  defines the compressibility in the elastic deformation regime.

Parallelism

3D simulation of coupled electromagnetic/solid mechanics/heat transfer problems with large displacements and large strains leads to wide number of increments and can be CPU time intensive. Therefore use of parallel clusters machine is strongly recommended for the resolution of such long processes.

The coupled MATELEC/FORGE® simulation strategy makes it easy to handle two levels of parallelism: the parallelization of the FORGE® thermo-mechanical Finite Element software which is already efficiently done [START_REF] Coupez | Parallel simulation of 3D forming processes[END_REF] and the parallelization of the electromagnetic algorithm itself [START_REF] Micheli | Towards the simulation of the whole manufacturing chain processes with FORGE®[END_REF]. In both cases the parallelization method is similar, based on a SPMD (Single Program on Multiple Data) strategy where each processor runs a full version of the code on a partitioned mesh. Within FORGE® the SPMD paradigm is based on the tetrahedral elements partitioning of the mesh billet because the thermo-mechanical computation is carried out only on the work piece.

The electromagnetic formulation consists in discretizing in space of the whole domain in one single mesh (parts, inductors, air) using Nedelec edge elements [START_REF] Nédélec | A New Family of Mixed Finite Elements in R 3[END_REF]. Thus the magnetic vector potential/electric scalar formulation described in the previous paragraph leads finally to resolve symmetric linear systems with advantageous small bandwidth. It ensures a good performance on parallel linear iterative solvers which are solved via the external PETSc library [START_REF] Balay | PETSc Users Manual[END_REF]. The parallel Finite Element assembly is based on tetrahedral elements partitioning of the whole mesh domain. A nodes partitioning data structure is built with a specific renumbering where all the nodes of a parallel sub-domain are on the same processor. It enables to minimize communications between processors when solving linear systems. A parallel data structure partitioning by elements, nodes and edges is implemented to update and synchronize global vectors and scalars between parallel sub-domains.

The edges formulation of the magnetic potential induces a large number of degrees of freedom in comparison to the parallel communication inter processors which is proportional to the number of nodes of the subdomains interfaces. Therefore, the parallelization strategy provides a good parallel efficiency on a large number of cores as well on entry level parallel computers.

Models

The first model consists in a magnetic pulse welding simulation of an aluminum tube with a copper cylinder as shown in Fig. 3. In the EM module we take into account the inductor, the work piece and a surrounding air box in order to allow the propagation of the electromagnetic waves. Besides the EM computation, we also define a coupling with a RLC circuit model in which the electromagnetic forming machine parameters are defined. The electric current evolution is then computed as a result from the coupling between the circuit and the FE model. In the thermo-mechanical module the work piece is defined again together with any tools/holders that may be needed for contact conditions. Since the meshes are in principle different, the communication is carried out by interpolation of the data at the integration points. The second model consists on an induction heating process of a steel bar for heat treating the outer surface by means of an AC current of 2.1kAmps as shown in Fig. 4. For this process, the EM module computes a complete period of the AC signal and the Joule heating transmitted to FORGE® corresponds to the rms value on the period. Then, FORGE® computes the evolution of the temperature field until a defined variation before the coupling is reestablished.

Results

For the first case corresponding to the magnetic pulse welding application, we had small displacements of about 1.5 mm but large deformations on the work part reaching 33% before impact. This condition allows taking full advantage of the algorithm presented Fig. 2. The small incremental displacements enable to use r-adaptation algorithm to follow the geometry evolution as can be seen in Fig. 5. At this stage, several elements have in surrounding air have withstood critical deformation. This condition automatically launches the remeshing module at the next electromagnetic time step increment.

The thermo-mechanical model also counts with its automatic remeshing/mesh adaptation module. While several options are available, we show in this study how it is possible to setup a "box" region near the impact zone. This is useful to create a fine mesh zone just prior the contact condition as can be seen in Fig. 6. While in EMF large deformation with small to medium displacements are expected, in induction heating application we will most likely face large displacements (of the inductor and/or the treating piece) with small to negligible deformations. In this case, the r-adaptation algorithm is not adapted to capture the movement. In such case, only the remeshing module is used. Using a predefined parameter, an influence zone is defined from the inductor. Within this zone, the material parameters are used to automatically define the mesh size in order to respect the skin depth effect at the surface of the work piece as is shown in Fig. 7.

Conclusions

One of the complexities associated to multi-physical simulation of electromagnetic forming or induction heating processes it the careful management of the mesh generation and the continuous tracking of displacements and deformations. In this work we have presented the algorithmic procedure that was defined within the FORGE®/MATELEC interface in order to automatically generate and adapt the global electromagnetic mesh due to the displacement and/or deformation of the immersed solids.

Fig. 2

 2 Fig. 2 describes the global work flow to accomplish a single time step increment. The coupling between the solvers is driven by the communication of: The Lorentz volumetric forces ( F ) and the Joule heating term ( Q  ) at each element of the solid body from the EM module. The temperature evolution ( T ) and the velocity field ( v 

Fig. 1 :Fig. 2 :

 12 Fig. 1: Scheme of the physics coupling and interaction

( 2 )

 2 Equations (3) represent the conservation of linear momentum (3.a) and the conservation of mass for a compressible body (3.b).

Fig. 3 :

 3 Fig. 3: Magnetic pulse welding of an aluminum A2024 tube with a copper cylinder. Fig. 4: Induction heating case of an steel A304L billet.

Fig. 5 :

 5 Fig. 5: Mesh adaptation of the air around the deformed billet. (Left) R-adaptation until t=3.2µs. (Right) Automatic mesh generation at t=3.4µs.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Mesh refinement in the thermo-mechanical computation to capture impact conditions. (Up) Before impact. (Bottom) After/during impact.


E Sauvage

P Brun

A Bonnetier

J Lacombe

R Didierlaurent

Numerical simulations in the development of the French radioactive waste vitrification processes using induction furnace

Keywords: cold crucible, simulation, glass, fluid flow, induction

For many years, the CEA (Commissariat à l'Énergie Atomique et aux Énergies Alternatives) Marcoule France has developed various processes dedicated to radioactive waste confinement, especially vitrification processes for HLLW. For 15 years now, the numerical simulation has become an important tool for research and developement projects held in the CEA-AREVA Joint Vitrification Laboratory (LCV). Induction heating, fluid mechanics and thermal simulations take part of all new R&D projects. The apports of such simulations are, first, the enhancement of the working knowledge of existing process. Those data are very useful to define optimisation choices, for example upgrades made on the hot metallic melter used since the 90s at LaHague facility. Second, the simulations are, of course, also used at the conception stage of new processes as a tool allowing wide ranges parametric tests. This has been extensively used in the design of the cold crucible inductive melter (CCIM) commissioned in 2010 at La Hague plant. Finally, it is a powerful and relatively cheap tool for prospective studies for processes of the future. Whatever the purpose, the potential benefits are gains on the reliability, the output capacity and the life time.

Introduction

Vitrification of high-level liquid waste is the internationally recognized standard to both minimize the impact to the environment resulting from waste disposal and the volume of conditioned waste. In France, the vitrification of highlevel liquid waste produced from nuclear fuel reprocessing has been successfully operating now for more than 35 years with three major objectives: durable containment of the long-lived fission products, minimization of the final waste volume and operability in an industrial context. As a result, CEA & AREVA have integrated a unique experience in the field of high level waste vitrification through the design and operation of facilities with high records of safety, reliability and product quality, in line with efficient reprocessing plants; continuous efforts to improve at the same time the technology (from hot to cold crucible) and the associated matrix formulations, with constant emphasis on quality and volume reduction, ended up with the design and qualification of the cold crucible melter (CCM) technology.

Historical of numerical simulation development

Numerical simulation has been used since the beginning of the process development in the 80s. At first it was mainly used for designing the induction generator, the induction coil and was quickly applied for the designing of the cold crucible itself. Thanks to these first tools, mainly based on strong analytical studies, various induction furnace (hot melter or cold crucible) were designed and then used to melt metals and various oxides. Then the need of thermal simulation has pointed out and was performed in coupling with induction simulation. But thermal repartition is greatly influenced by the fluid flow. At this time (90s), computation capacity was too weak to perform full 3D computation and the development of numerical simulation has slow down. The figure 1 shows typical results obtained in simulation of natural convection driven flow for a glass melt in cold crucible. Simulations were only 2D axisymetrical with a coupling between the velocity and temperature field due to the variation of the viscosity with the temperature. In the early 2000s, commercial softwares specialized in induction and CFD become more efficient and thanks to the great increased of computation capacity, 3D simulations were successfully conducted. The coupling between 3D-CFD with 2D-induction has been done (Jacoutot 2006 [1]) followed with a 3D-CFD/3D-induction coupling (Sauvage 2009 [2]). Diphasic simulations of glass flow were added (Lima Da Silva 2014 [START_REF] Lima Da Silva | HES[END_REF]).

Modelisation description

In this section we will briefly described the simulation performed in our laboratory. All details will not be given here but are available in the cited reference [1,2,3]. Our tools are able to simulate induction heating from low to high frequency coupled to thermal simulation with the Joule power density in the metallic part of the process as well as in the molten glass itself. The flow simulation is also coupled and takes into account the mechanical stirring, gas injections and even Lorentz forces if necessary. The radiative heat transfer and the Marangoni convection are simulated as well. All these phenomena are computed on the full 3D real geometry of the furnace thanks to the use of supercomputing center as the TGCC (Très Grand Centre de Calcul) of the CEA. Software used was formerly academic software then commercial or open-source : Ansys-Fluent® for CFD , Flux® (Cedrat) for induction, OpenFoam®. The figure 2 shows an example of simulation of molten glass in a cold crucible with a cooled mechanical stirrer. quantification of the differences of potential and powers joules dissipated in the structures simulation used in testing new concept of furnaces of higher capacity.

The benefits are a higher power efficiency of the process, a longer life-time of optimized equipment, contribution to ensure the glass quality and of course conception of futur vitrification process.

Main issues of simulation of vitrification processes

Fluid flow simulation of glass has one major advantage in common fluid dynamics problems: it is only laminar flow. But in return, all the physical properties of the glass varies strongly with the temperature. Two problems arise: (i) a coupling between velocities and induced electrical current appears with temperature (ii) the law of the physical properties of the glass with the temperature are overriding.

The first point leads to develop numerical strategy to solve correctly the coupled physics especially between the fluid flow simulation and the electromagnetic simulation which are two physics solve by two different software (most of the time) and even with different numerical scheme (Finite Element and Finite Volume). This coupling has been achieved quite early in two dimensional simulation but the coupling between two 3D software has been reached only in 2009 [3].

The second point is the acquisition of physical properties of the glass. In our laboratory, characterizing technique of glass have been extensively developed from the early beginning of vitrification (in the 80s) but for other purpose than numerical simulation feeding. As a consequence, a lot of data of this time are inadequate to perform numerical simulation. For example, it was common to measure the viscosity of glass only at the temperature of elaboration whereas we need the full variation of the viscosity with the temperature. So effort have been done to adapt the existing characterization method of glass for numerical simulation purpose or even to develop new characterization technique to measure new properties (surface tension and thermal conductivity for example).

Validation work

As the simulation is often used in the development work, the need of validation of numerical results is very strong. Unfortunately, in-situ glass instrumentation at 1200°C is very limited. We perform a step by step validation of the simulation : (a) For the fluid flow simulation of glass we use the hydraulic similarity technique which consists in replacing the glass by oil with the same kinematic viscosity. As a result, the flow pattern of oil is similar to the glass one and instrumentation in oil is much easiest. This similarity technique allow us to perform quantitative comparison between CFD and experiment in the specific domain of highly viscous liquid which is not really widespread in the literature. (b) Electromagnetic simulation has been developed in the same time than cold crucible technology itself. The validation of such validation is made by magnetic field and induced current measurement and Joule losses power measurement in the cooled metallic part of the furnace (cf. figure 3). The main limitation of this validation is the precision of experimental measurement due to the high frequency and high voltage of induction generator. (c) For the full coupled simulation (CFD+Induction+Thermal) experimental comparison can be done only by direct comparison with glass experiment. To do so, we used specific devices to measure the thermal power coming from the glass to a cooled device. Precision of this kind of probing is quite low (~10%) and limit the validation of our simulation.

In order to optimize the effort of characterization enhancement, sensitivity analysis is performed with our simulation to target the physical properties having the most importance in the prediction of a specific result (wall heat flux, skullmelter thickness etc.). This kind of work may also be used as an uncertainty propagation method to determine the uncertainty of prediction of wall heat flux (for example) as a function of the uncertainties of physicals properties due to the acquisition method. This work helps us to assess the degree of confidence of our tool.

Modelisation development

As the simulation capability grows, the model is upgraded with new phenomena taken into account. The aim is to enhanced the prediction capability of simulation and widen the fields of its relevance. Different models are under development like a red-ox model to predict the fugacity of glass or an elaboration model for the glass from his raw materials inspired from Hrma et al. [4]. The figure 3 shows new developments in the simulation of bubbling in viscous liquid with the Volume Of Fluid (VOF) method compared to an experimental visualization of the same flow with oil.

Conclusions

The numerical simulation of vitrification furnaces has started in the same time of the first development of the process in the 80s. Thanks to the software and hardware development, the simulation takes a more and more important part in the current projects. Contributions of the simulation become valuable but simulation itself is not so easy. Physical properties measurement and data acquisition for validation are essential and specific development to feed the simulations are required. Also, even if commercial software becomes quite easy-to-use, the savoir-faire remains very important.

Fig. 1 .

 1 Fig. 1. Simulated velocity (m/s) and temperature field (K) of a natural convection driven flow of a glass load melt in cold crucible. 2D simulation with Flux-Expert® (1990).

Fig. 2 .

 2 Fig. 2. 3D simulation of temperature on the left (K) and right Joule power density (W/m3) in CCIMExamples of the contribution of numerical simulationNumerical simulation first objective is to help to understand physical phenomena within the molten glass. In fact, due to high temperature and high corrosive properties of the melt, accurate experimental data are very difficult to obtain. Thereafter, the simulation can help to optimize the design of structures or the high frequency power lines distribution. But also, the simulation can give an evaluation of thermal and chemical homogeneity and the impact of the glass properties on this homogeneity. Here is a list of applications: -optimization of internal water cooling of parts of the cold crucible melter, -optimization of glass mixing to assure a perfect thermal and chemical homogeneity, -studies to increase the throughput of glass production, -determination of residence time of material in the furnace, -induction optimization of inductor design, melting pot,

Fig. 3 .

 3 Fig. 3. Validation example of induction simulation of an empty cold crucible. The Joule losses in the parts of the crucible are compared for the same inductor current. On the left, magnetic field in a vertical cross plane of the cold crucible.

Fig. 4 .

 4 Fig. 4. Diphasic simulation of the interaction of two bubble columns in a viscous fluid representative of molten glass. Left, experimental view of oil, right CFD with Openfoam.

Acknowledgment

The authors thank AREVA for their financial support.