P Beckstein
email: p.beckstein@hzdr.de

V Galindo

G Gerbeth

A Schönecker

Free-surface dynamics in the Ribbon Growth on Substrate (RGS) process

Keywords: Ribbon Growth on Substrate, Semi-conductor processing, AC magnetic fields, Magnetic retention, OpenFOAM extensions, COMSOL Multiphysics, Moving mesh, Surface-tracking, Free-surface, Dome-shaping

The cost efficient, high throughput production of metal-and semiconductor alloys is the foundation of many advanced technologies. With the development of the Ribbon Growth on Substrate (RGS) technology, a new crystallization technique is available that allows the controlled, high crystallization rate production of silicon wafers and advanced metal-silicide alloys. Compared to other crystallization methods, such as melt spinning, the RGS process allows better crystallization control, high volume manufacturing and high material yield due to the substrate driven process. In order to optimize RGS further, insights from modelling the liquid metal in the casting frame under electromagnetic fields are very desirable. We performed numerical investigations in order to study the involved AC magnetic fields, which are an essential part of the RGS process to realize a magnetic retention effect. Our simulation results demonstrate the effect of the applied AC magnetic fields on the silicon melt flow. The main focus is thereby devoted to the simulation of the melt surface deformation based on a complex modelling approach. This time-dependent free-surface flow under the influence of magnetic forces is the key for optimizing the RGS process.

Introduction

The basic idea of the RGS process [1,2] for the production of silicon wafers is a continuous feeding of molten silicon into a bottomless casting frame, while a solidified silicon foil is extracted sidewise on a sub-cooled moving substrate underneath. Fig. 1 depicts this principle. The sketch represents one half of the core process viewing towards a central cutting section.

The AC magnetic fields which are used for inductive heating bear a two-fold meaning. That is, a single excitation coil provides both heat and a kind of magnetic valve. The latter actively prevents leakage in the slit regions and reduces oscillations at the extraction site of the silicon foil through electromagnetic forces. This Lorentz force field is acting to counter the gravitational forces on the melt, which only works in conjunction with surface tension.

We demonstrate the effect of the applied AC magnetic fields on the silicon melt in the RGS process, taking into account the timedependent coupling to the flow in the melt.

Fig. 1: Scheme of the RGS process with excitation coils [1].

Numerical Model

The casting region in the real RGS prototype machine [1] is a complex and detailed system. Therefore, a simplified set of modeling parameters and casting environment geometry was derived. It is expected that this numerical model allows the calculation of the main effects of the liquid silicon behavior and its results should still allow a comparison to real process data from the RGS process. For this simplified model the RGS wafer size is given by 156 𝑚𝑚 × 156 𝑚𝑚 × 0.4 𝑚𝑚. The mean melt level height inside the casting frame is assumed to be ℎ = 20 𝑚𝑚, the length of the melt region in process direction is 𝑙 = 70 𝑚𝑚 and the width 𝑤 = 150 𝑚𝑚. A typical set of process parameters, which evolved mainly from empirical analysis, is a RMS-current of 𝐼 𝑅𝑀𝑆 = 1000 A at a frequency of 𝑓 = 10 kHz to feed the excitation coils in combination with a substrate velocity of 𝑢 S = 0.1 m/s in process direction. The most important properties of materials shown in Fig. 1 The physical description of the magnetic fields in our model is based on the 𝑨-𝑉-formulation of the quasi-static Maxwell-Equations with implied MHD approximations. In addition, assuming only a very small magnetic Reynolds number and neglecting magnetization effects is valid for our purposes with sufficient quality. Introducing the velocity field 𝒖 and the physical time 𝑡 we may summarize these assumptions as follows:

𝑩 = 𝛁 × 𝑨 ; 𝑬 = -(∂ 𝑡 𝑨 + 𝛁𝑉); ‖𝒖 × 𝑩‖/‖𝑬‖ ≪ 1 . ( 1
)
Using the magnetic vector potential 𝑨 with applied Coulomb-Gauge (∇ • 𝑨 = 𝟎) and the electric scalar potential 𝑉 to describe the magnetic 𝑩 and electric field 𝑬 allows us to explicitly introduce a source current density 𝒋 E (𝐼 𝑅𝑀𝑆 ) term into the system which represents the effect of the excitation coil [START_REF] Binns | The Analytical and Numerical Solution of Electric and Magnetic Fields[END_REF][START_REF] Moreau | Magnetohydrodynamics[END_REF]. For constant electrical conductivity the scalar potential may be incorporated into a modified magnetic vector potential according to 𝑨 ̂= 𝑨 + ∫ 𝛁𝑉dt [START_REF] Morisue | [END_REF]:

𝛁 × 𝛁 × 𝑨 ̂+ 𝜎𝜇 0 ∂ 𝑡 𝑨 ̂= 𝜇 0 𝒋 E ; 𝒋 = 𝒋 I + 𝒋 E ; 𝒋 I = -𝜎 ∂ 𝑡 𝑨 ̂ . ( 2
)
The field 𝒋 𝐼 represents the purely induced current density. To apply suitable boundary conditions, the numerical domain is divided into several regions with constant 𝜎. A detailed mathematical description of all boundary/transition conditions can be found in [START_REF] Binns | The Analytical and Numerical Solution of Electric and Magnetic Fields[END_REF]. By introducing complex-valued sinusoidal fields, equation ( 2) can be transformed into the complex domain for an angular frequency of 𝜔 = 2𝜋𝑓. This approach leads to a quasi-stationary problem. For the momentum balance of the fluid, only the time-averaged Lorentz force 𝑭 𝑳 = 1/𝜌〈𝒋 × 𝑩〉 t is important. Supposing a negligible influence of natural convection, the fluid dynamics of the liquid melt (cf. Fig. 1) are governed by the principle of conservation of mass and momentum in form of the incompressible, isothermal Navier-Stokes-Equation [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] for a Newtonian fluid with additional terms for gravity and the time-averaged Lorentz-Force as described above:

𝜌[∂ 𝑡 𝒖 + (𝒖 • 𝛁)𝒖] = 𝛁 ⋅ 𝛕 + ρ𝒈 + 𝑭 𝐿 ; 𝛁 • 𝒖 = 0 ; (3) 𝛕 = 𝜂[𝛁𝒖 + (𝛁𝒖) 𝑇 ] -𝑝𝑰 . ( 4
)
Here 𝛕 represents the total stress tensor including the diagonal fluid pressure 𝑝. For our model only the fluid flow in the bulk region (cf. Fig. 1) is considered. The corresponding numerical domain is hence bounded by the top free-surface, stationary side walls and the liquid/solid interface (cf. Fig. 1). The latter is assumed to be planar and being translated with the casting velocity 𝒖 𝑆 of the moving substrate.

Surface tension is acting on the free-surface boundary. To keep our model simple, it is assumed to be constant and given by a constant coefficient 𝛾. It is thereby worth to mention that the viscosity of the external atmosphere, which is in contact with the liquid melt at the free-surface, is several orders of magnitude smaller than the viscosity of the melt itself. Thus, the fluid boundary condition at the free-surface can be modeled using a simplified Young-Laplace-Equation [START_REF] Scardovell | [END_REF] neglecting tangential drag forces from the atmosphere. Stationary walls are modelled with no-slip boundary conditions, whereas for the bottom wall of the fluid domain an inhomogeneous Dirichlet boundary condition is applied to simulate the substrate movement.

The free-surface movement and flow calculation is represented by means of a three-dimensional finite-volume surface tracking method based on OpenFOAM [8] with a dynamic moving mesh [START_REF] Scardovell | [END_REF][START_REF] Jasak | AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition[END_REF]. This utilized approach is known as Arbitrary Lagrangian-Eulerian formulation (ALE). In simplified terms, the essential idea of ALE is to allow the gridwhich is used for discretization -to move independently from the fluid flow, except for the domain boundaries. The free-surface is under constraint, such that the fluid velocity equals the mesh velocity 𝒖 𝑀 in normal direction there. For all other boundaries 𝒖 M is restricted according to the meaning of the boundary for the mesh (e.g. slip/no-slip). The independent mesh-movement away from the boundaries allows a free and preferably smooth mesh point distribution. In our case a Laplace-smoothing for 𝒖 m was utilized [START_REF] Jasak | Transactions of FAMENA[END_REF]. In the sense of ALE, equations ( 2) and (3) have to be considered as relative to the mesh motion 𝒖 m in a Lagrangian manner.

In order to realize the coupling of the magnetodynamic and hydrodynamic effects, a combination of a modified solver of OpenFOAM extensions [START_REF] Tukovic | [END_REF] and the finite-element software COMSOL Multiphysics [12] is used. The Lorentz force calculation with interpolated free-surface location-data is achieved with COMSOL Multiphysics which is iteratively coupled to OpenFOAM for the flow calculation with interpolated Lorentz force data. The backend-scripting for automation, interpolation and smoothing is written in Python and PyFoam. The surface-tracking approach on a moving computational domain allows us to "carry" the Lorentz force distribution 𝑭 𝐿 with the mesh for a defined and reasonable small simulation time, which allows us to greatly improve the overall performance.

Simulation results

Previous studies [START_REF] Beckstein | Proceedings of the 9th International PAMIR Conference on Fundamental and Applied MHD[END_REF][START_REF] Beckstein | Magnetohydrodynamics[END_REF] provided an overview of magnetohydrodynamic effects in the RGS process at its characteristic process parameters as given above. We have already successfully performed 3D-simulations with fixed melt geometry to numerically confirm the functioning melt retention based on magnetic fields. It has been illustrated that a retention effect is correlated with a strong forced fluid flow. A parameter study revealed the total system inductivity as a function of the melt level, which might be usable for a fill-level sensing of the silicon melt. In particular, we also demonstrated that the surface deformation is substantially important for a satisfactory model. The Lorentz force acting on a conducting liquid inside a coil is roughly pointing towards the center of this coil. Consequently, the shape of the deformed fluid domain for these and similar cases looks similar to a dome. This process is often referred to as "dome shaping" [START_REF] Moreau | Magnetohydrodynamics[END_REF]. For the RGS process parameters we revealed that the influence of the Lorentz force on the fluid flow is much more intense than the drag force from the moving substrate wall. Due to the challenging numerical conditions, several numerical tests were necessary in order to ensure stability of our modified OpenFOAM-solver. In the frame of this investigation we had to optimize our selection of finite-volume discretization schemes to account for possible high mesh skewness and a bounded solution, implement an automatic time step control, introduce a volume correction algorithm to obey mass conservation, improve parts of the surfacetracking procedure and fix bugs in OpenFOAM extensions. To increase solver performance and to allow unsteady RANS-simulations on coarse meshes, the usage of the class for turbulence and sub-grid-scale-models of OpenFOAM has been implemented in our solver, too. The 𝑘-𝜔-SST turbulence model [START_REF] Menter | [END_REF] was applied for all flow simulations.

As a proof of concept, Fig. 2 shows a time-series of the free-surface flow at the central section of the fluid domain with its evolving dome-shape. The magnetic field was applied at 𝑡 = 0 𝑠 for a constant surface level of ℎ = 20 𝑚𝑚 (Process direction is widthwise from left to right). In the left column one can comprehend how the Lorentz force density 𝑭 𝐿 /𝜌[𝑚/𝑠 2 ] is changing during time. The right hand side illustrates the corresponding velocity magnitude ‖𝒖‖ [𝑚/𝑠] and flow direction. For this calculation, the Lorentz force was updated every ∆𝑡 𝐿 = 0.01 𝑠. Since the magnetic vector potential as given by equation ( 2) is calculated for the whole simplified model (cf. Fig. 1) including casting frame, substrate and wafer, the result is not symmetric. This asymmetry is obvious for advanced time instances (𝑡 > 0.5 𝑠).

Even though one could spuriously hold the moving substrate wall at the bottom responsible for the large developing vortex flow in the central part, it is mainly induced by the asymmetric Lorentz force distribution.

Conclusion

The RGS process is a promising technology for the production of silicon wafers and advanced metal-silicide alloys. In order to optimize RGS, it is desired to deliver a better insight into the free-surface dynamics. We have developed a new solver concept to model the flow of conducting fluids under the influence of AC magnetic fields. On this basis we have successfully conducted unsteady simulations of the free-surface dynamics under typical conditions of the RGS process.

There are still a number of remaining problems to be addressed, especially for large 3D-simulations concerning numerical stability. This is subject of our current and future investigation. To validate our newly developed code, future work will be devoted to prepare different experimental test cases as a basis for quantitative comparison in the following step. There is also much room for improvement regarding the overall performance of the solution algorithm. Although we have already put much effort in reducing computational overhead due to the coupling of OpenFOAM and COMSOL Multiphysics to a minimum, it cannot be avoided completely with our current implementation. In this regard, we will develop a concept for assessing the feasibility of solving the partial differential equation ( 3) for the magnetic vector potential in OpenFOAM, too. This would allow us to solve all physics in just one framework.

Fig. 2 :

 2 Fig. 2: Evolution of the free-surface flow at the central section of the fluid domain with its dome-shape after applying the magnetic field for a constant surface level at 𝑡 = 0 𝑠 (left: 𝑭 𝐿 /𝜌[𝑚/𝑠 2 ]; right: ‖𝒖‖ [𝑚/𝑠]). The Lorentz force was updated every ∆𝑡 = 0.01 𝑠.

Table 1 :

 1 are listed in Table1. Material properties for different materials: Density 𝜌, kinematic viscosity 𝜂, electrical conductivity 𝜎, surface tension 𝛾 and skin depth 𝛿 = √1/(𝜋𝑓𝜇 0 𝜎) for a frequency of 𝑓 = 10 𝑘𝐻𝑧.

	Material	𝜌 [kg/m 3 ]	𝜂 [Pa ⋅ s]	𝜎 [S/m]	𝛾 [N/m] (Air)	𝛿 [mm]
	Liquid silicon	2580	0.86 x 10 -3	1.20 x 10 6	0.733	5.0
	Solid silicon	2330	-	8.30 x 10 4	-	17.0
	Graphite (casting frame)	1880	-	1.25 x 10 5	-	14.0
	Copper (inductor)	8960	-	6.00 x 10 7	-	0.7

Acknowledgements

Financial support for this research from the German Helmholtz Association in frame of the Alliance "Liquid Metal Technologies (LIMTECH)" is gratefully acknowledged. We would also like to thank our colleagues from "RGS Development B.V." for all their generous support and lively discussions.


Marie Lasserre

Stéphanie Bidon

Olivier Besson

MS Franc ¸ois Le Chevalier

Bayesian Sparse Fourier Representation of Off-Grid Targets

We consider the problem of estimating a finite sum of cisoids via the use of a sparsifying Fourier dictionary (problem that may be of use in many radar applications). Numerous signal sparse representation (SSR) techniques can be found in the literature regarding this problem. However, they are usually very sensitive to grid mismatch. In this paper, we present a new Bayesian model robust towards grid mismatch. Synthetic and experimental radar data are used to assess the ability of the proposed approach to robustify the SSR towards grid mismatch.

I. INTRODUCTION

In many radar applications, the received signal is conventionally described by a linear model where the signal of interest, namely the target signal, is represented by a sum of cisoids embedded in additive noise, i.e., y = N n=1 α n a n + n with [a n ] m = exp{j2πf n m} [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF] where y ∈ C M is the observation vector and M is the size of the observation space; α n , a n are respectively the complex amplitude and the steering vector with frequency f n of the nth target signal; n is the noise vector.

Several approaches can be used to estimate the target scene {(α n , f n )} related to the measurement y. These approaches are usually distinguished into two classes according to whether the method assumes, or not, a specific model about the noise covariance matrix. In the latter case the technique is said to be non-parametric (e.g., Fourier transform, Capon's method [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF], APES [START_REF] Li | An adaptive filtering approach to spectral estimation and SAR imaging[END_REF]) while in the former case it is called parametric (e.g., subspace method [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF], autoregressive model [START_REF] Hayes | Statistical digital signal processing and modeling[END_REF]). Over the last few decades, a new estimation paradigm called sparse signal reconstruction (SSR) has emerged. It has been applied to many signal processing applications and, in particular, to the target estimation problem (1), e.g., [START_REF] Anitori | Design and analysis of compressed sensing radar detectors[END_REF]. SSR aims at describing the signal as a linear combination of a few atoms from a (possibly pre-defined) dictionary. Towards this end and given

The work of M. Lasserre, S. Bidon and O. Besson is supported by the Délégation Générale de l'Armement under grant 2012.60.0012.00.470.75.01. the problem at hand [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF], a natural sparsifying dictionary is the Fourier basis so that (1) can be reformulated as

y = F x + n (2)
with

F ∈ C M × M the Fourier dictionary of size M × M where usually in SSR M ≥ M ; x ∈ C M
the sparse vector having ideally exactly N nonzero components.

The literature describing methods solving such SSR problem is nowadays prominent, e.g., ℓ 1 penalized least squares formulations [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] or Bayesian models [START_REF] Ji | Bayesian compressive sensing[END_REF]. However, most of these techniques are very sensitive to grid mismatch [START_REF] Herman | General deviants: An analysis of perturbations in compressed sensing[END_REF], [START_REF] Chi | Sensitivity to basis mismatch in compressed sensing[END_REF].

In the context of (2), grid mismatch occurs whenever the frequency of a target signal does not belong to the frequency grid associated with the sparsifying dictionary F , i.e., if f n / ∈ {0, 1/ M , . . . , 1 -1/ M }. Several attempts have been made at robustifying SSR technique towards grid mismatch. The most natural way to deal with grid mismatch is to refine the grid [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. However, most of the techniques choose to represent grid mismatch via a perturbation matrix E added to F [START_REF] Zhu | Sparsity-cognizant total leastsquares for perturbed compressive sampling[END_REF]. E usually stems from a first order Taylor expansion [START_REF] Yang | Off-grid direction of arrival estimation using sparse Bayesian inference[END_REF]- [START_REF] Tan | Sparse direction of arrival estimation using coprime arrays with off-grid targets[END_REF]. Another strategy is to consider that F is parameterized by a frequency grid that is estimated jointly with x [START_REF] Hu | Compressed sensing of complex sinusoids: An approach based on dictionary refinement[END_REF]. In this paper, we adopt the latter strategy and propose a hierarchical Bayesian model which takes into account the possible grid mismatch in a Fourier dictionary. The model assumes a white noise background and is an extension (aside from the dictionary) of a non-robust SSR method [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF]. The advantage of the proposed formulation over that of [START_REF] Yang | Off-grid direction of arrival estimation using sparse Bayesian inference[END_REF], [START_REF] Hu | A fast and accurate reconstruction algorithm for compressed sensing of complex sinusoids[END_REF], [START_REF] Hu | Compressed sensing of complex sinusoids: An approach based on dictionary refinement[END_REF] is that it enforces more sparsity via the use of a hierarchical Bernoulli-complex Gaussian prior on x. But mostly it is specifically designed to be used, in future work, in other estimation schemes recently developed by the authors.

The remaining of the paper is organized as follows. The proposed Bayesian model and its associated estimation scheme are described respectively in Section II and III. Numerical simulations are conducted in Section IV to assess the benefit of the robustification on both synthetic and experimental radar data.

II. BAYESIAN MODEL

The hierarchical Bayesian model considered in this paper is represented graphically in Fig. 1 and detailed herein. Compared to that of [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF], the technical novelty resides into the modeling of the grid mismatch.

A. Observation model 1) Modeling of grid mismatch:

We consider the observation model y = F x + n described in (2) with M fixed. In absence of grid mismatch, an appropriate choice for the Fourier dictionary F yields

[f m] m = 1/ √ M exp{j2πm( m/ M )}
where f m is the mth column of F . We propose to model the possible grid error by introducing a perturbation vector on the frequency axis denoted as

ε = [ε 0 . . . ε m . . . ε M -1 ]
T so that the Fourier dictionary is parameterized by ε as follows

F F (ε) = [f 0 (ε 0 ) . . . f m(ε m) . . . f M -1 (ε M -1 )]
where the mth column of F (ε) is now expressed as

[f m(ε m)] m = 1/ √ M exp{j2πm( m + ε m)/ M } .
To avoid overlapping between the frequency bins of F , the amplitude of the grid errors are assumed bounded such that for m = 0, . . . , M -1, ε m ∈ [-0.5, 0.5).

2) Likelihood: As stated in the introduction, an additive white noise background is considered. More precisely n is assumed to be centered Gaussian with power σ 2 , which is denoted as

n|σ 2 ∼ CN M 0, σ 2 I ( 3
)
where I is the identity matrix. The likelihood function is thus given by

f (y|ε, x, σ 2 ) = 1 π M σ 2M exp - y -F (ε)x 2 2 σ 2 . (4)
In ( 4), the vectors of interest are x and ε while σ 2 is a nuisance parameter. To estimate the target scene x, ε a Bayesian framework is chosen where each unknown parameter is modeled by a random variable with a given prior probability density function (pdf). In what follows, each prior density is designed to find a convenient balance between 1) mathematical tractability when performing the estimation 2) guarantee to preserve some physical sense to the hierarchical model.

B. Prior pdfs of the parameters

1) Target amplitude vector: Ideally in SSR the vector x introduced in (2) has exactly N nonzero elements whose value represents the post-integration amplitude of the target signal. As in [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF] a Bernoulli-complex Gaussian prior is chosen to actually enforce sparsity in x. More precisely, the elements x m [x] m of the amplitude vector are assumed independent and identically distributed (iid) according to the following mixed type pdf

f (x m|w, σ 2 x ) = (1-w)δ(|x m|)+w 1 πσ 2 x exp - |x m| 2 σ 2 x . (5)
Using the prior (5), denoted as x m|w, σ 2 x ∼ BerCN w, 0, σ 2 x , leads to considering that a target with power σ 2

x is present at the mth frequency with probability w.

2) Grid errors:

In this paper, we propose to define the prior pdf of the grid error ε m conditionally to the magnitude of x m. The idea behind this approach is that it may be unnecessary to estimate a grid error if no target signal is present at the corresponding frequency bin. More specifically, we assume that the ε m|x m are iid with pdf

f (ε m|x m = 0) = δ(ε m) (6a) f (ε m|x m = 0) = I [-0.5,0.5] (ε m) (6b)
where I A (.) is the indicator function of the set A. Note that the uniform distribution (6b) over [-.5, .5] is equivalent to a Betalike distribution of parameters (1,1) over [-.5, .5], denoted as Be [-.5,.5] (1, 1).

3) Noise power: A suitable prior for the white noise power σ 2 is an inverse-gamma distribution since it is conjugate to the likelihood (4). The prior pdf of σ 2 can therefore be expressed as

f (σ 2 |γ 0 , γ 1 ) ∝ e -γ1/σ 2 (σ 2 ) γ0+1 I [0,+∞) (σ 2 ) (7)
where γ 0 , γ 1 are respectively the shape and scale parameters. The distribution ( 7) is denoted as σ 2 |γ 0 , γ 1 ∼ IG (γ 0 , γ 1 ). Note that by tuning adequately the shape and scale parameters (γ 0 , γ 1 ), the prior can be made very informative or on the contrary flat. In radar applications, the thermal noise power is usually well known so that a moderately informative prior can be favored.

C. Prior pdfs of the hyperparameters

Since the probability w and the target signal power σ 2 x are both unknown, another level is added to the hierachical model.

1) Target signal power:

Similarly to the case of σ 2 , an inverse-gamma prior is chosen for the target signal power σ 2 x and is denoted as σ 2

x |β 0 , β 1 ∼ IG (β 0 , β 1 ). Nonetheless note that the shape and scale parameters β 0 , β 1 must be chosen to obtain this time a not so informative prior since the target signal amplitudes may vary significantly from one to another.

The MMSE estimator of x was intractable to derive analytically in [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF] when the grid mismatch was not taken into account. Therefore, the analytic calculation of the MMSE estimators (8a) and (8b) seem all the more trying. A Monte-Carlo Markov Chain (MCMC) is thus implemented [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. The MCMC algorithm simulates iteratively samples σ 2 (t) , ε (t) , x (t) , w (t) , σ 2

x (t) according to their conditional posterior distribution f (θ i |y, θ -i ) where θ = σ 2 , ε T , x T , w, σ 2

x T and θ -i is the vector θ whose ith element has been removed. After a burn-in time N bi , the samples are distributed according to their posterior distribution f (θ i |y). When a sufficient number of samples N r is collected, conventional Bayesian estimators can be built empirically

θiMMSE = N -1 r Nr t=1 θ i (t+N bi ) . (9)
The conditional posterior distributions are obtained from the joint posterior pdf of

σ 2 , ε, x, w, σ 2 x |y f (σ 2 , ε,x, w, σ 2 x |y) ∝ f (y|ε, x, σ 2 )f (ε|x)f (x|w, σ 2 x )f (w)f (σ 2 x )f (σ 2 ).
In particular, both vectors x and ε are sampled elementwise. The conditional posterior distributions of ε m and x m are derived from their conditional joint posterior distribution

f (ε m, x m|y, ε -m, x -m, σ 2 , w, σ 2 x ) ∝ exp -σ -2 |x m| 2 -x * mf m(ε m) H e m -x me H mf m(ε m) ×f (ε m|x m)f (x m|w, σ 2 x ) (10)
where e m = yi = m f i (ǫ i )x i .

A. Sampling of x

Following [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF], x is sampled element-wise. The mth element of x follows the distribution BerCN w m, µ m, η 2 m with

η 2 m = 1 σ 2 + 1 σ 2 x -1 µ m = η 2 m σ 2 f m(ǫ m) H e m w m = w η 2 m σ 2 x exp |µ m | 2 η 2 m 1 -w + w η 2 m σ 2 x exp |µ m | 2 η 2 m .

B. Sampling of ε

As for the vector x, the parameter ε is sampled elementwise. Using [START_REF] Chi | Sensitivity to basis mismatch in compressed sensing[END_REF], the conditional posterior distribution of ε m is calculated and rearranged to obtain

f (ε m|y, ε -m, x, σ 2 )∝ exp M -1 m=1 κ m cos 2π ε m M m -φ m ×f (ε m|x m) (11)
where

κ m = 2 σ 2 √ M ×|b m | and φ m = ∠b m . b has been defined as b = x * mu * m ⊙ e m, (12)
where u m = exp{j2πm( m/ M )}. When x m = 0, f (ε m|x m) = I [-0.5,0.5] (ε m). Thus, we recognize from (11) a dilated and truncated generalized von Mises distribution [START_REF] Gatto | The generalized von Mises distribution[END_REF]. Such a distribution can be troublesome to sample, so a Metropolis-Hastings (MH) algorithm is used [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. This algorithm is based on a proposal distribution that is easy to simulate from, and is chosen as close as possible from the target distribution. The conditional posterior distribution of ε m is represented in Fig. 2 for different values of mismatch and post-processing SNR (defined in ( 16)). This figure suggests that a flat proposal would be appropriate in the case of low power, and a Gaussian distribution in the case of high power. Thus, in our MH algorithm we switch from a flat proposal to a Gaussian proposal (and vice-versa) depending on the estimated target power. This adaptive scheme should be employed only during a burn-in period in order to preserve the convergence properties [START_REF] Robert | Monte Carlo Statistical Methods[END_REF], but in fact it does not damage the performance when used in the whole process.

C. Sampling of w, σ 2

x and σ 2 Like in [START_REF] Bidon | Bayesian sparse estimation of migrating targets for wideband radar[END_REF], the parameter σ 2 and hyperparameters w and σ 2

x are sampled according to their conditional posterior distribution

σ 2 |y, x, ε ∼ IG γ 0 + M, γ 1 + y -F (ε)x 2 2 (13) w|x ∼ Be (1 + n 1 , 1 + n 0 ) (14) σ 2 x |x ∼ IG β 0 + n 1 , β 1 + x 2 2 ( 15
)
where n 1 is the number of nonzero elements of x and

n 0 = M -n 1 . í f (ǫ m|y, ǫ -m, x, σ 2 ) ǫ m (a) í f (ǫ m|y, ǫ -m, x, σ 2 ) ǫ m (b) í f (ǫ m|y, ǫ -m, x, σ 2 ) ǫ m (c) í f (ǫ m|y, ǫ -m, x, σ 2 ) ǫ m (d)

IV. NUMERICAL SIMULATIONS

We now present various numerical examples illustrating the performance of the proposed SSR algorithm. Note that for each scenario, the constant hyperparameters (γ 0 , γ 1 ) and (β 0 , β 1 ) are calculated to obtain the desired values of mean and variance of the prior distributions of σ 2 and σ 2

x . Indeed for an inverse gamma distribution g ∼ IG (ν 0 , ν 1 ) the mean and variance are respectively

m g = ν 1 ν 0 -1 , ν 0 > 1 var g = ν 2 1 (ν 0 -1) 2 (ν 0 -2)
, ν 0 > 2 so the parameters are calculated as

ν 0 = m 2 g var g + 2 ν 1 = m g m 2 g var g + 1 .

A. Synthetic data

To begin with, synthetic data are generated according to the model described by ( 1) and (3).

1) Results after one realization: Three targets with possible grid mismatch are simulated in a white noise background with unit power. They all have a post-processing signal-to-noise ratio (SNR) of 25 dB, defined as

SN R n = E{|α n | 2 } σ 2 × M. (16)
In Fig. 3 are plotted the post-processing amplitude of the target vector xMMSE and the grid mismatch εMMSE . These estimates are compared with the true target scene x, ε and with the target scene estimated by the non-robust method (the proposed SSR technique when the grid mismatch is ignored, i.e., with ε = 0). From Fig. 3, we see that the proposed SSR algorithm estimates correctly both x and ε. On the contrary, if the grid error is ignored, the target energy tends to be spread over the adjacent frequency bins especially when the mismatch is important, i.e., |ε m| → .5. As as result it invalidates the sparsity of the reconstruction (this performance loss is well known).

2) Performance: A simple scenario is considered with a single target at zero velocity for different values of SNR and mismatch ǫ 0 . The performance of the proposed estimators ( 8) is assessed after 500 Monte-Carlo simulations through the calculation of the mean square error (MSE) of the reconstructed estimated target scene F (ε MMSE )x MMSE i.e.

MSE(F (ε

MMSE )x MMSE ) = (17) 1 N mc Nmc n=1 F (ε (n) MMSE )x (n) MMSE -F (ε)x 2 2 .
This metric seems more relevant than the MSE of xMMSE and εMMSE . Indeed, the sparse representation sometimes induces ambiguity about the position of the target: is it in the mth frequency bin with a mismatch of .5 (ε m = .5) or in the m + 1th frequency bin with a mismatch of -.5 (ε m+1 = -.5)?

With such an ambiguity, two problems might arise: the target might be shifted or split. Let us consider the case of a target in the mth frequency bin with a mismatch of .5 (ε m = .5).

First, the analysis may result in a shifted target: the target is estimated in the next range bin with a mismatch of -.5. The analysis may also result in a split target: a target is estimated in the mth frequency bin with a positive mismatch (ε m → .5), and another in the m + 1th frequency bin with a negative mismatch (ε m+1 → -.5). In both cases, the MSE of xMMSE and εMMSE will be high. On the other hand, the MSE of the reconstructed estimated target scene F (ε MMSE )x MMSE will be lower, thanks to the estimation of εMMSE . This metric may be more representative of the estimation quality.

We represent on Fig. 4 the MSE of F (ε MMSE )x MMSE as a function of grid mismatch for a SNR of 10 and 20 dB. We can see from Fig. 4 that for a SNR of 20 dB, the robust analysis outperforms the non-robust analysis as soon as ǫ 0 > .02 (ǫ 0 > .75 for a SNR of 10 dB) meaning that it is almost always profitable to estimate ε. We can also see that the benefits of the robust analysis compared to the non-robust analysis increase with the SNR. The clairvoyant case with respect to (wrt) grid mismatch is also represented on Fig. 4 and corresponds to the robust analysis when ε is known. It is a reference and it shows that when we know the true value of the grid mismatch, it is worth taking it into account for the estimation of the target scene. It also indicates how reliable the sampling of ε is in the proposed SSR method: it is highly efficient.

B. Experimental data

Let us now discuss the performance of the SSR algorithm in a practical case. To that end, data collected from the PARSAX radar [START_REF] Krasnov | Basics and first experiments demonstrating isolation improvements in the agile polarimetric FM-CW radar -PARSAX[END_REF] on November 2010 are considered. For this dataset the radar was pointing on a freeway during a heavy traffic time. The exact number of targets as well as their amplitude and location in the range-velocity map are unknown. Note that compared to the synthetic case, the target amplitudes are very high, so the hyperparameters (β 0 , β 1 ) are set such that (m σ 2

x , √ var σ 2 x ) = M dB + (60, 3.5) dB. In Fig. 5 are depicted the elements of xMMSE whose magnitude is greater than the lowest value of the colorbar. For the sake of comparison, the amplitude estimated by a conventional Capon's method is represented as a transparent background. (Note however that the Capon's output gives only an estimate of the target scene and should not be taken has the ground truth.) We see again the benefit of estimating the grid mismatch to avoid target spreading and/or miss detection. It is worth noticing that the target at range bin #4 and velocity ≈ -17 m/s is split by the robust analysis. This problem has been discussed before, and we can see that the non-robust estimation of this target is even worse. Besides, this anomaly is not inherent to our model and has been identified in [START_REF] Tan | Sparse direction of arrival estimation using coprime arrays with off-grid targets[END_REF] where a perturbation matrix is used.

V. CONCLUSION

In this paper we have described a new Bayesian algorithm for the sparse representation of targets in Fourier basis. Particularly, an error vector is introduced to model the possible mismatch between the target frequency and the nearest frequency point in the Fourier basis. The method, though computationally intensive, allows the sparsity of the scene to be preserved in case of grid mismatch. However, a problem persists in the limiting case when |ε m| → .5 but the proposed estimation method still remains more reliable than the nonrobust estimation. A number of parameters need to be set by the radar operator, and they can change significantly the performance of the analysis (in particular the parameters of the target signal power). This point will be further investigated in the near future. Furthermore, the proposed SSR technique will be integrated into more advanced hierarchical model.

Fig. 1 .

 1 Fig. 1. Directed acyclic graph associated with the hierarchical Bayesian model. The dashed lines indicate that (β 0 , β 1 ) and (γ 0 , γ 1 ) must be set by the operator.

2 )

 2 Level of occupancy: If no information is available to the radar operator about the sparsity level of the target scene, a convenient prior is a uniform pdf over the interval [0, 1], i.e., w ∼ U [0,1] . III. BAYESIAN ESTIMATION Herein we propose an estimation scheme of the target scene x, ε according to the Bayesian hierarchical model described in Section II. Particularly in what follows, we propose to study the minimum mean square error (MMSE) estimators xMMSE = xf (x|y)dx, (8a) εMMSE = εf (ε|y)dε.

Fig. 2 .

 2 Fig. 2. Shape of the conditional posterior distribution f (ε m|y, ε -m, x, σ 2 ). (a) SNR=5dB, ε m = .15. (b) SNR=5dB, ε m = .45. (c) SNR=15dB, ε m = .15. (d) SNR=15dB, ε m = .45.

Fig. 3 . 2 x, var σ 2 x

 322 Fig. 3. SSR of synthetic target scene: M = 16, N = 3, σ 2 = 1, M = M , (m σ 2 x , var σ 2 x ) = M dB + (0, 3.5) dB where M dB = 10 log 10 (M ), (m σ 2 , √ var σ 2 ) = (0, 2.4) dB. (a) Grid errors ε. (b) Target amplitude vectorx.

Fig. 4 . 2 x, var σ 2 x

 422 Fig. 4. Performance on synthetic target scene: M = 32, N = 1, σ 2 = 1, M = M , (m σ 2 x , var σ 2 x ) = M dB + (0, 3.5) dB, (m σ 2 , √ var σ 2 ) = (0, 2.4) dB. Comparison with the clairvoyant wrt ǫ and non-robust estimation. (a) SNR=10dB. (b) SNR=20dB

Fig. 5 . 2 x, var σ 2 x

 522 Fig. 5. SSR of experimental target scene (PARSAX data): range resolution of 1.5 m, M = 16, N unknown, σ 2 ≈ 1 (pre-normalization of the data), M = M , (m σ 2 x , var σ 2 x ) = M dB + (60, 3.5) dB, (m σ 2 , √ var σ 2 ) = (0, 2.4) dB. (a) Target amplitude vector xMMSE (diamond marker). (b) Target amplitude vector xMMSE when ε = 0 (diamond marker).

ACKNOWLEDGMENT

The authors would like to thank O. Krasnov at TU Delft for kindly providing the PARSAX experimental data.


D Krasnov
email: dmitry.krasnov@tu-ilmenau.de

Yu Kolesnikov
email: yuri.kolesnikov@tu-ilmenau.de

T Boeck

Numerical simulation of a round jet in the presence of an axial magnetic field

Keywords: direct numerical simulation, jet flow, flow instability

We report numerical simulations of a submerged round jet of conducting liquid issuing into a square duct in the presence of a uniform axial magnetic field. Starting from a quiescent initial condition, the spatial development of the jet is simulated until the transition or breakup zone is no longer evolving in space or the jet has reached the outlet while remaining intact. By increasing the magnetic field strength, the transition length of the jet grows and breakup can eventually be completely suppressed within the computational domain. For even stronger fields the laminar jet becomes unsteady again, which is presumably caused by a MHD pinching mechanism.

Introduction

Jet flows through nozzles occur in a number of metallurgical processes, e.g. in continuous casting machines. Electromagnetic brakes that are routinely applied in these machines also serve to control the jets from the submerged entry nozzle and the resulting flow pattern in the mold. How jets behave in the presence of magnetic fields is therefore not only of fundamental interest but also of practical importance. Generally, a magnetic field leads to flow anisotropy due to selective damping of flow structures. In a parallel flow, a transverse field gives rise to electromagnetic boundary layers on the walls, i.e the Lorentz force modifies the mean flow directly. A streamwise or axial field in a parallel flow only affects the fluctuations but not the mean flow. In particular, the formation of such electromagnetic boundary layers is avoided. In the present work we, therefore, focus on a streamwise (axial) field in a duct geometry. Specifically, we study the evolution of a jet that issues into a duct with initially quiescent liquid by means of direct numerical simulations. We are interested in the stabilizing effect (suppression of shear instability) of the magnetic field and in the behavior in a strong field. Theoretical and experimental investigations on submerged jets with an axial magnetic field have already been performed during the 1960s. However, they only give limited information. The theoretical investigations [1][2][START_REF] Plaschko | On the MHD-Instability of slowly spreading jets[END_REF] deal with the laminar flow and the behavior of small perturbations. Experiments [START_REF] Sajben | Measurement of the growth of a turbulent mercury jet in a coaxial magnetic field[END_REF][START_REF] Ievlev | Laminarization of a submerged jet of electrically conducting fluid by means of a longitudinal magnetic field[END_REF] are limited by the problematic flow measurement in the opaque and chemically aggressive liquid metals. Although simulations cannot reach the high Reynolds and Hartmann numbers of experiments, they provide full details on the flow structure and evolution. This information should be valuable for a better understanding of the complex dynamics of turbulent jet breakup and its modification by the magnetic damping.

Numerical procedure and parameters

We consider the flow of an incompressible, electrically conducting fluid (e.g. liquid metal) in a square duct subjected to a uniform magnetic field B, applied in the streamwise direction x. MHD effects are considered in the limit of small magnetic Reynolds number, i.e. the quasi-static approximation [START_REF] Davidson | An Introduction to Magnetohydrodynamics[END_REF] is applied. The governing non-dimensional equations are the Navier-Stokes equations with an additional Lorentz force term, which are given in details in [START_REF] Zikanov | [END_REF]. The non-dimensional parameters are the Reynolds number Re a =U q a/ and the magnetic interaction parameter N=Ha a 2 /Re a , where Ha a =Ba(/) 1/2 is the Hartmann number. Here U q is the mean flux velocity, a is the half-height of the duct,  is the electrical conductivity,  and  are, correspondingly, the density and kinematic viscosity. The governing equations are solved numerically by our in-house DNS solver, implemented for rectangular geometry and based on finite-difference method described in [8]. The spatial discretization of 2 nd order is performed on a rectangular structured grid, the time integration uses projection-correction procedure to satisfy incompressibility constraint. The computational grid is uniform in the streamwise x-direction and can be clustered in the wall-normal y,z-directions. The walls are considered perfectly insulating, i.e the electric potential /n=0 at the walls.

To address specific problems of jet flows, the non-periodic solver [START_REF] Zikanov | [END_REF] has been modified. First, the inlet conditions have been set to form a laminar round jet with the parabolic Poiseuille velocity profile that enters the computational domain. In our simulations the jet diameter d at the inlet is fixed and amounts to 1/8 of the full duct height. The jet is centered on the square inlet boundary. To characterize the jet we define additional Reynolds and Hartmann numbers based on the jet mean velocity U j and jet diameter d as characteristic scales. These values are commonly accepted in studies of jet flows and will be used further in the text. The second modification concerns boundary conditions at the exit.

The unsteady convective outflow boundary condition u/t+U q u/x=0 used in [2] was found to produce numerical artefacts (e.g. "zig-zag" or piling-up behavior) once the jet reached the exit section. The outflow condition has been replaced by a simplified form u/x=0, which has been tested and shown to produce numerically stable and smooth solutions of the velocity field.

A non-uniform grid with a modified Gauss-Lobatto clustering of points z=sin(/2) 2 sign() has been applied, where  and z are the uniform and wall-normal coordinates, correspondingly. Such type of grid-stretching provides tight clustering of the points not only at the walls, but also in the middle of the duct, where the jet is specified. In our grid-sensitivity study we have identified that this refinement of the core flow is really essential for the accurate resolution of the jet. The simulations have been conducted for the streamwise domain length L x =16 and the resolution of 2048x96x96 points. Further increase of the resolution did not change the results within considered Re and Ha. Finally, we remark that at each time-step the inlet jet-profile was perturbed with random noise with the amplitude 10 -3 in respect to the mean jet-velocity. This procedure has helped to regularize the process of jet instability. Otherwise, given the idealized numerical conditions, the jets could extend very far and remain laminar, which would require extremely long domain sizes and, correspondingly, computational time.

Results and discussion

In the absence of a magnetic field, jets are unstable already at fairly low Reynolds numbers (on the order of Re≈10) [START_REF] Danaila | Coherent structures in a round, spatially evolving, unforced, homogeneous jet at low Reynolds numbers[END_REF]. In our study the Reynolds number Re = U j d/ = 1800 exceeds such values significantly, and one observes a rapid transition to turbulence (Fig. 1a). At Ha=0, growing perturbations can already be observed at a distance of 2 diameters. Further downstream, these perturbations are destabilized and the jet expands rapidly. At the distance of 16d it has become fully turbulent and occupies the entire area of the duct. Further downstream, the flow reverts to the laminar state on account of the low global Reynolds number Re a ≈90, where turbulence in a duct flow cannot be sustained. In contrast to Fig. 1a, which shows a developed state of jet breakup after initial transients, Fig. 1(b-d) show the temporal evolution of the jet at Ha=25. The injection of the jet leads to a rapid development of a turbulent region (length of about 10-12d), which is advected downstream (Fig 1b). Behind it, the magnetic field suppresses instabilities. As time progresses, the turbulent region is increasingly damped and advected further (Figs. 1c,d) while it largely retains its length. The damping is apparent by the loss of small-scale features. In the final developed state of the jet, the turbulent region persists close to the exit. For Ha=50, the jet becomes fully laminar inside the computational domain in the developed state. The visible undulation in Fig. 2a is a transient feature. At Ha=125, the jet is stable and straight (Fig. 2b). Upon increasing the Hartmann number further this stable flow structure persists up to a certain value of Ha. For Ha=250 longitudinal waves appear on the jet as can be seen from the two snapshots of Fig. 2(c,d) that correspond to different times during jet development. The cause of these waves in a longitudinal magnetic field is related to the action of Lorentz forces. Due to the tendency of jet to expand a radial velocity component develops, which is accompanied by an annular electric current in the crosssection. Direction of the electric current depends on direction of radial velocity component. Interaction of the magnetic field with the current generates Lorentz force that counteracts the expansion or compression of jet cross-section. The evolution and expansion of the jet is further illustrated by Fig. 3, which shows snapshots of the distribution of axial velocity along the central axis of the jet, which coincides with the central axis of the computational domain. The velocity is displayed in units of the mean velocity U q. The curve for Ha=0 in Fig. 3a illustrates the developed state where the jet breaks down rapidly and the velocity reaches a constant value after a certain distance. For the higher values of Ha in Fig. 3a, the velocity distributions correspond to a transient state with approximately identical times from the beginning of the injection. One can clearly identify the turbulent region at the front at Ha=25. For Ha=50 and Ha=125 the jet has penetrated increasingly further and its front has become smoother. The front or leading edge visually resembles a running wave of the soliton type [START_REF] Davidson | An Introduction to Magnetohydrodynamics[END_REF] for flow regimes at Ha=50 and 125. The spreading of the jet, which is apparent from the decreasing velocity with the distance, is also reduced by increasing Ha from Ha=25 to Ha=125. Fig. 3b displays different snapshots of the axial velocity distribution for different times. One can see an early stage where the jet has only penetrated a certain distance (plot 1) and a later stage where the jet has almost reached the outlet (plot 2). It subsequently develops the oscillatory behavior (plots 3-5) described earlier with stronger variations in velocity for x<20.

Conclusions

Our results are in qualitative agreement with an experimental study of jet stabilization under longitudinal magnetic field that was conducted with the liquid metal alloy Galinstan at Re =3300 and different Hartman numbers around 100 [START_REF] Ievlev | Laminarization of a submerged jet of electrically conducting fluid by means of a longitudinal magnetic field[END_REF]. In this work it was found that laminarization of the jet also allowed its propagation on large distances. More detailed analyses of the expansion of the jet and its spatio-temporal dynamics for different values of Ha will be the subject of future work.

Fig. 1 :

 1 Fig. 1: Results of direct numerical simulations of jet flow at Re=1800 under streamwise magnetic field at Ha=0 (a) and three phases of evolution at Ha=25 (b,c,d). Shown is the streamwise velocity in a (x,y)-mid-plane at z=0.

Fig. 2 :

 2 Fig. 2: Results of direct numerical simulations of jet flow at Re=1800 under streamwise magnetic field at Ha=50 (a), 125 (b) and two-phases at Ha=250 (c,d). Streamwise velocity is visualized in a mid-plane, the radial "pinch" instability is visible at Ha=250.

Fig. 3 :

 3 Fig. 3: Variation of the centerline velocity along the streamwise direction for the cases of Ha=0, 25, 50 and 125 (a) and Ha=250 (b). The results at Ha=250 show five different phases of temporal evolution denoted by the indices running from 1 to 5 and corresponding to progress in time.

Acknowledgment D.K. and T.B. acknowledge financial support by the Helmholz alliance "Liquid Metal Technologies" (Limtech). Computer resources have been provided by the Forschungszentrum Jülich (Juropa supercomputer).


A Jakovics

S Pavlovs

E Baake

T Steinberg

V Sushkovs

A Tatulcenkovs

Numerical study of gas bubbles collective dynamics in the melt influenced with magnetic field

Keywords: bubbles reactor, Euler-Euler approach, Boltzmann approach, multiphase flow, magnetic field, LES

The paper presents numerical study of methane (CH 4 ) bubbles flow through molten tin (Sn) in electromagnetic (EM) field for cylindrical bubble reactor, which produces hydrogen (H 2 ) without carbonic acid (CO 2 ) emission. Eulerian-Eulerian approach, realized in commercial packages ANSYS CFX and ANSYS Fluent, as well as Boltzmann approach, realized as own code, are used for multiphase flowliquid metal (continuous phase) and gas bubbles (dispersed phase). Magnetic field is uniform. Transient distributions of methane flowing in liquid tin, influenced with magnetic field, are discussed.

Introduction

The production of hydrogen H 2 without generation of carbonic acid CO 2 may be performed with decomposition of methane CH 4 → C + 2H 2 . The reasonable rate of this endothermic reaction (Δh 0 = 74.85 kJ/mol) is reached at temperature values above 600ºC [1]. The main problem of such process is the decreasing of reactor life. The cause is deposition of carbon particles on the reactor walls and catalyst surface; as the consequence is the pressure drop and increase of thermal resistance. A possible way to increase efficiency of production and to control the process is the utilization of magnetic field, produced by direct (DC) or alternating (AC) current source, on electrically conductive liquid, when methane bubbles are flowing through it. The mentioned above process is realized with methane-tin cylindrical bubble reactor (CBR) [START_REF]Hydrogen from methane without CO 2 emissions[END_REF], the scheme of installation is shown in Fig.

1 [START_REF] Steinberg | [END_REF]. Main parameters of methane and tin [START_REF] Buhmirov | Reference materials for problems solution for "Heatand-mass-exchange" course[END_REF] are gathered in Tab. 1.

Previous consideration of single bubble dynamics [START_REF] Tucs | Proceedings of the 9 th PAMIR Int. Conference[END_REF], which covered both models implementation and experimental verification, are continued in current paper with numerical study of bubbles collective dynamics.

Multiphase Euler-Euler approach for methane bubbles flow through liquid tin

The systems of hydrodynamic equations for two-phase system of gas and liquid (denoted "g" and "l") consist of continuity, volume conservation and momentum equations according to Euler-Euler approach, which is realized in commercial packages -ANSYS CFX [START_REF]ANSYS CFX-solver theory guide[END_REF] and ANSYS Fluent [START_REF]ANSYS Fluent theory guide[END_REF]. The closure of equation system is performed with constrain on the pressure, namely two phases share the same pressure field. Static pressure in liquid depends on height z from the CBR bottom [START_REF] Lou | [END_REF]

  K kmol kJ  .
The flow is assumed isothermal at the operation temperature op T T 

. The liquid density is uniform constant.

In momentum equations the following interfacial forces, acting on one phase due to the presence of other phase [START_REF]ANSYS CFX-solver theory guide[END_REF][START_REF]ANSYS Fluent theory guide[END_REF], are taken into account: buoyancy force; drag force; lift force; virtual mass force; wall lubrication force.

The computations for tuning of numerical model show that only buoyancy and drag forces may be taken into account because changes of distribution of methane volume fraction due to other forces' contribution does not exceeded 5%.

The interphase turbulent dispersion force is not applicable for Large Eddy Simulation (LES) model of turbulence. This force may be used along with models of isotropic turbulence such as k-ω Shear Stress Transport (SST) model.

For computations of two-phase gas-liquid system commercial software packages are used, where the following discretization schemes for momentum, continuity and volume conservation equations are used:

• finite element method (FEM) in ANSYS CFX • volume of fluid (VOF) method in ANSYS Fluent, which is suitable for evaluation of liquid and gas volume fractions.

The test computations for the air-water flow in the rectangle bubble column have been performed. The characteristic ("snake-like") distributions of air volume fraction and velocity vectors field [START_REF] Steinberg | [END_REF][START_REF] Tucs | Proceedings of the 9 th PAMIR Int. Conference[END_REF] are obtained in 2D (mesh ~2 500 elements) and 3D (mesh ~32 500 elements) computations using steady-state or transient SST model of turbulence.

The numerical study of the methane-tin flow in CBR have been performed with LES models of turbulence using the structured 3D mesh (~400 000 elements) with the inflation in vertical direction from the inlet zone to the outlet zone.

To ensure the stability of computations for the transient flow regime the time step is chosen in the range 0.001-0.005 sec, which corresponds to the Courant number value less than the unity. When flow changes in time are extremely fast (magnetic field "switching on"), to gain the computations stability the time step is chosen < 0.0001 sec.

A lattice Boltzmann model for numerical simulation of dynamics of bubble flow

The lattice Boltzmann method (LBM) is a relatively recent computational approach based on kinetic theory for solving fluid mechanics and other physical problems [9,10]. In brief, the LBM consists of solving the lattice Boltzmann equation for the evolution of a distribution function

) , v , ( t x f  
as they move and collide on a lattice. LBM is originated from the Boltzmann equation:

) ( f f a f t f x            

Electromagnetic force

In the case of electrically conductive liquid, which flows with velocity l v  in external magnetic field with induction B  , the dynamics of two-phase system depends on EM force EM l f  , which is the results of interaction of electrical current, induced by liquid motion, and magnetic field,

B j f l EM l        B j l EM l l l         v    (1)
with volume fractions of liquid l  , electrical conductivity of liquid EM l  and electrical potential  .

Computations results of methane-tin flow in CBR

Main parameters of installation are gathered in Tab. 2. In the case of absence of magnetic field the methane jet has the instantaneous "snake-like" shape in the zone near inlet (Fig. 4 The methane jet is essential only near inlet with the maximum height of the jet penetration only 2.5% of the reactor height H CBR . The diameter of jet is comparable with the diameter D inlet of the inlet opening (Tab. 3). Above the inlet zone the distribution of methane is almost uniform. In the case when the constant vertical magnetic field is applied, the impact of Lorentz force (the maximum value EM f m ax ~3.3 •10 5 N/m 3 , the vectors are located in the horizontal plane and are perpendicular to the magnetic field induction) results in thinning-down of the methane jet (Fig. 4 (b)), which the averaged diameter may be estimated as 10% of the diameter D inlet of the inlet opening (Tab. 3). The penetration of extremely thin methane jet is very fast (archived flow time is approximately 3.5 sec) with timeaveraged speed aver ax meth m , v ~2.03 m/s, which is ~40.5 times greater than the methane speed in the inlet opening inlet meth v = 0.05 m/s (Tab. 3). The height of methane jet penetration is comparable with the reactor height H CBR (Tab. 3). In the case of the constant horizontal magnetic field, the vectors of Lorentz force (the maximum value EM f m ax ~5.0•10 4 N/m 3 ) are located in the vertical plane. The prevailing direction of Lorentz force is opposite to the tin flow direction, that is, the magnetic field decelerates the flow of tin and as sequence it decelerates the flow of methane. The diameter of methane jet is comparable with the inlet opening diameter D inlet (Fig. 4 (c), Tab. 3). This thick and straight methane jet penetrates to approximately 40% of the reactor height H CBR for archived flow time (10 sec) with the time-averaged speed aver ax meth m , v ~0.47 m/s, which is ~9.5 times greater than the methane speed in the inlet opening. • based on lattice Boltzmann method (LBM). √ Computational results for the methane and tin volume fractions and the velocity vectors in high temperature reactor without magnetic field are of a good qualitative agreement with the experimental and predicted results for the water-air two phase system in the rectangular bubble reactors. The further tuning of numerical approaches will be available on the base of the experiment data in gas and liquid metal system. √ The constant magnetic field is taken into account using simplified approach in order to estimate the impact of Lorentz force on the methane jet penetration into the liquid tin. The obtained results for distributions of the methane volume fraction show that the vertical magnetic field is narrowing the methane jet and consequently accelerates it. The horizontal magnetic field decelerates the jet and keeps the size of its diameter. The methane jets influenced with the magnetic field are straight, in contrast to the "snake-like" distributions of volume fractions and velocity vectors in the case without magnetic field.

Fig. 2 .

 2 simulation using the LBM model of the rising process of randomly generated bubbles in rectangular domain with 250×1250 lattice is carried out (Figs. 2, 3), the non-slip boundary conditions are imposed in all directions. The density ratio g l   between the fluid and gas is 1000 and viscosity ratio Rectangular column scheme. Fig. 3. Randomly generated bubbles in rectangular column x L = 0.5 m, y L = 2.5 m: t = 21.5 s; (b) t = 48.0 s; (c) t = 90.0 s.

  (a)). Time-averaged (flow time 40 s has been reached) distribution of methane volume fraction is almost axissymmetrical. Maximum time-averaged speed of methane is aver ax meth m , v ~0.33 m/s, which is ~6.5 times greater than speed in inlet opening inlet meth v = 0.05 m/s (Tab. 3). The instantaneous speed of methane may exceed inst are of the same range.

Fig. 4 :

 4 Fig. 4: Methane volume fraction instantaneous distributions in vertical cross-section of CBR inlet zone: (a) t = 40 s, B=0; (b) t = 13.5 s, B z = 0.5 T ("switched on" for t = 10-13.5 s; (c) t = 20 s, B x = 0.5 T ("switched on" for t =10-20 s).

Table 1 :

 1 Parameters of gas and liquid

	Parameter	Unit	Gas	Liquid
	Materials	-	methane (CH 4 )	tin (Sn)
	State	-	ideal gas	molten metal
	Operation temperature	°C	600	600
	Molar mass	kg/kmol	16.04	118.71
	Density	kg/m 3	0.224-0.353	6720
	Thermal expansivity	K -1	-	2.24•10 -5
	Dynamic viscosity	kg/(m•s)	1.42•10 -5	1.6•10 -3
	Gas-liquid surface tension coefficient	N/m -1	0.518

Fig. 1: Liquid metal bubble column for H 2 production without CO 2 emission.

Table 2 :

 2 Geometrical and operational parameters

	Geometrical parameters	Unit	Value	Operational parameters	Unit	Value
	Height of column	m	1.15	Gas speed at inlet	m/s	0.05
	Filing level of the melt	m	1.05	Gas volume flow	m 3 /s	2.0•10 -6
	Diameter of column	m	0.036	Gas mass flow	kg/s	7.5•10 -7
	Diameter of gas inlet	m	0.0072	Induction of magnetic field	T	0.5
	Diameter of gas bubbles	m	0.0001

Table 3 :

 3 Computed characteristic parameters of methane-tin flow

	magnetic induction [T]	flow time [s]	EM f m ax [N/m 3 ]	v	inst meth m ax , [m/s]	v	aver meth m ax , [m/s]	methane jet average diameter [m]	methane jet maximum height [m]
	B =0	40	0		1.16		0.33	~0.35D	inlet	~0.025 H	CBR
	B =0.5	3.5	3.3•10 5		3.26		2.03	~0.1D	inlet	~HCBR
	z
	x B =0.5	10	5.0•10 4		1.18		0.47	~Dinlet	~0.4H	CBR

Acknowledgment

The current research is supported by the Helmholtz Association of German Research Centers, within the scope of the Helmholtz Alliance -Liquid Metal Technologies (LIMTECH).


Ma Yong-Lin

Chen Yan-Jiang

Jiang Tao

Xing Shu-Qing

Design of a small Annular Linear Induction Electromagnetic Pump -Simulation and Experiment

Keywords: annular linear induction, electromagnetic pump, electromagnetic force, developed pressure, simulation

An annular linear induction electromagnetic pump (ALIP) is used for pumping liquid tin alloy with a flow rate and developed pressure of 17l /min and 0.5bar respectively. The electromagnetic force, flux density and inducted current density are calculated using a by finite element method (FEM) the various design variables. The performance of ALIP is tested with liquid tin for various melting point temperatures. Results show that pressure increases as electromagnetic forces. Comparing experimental and measured values, the relative error of simulation results was lower 15%.

Introduction

To pump molten metal annular linear induction electromagnetic pump (ALIP) has been widely investigated to transmit coolant, toxic heavy metals such as mercury, lead and so on. Early electromagnetic pumps (EMP) were applied to circulate sodium for nuclear reactor [START_REF] Baker | Handbook of Electromagnetic Pump Technology[END_REF] . Since the invention of the Cosworth Process able to put EMP into use for molten aluminum alloy in low pressures casting, the exploration and development of EMP was never stopped [START_REF] Campbell | 79 th Annual Conference[END_REF] . The planar DC EMP was applied to low pressure and casting of aluminum alloys. Advantages of EMP, are the simple construction, no mechanical motion which minimize oxidation of liquid metal and corrosion of the mechanical components and no maintenance that steadily carries out for a long period of time. Those reasons conduct ALIP to be the only option for a nuclear liquid metal reactor (LMR), rather than planar induction electromagnetic pump because of electrodes easy corrosion [START_REF] Saito | Electromagnetic pump type automatic molten-metal supply apparatus[END_REF][START_REF] Cheng | Planar DC electromagnetic pump for casting of Al-alloy[END_REF] . Simulation with 3D transient magnetic field, changing geometrical variables and load current is used for ALIP. In this paper, a small annular linear induction electromagnetic pump with nominal flow rate and static pressure of 17 l/min and 0.5bar respectively has been build with the help of ANSYS for for optimizing the design [5,[START_REF] Karasev | Magnetohydrodynamics[END_REF] .

Design analysis

ALIP could be divided into two parts, pump body serving as primary part of linear induction motor and the duct with inner core as the path of transmitting molten metal. The geometrical variables of ALIP, tooth width of outer core, pole pitch affecting the length of magnetic yoke iron core were considered according to 2p τ L  ; groove width and depth which affect coil number of turns were analyzed in the numerical simulation. The load variables such as current, voltage and frequency that change magnetic field intensity directly were determined by simulation. The developed pressure and flow rate were derived with electromagnetic force getting from finite element results [START_REF] Karasev | Magnetohydrodynamics[END_REF][START_REF] Kirillov | [END_REF] . The symbols definition in this paper is listed in Table 1.

Those design variables which affect developed pressure and mechanical efficiency are optimized using ANSYS software. The three-phase AC loading was employed and the edge method to settle 3D magnetic field of transient state was adopted during the analysis of ALIP. The differential form of Maxwell's equations, have been employed [START_REF] Kirillov | [END_REF]8] .

According to the requirement of geometrical variables, model dimensions of various parts are shown in Fig. 2. The transient analysis based on edge FEM (finite element modeling) is applied to reveal variation and spatial distribution of the magnetic flux, the induced currents and electromagnetic forces in the liquid tin. The simulation analysis is divided sixteen steps in a cycle. Each step could reveal the Z-axis of electromagnetic force (FMAGZ) by transient analysis. Here, the spatial distribution of magnetic flux density changes in liquid tin in the duct is calculated when inducting currents and frequency equal to 20A and 50Hz respectively [9] .

There are four frames corresponding to the vector distribution of the magnetic flux density at time T/4, T/2, 3T/4 and T in a calculation cycle.

Results of experimental analysis and comparison

Fig. 5 shows the ALIP fabricated combined with the design factors discussing above and the tin experimental facility for testing the designed ALIP. The experiment is conducted when the tin fluid temperature is 350°C. The measurement of the developed pressure was obtained by the raising liquid weight [9,10] . From Fig. 6, the FMAGZ and the raised liquid height( ΔH ) increased with the input voltage multiplied current (VI) in approximate linear relationship. The error between the experimental and numerical simulation values is equal to 11% in average. The difference is decreasing when the input VI voltage is increasing.

Conclusion

An annular linear induction electromagnetic pump has been fabricated based on the FEM calculations. The spatial distribution and evolution of magnetic flux, the induced current and electromagnetic force in liquid tin were investigated by a transient time analysis. The magnetic flux is mainly perpendicular to the duct and changed direction to opposite direction every interval of half cycle. The electromagnetic force is in the liquid flowing direction all time. In the operation temperature of 350℃, the desired performance of the ALIP was available and experimental was consistent with the simulation data.

Table 1 :Fig. 3 :

 13 Fig.2 : 3D simulation model of ALIP for the supply of liquid tin (a) Mesh model of mine components tetrahedron ; (b) Boundary conditions.

Fig. 4 (

 4 a) shows the magnetic flux density through the liquid tin which is generated by the outer core. The direction is changing every interval of T/2 and is always perpendicular to the z-axis . At any time, two opposite directions of magnetic flux density alternating distribution are existing along the axial .The spatial distributions of electromagnetic force with four different times are shown in Fig.4(b). Moreover, occupied most of the electromagnetic resultant force (FMAG) , FMAGZ, along with +Z direction, forced the liquid tin to move.So it is the effective component. FMAGX and FMAGY, the electromagnetic component forces along with + X and +Y respectively are not contributing to the driving force. The spatial distributions of FMAG components were reproduced each T/2 , i.e., the periodic transformation of FMAG was half that of AC power loaded to the coil.

Fig. 4 :

 4 Fig.4 : Spatial distribution of magnetic flux density and Lorentz force in liquid tin in the duct (20A, 50Hz) (a) The magnetic flux density ; (b) The Lorentz force

Fig. 5 :Fig. 6 :

 56 Fig.5 : Tin experimental facility for testing the designed ALIP

Acknowledgment

Financial support from the Innovation Foundation (Grant No.2009NC011) and the Foster Foundation (Grant No.PY-201004) are gratefully acknowledged. I am also deeply indebted to all the other researchers and teachers in Inner Mongolia University of Science and Technology for their direct and indirect help to me.