# Model-free Load Control for High Penetration of
Solar Photovoltaic Generation

Ouassim Bara1, Mohammad Olama4, Seddik Djouadi1, Teja Kuruganti4, Michel
Fliess25 and Cédric Join356  1Department of Electrical Engineering & Computer
Science, University of Tennessee, Knoxville, TN 37996, USA.
Email: {obara<EMAIL_ADDRESS>4Oak Ridge National Laboratory, Oak Ridge,
TN 37831, USA. Email:{olamahussemm<EMAIL_ADDRESS>2LIX (CNRS, UMR
7161), École polytechnique, 91128 Palaiseau, France. Email:
<EMAIL_ADDRESS>3CRAN (CNRS, UMR 7039), Université de
Lorraine, BP 239, 54506 Vandœuvre-lès-Nancy, France.
Email<EMAIL_ADDRESS>5AL.I.E.N., 7 rue Maurice Barrès, 54330
Vézelise, France. Email: {michel.fliess<EMAIL_ADDRESS>6Projet
Non-A, INRIA Lille – Nord-Europe, France

###### Abstract

This paper presents a new model-free control (MFC) mechanism that enables the
local distribution level circuit consumption of the photovoltaic (PV)
generation by local building loads, in particular, distributed heating,
ventilation and air conditioning (HVAC) units. The local consumption of PV
generation will help minimize the impact of PV generation on the distribution
grid, reduce the required battery storage capacity for PV penetration, and
increase solar PV generation penetration levels. The proposed MFC approach
with its corresponding intelligent controllers does not require any precise
model for buildings, where a reliable modeling is a demanding task. Even when
assuming the availability of a good model, the various building architectures
would compromise the performance objectives of any model-based control
strategy. The objective is to consume most of the PV generation locally while
maintaining occupants comfort and physical constraints of HVAC units. That is,
by enabling proper scheduling of responsive loads temporally and spatially to
minimize the difference between demand and PV production, it would be possible
to reduce voltage variations and two-way power flow. Computer simulations show
promising results where a significant proportion of the PV generation can be
consumed by building HVAC units with the help of intelligent control.

###### Index Terms:

HVAC, distributed energy resources, building load control, solar variability,
model-free control, intelligent proportional controllers.

##  I Introduction

A large introduction of solar photovoltaic (PV) and other renewable resources
may significantly alter the stability and balance of the power grid, through
voltage fluctuations and the ability to maintain the required network
frequency [1]. A viable solution to this problem, knowing the widespread use
of HVAC systems in residential and commercial buildings, is to use the idea of
load shaping. By introducing responsive loads (buildings), the control of HVAC
units will act as a counter to undesirable fluctuations of power generated by
the building sized PV arrays, while not sacrificing some level of temperature
comfort.

The motivation of not relying on model-based control approaches comes from the
daunting task of devising a reliable building model, either via physical laws
or via blackbox identification, and also, in the fact that it is difficult to
account for all the uncertainties and unknown disturbances, especially with
respect to strong weather disturbances and changes in occupancy pattern (see,
e.g., [2, 3, 4, 5, 6] and the references therein). A good model often requires
the installation of too many sensors necessary for the control algorithm
decision making part. Not to mention the financial cost that this will incur.
††This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy.

We are therefore proposing to use a new model-free control (MFC) approach with
its corresponding intelligent controllers [7], where the need of any precise
modeling disappears. This setting

  * •

is _data-driven_ , i.e., only the input and output data are used. Mathematical
physical laws and the associated differential equations are ignored.

  * •

is not limited to finite-dimensional linear systems. It has been successfully
applied in much more involved situations with nonlinearities and where
modeling via partial differential equations has been proposed.111Employing
PDEs is quite common in the air conditioning of buildings. See, e.g., [8] and
the references therein.

  * •

is easy to implement and quite robust with respect to internal and external
disturbances, even when compared to the most popular Proportional Integral
Derivative controller (see, e.g., [9]);

  * •

has already been successfully employed all over the world and gave rise to
some patents: see, e.g., the references in [7], [10] and the references
therein, and [11, 12, 13, 14, 15, 16, 17]. For obvious reasons let us
emphasize here the heating and the humidification of an agricultural
greenhouse [18], and the heating of a single building [19].

Although the idea of using HVAC unit systems to provide ancillary services to
the grid is not new (see [20], [21]), the main contribution of this paper is
on designing a MFC framework in order to consume most of the generated PV
power locally by building loads, while simultaneously ensuring a certain
temperature comfort inside the buildings. How to solve the actual problem
turns out to be intuitive, in the sense that one can consider the total energy
provided by the PV profile and share it amongst the available building load. A
straightforward solution was then to consider at each instant of time the
available energy coming from the profile as a constraint for each HVAC system.
Note that the difficulty in trying to implement such procedures is the
necessity for the control to follow two reference trajectories simultaneously
(see also [22]). One corresponding to the interior temperature of the
buildings (set in advance within a certain comfort zone), and the other,
corresponding to the generated PV profile. It is clearly a hard requirement to
satisfy if one wants to avoid optimization problem, where it is easier to
introduce those constraints. The price to pay for such optimization will
result in higher computational power, even if sub-optimality mathematical
techniques are introduced (see, e.g., [8, 23, 24, 25] and the references
therein). A real-time implementation would thus become problematic.

The remainder of the paper is organized as follows. Section II provides a
short introduction to model-free control control and to the corresponding
intelligent controllers.222See [7] for full details. A simplified physical
model is presented in Section III together with the formulation of our
approach and the simulation results. Section IV summarizes the paper and
presents the conclusions.

##  II Model-free control and intelligent controllers

Without any loss of generality we restrict ourselves to Single-Input Single-
Output systems. Instead of trying to write down a complex differential
equation, introduce the _ultra-local_ model

${\dot{y}=F+\alpha u}$ (1)

where

  * •

$u$ and $y$ are the input (control) and output variables,

  * •

the derivation order of $y$ is $1$, like in most concrete situations,

  * •

$\alpha\in\mathbb{R}$ is chosen by the practitioner such that $\alpha u$ and
$\dot{y}$ are of the same magnitude.

The following explanations on $F$ might be useful:

  * •

$F$ subsumes the knowledge of any model uncertainties and disturbances,

  * •

$F$ is estimated via the measures of $u$ and $y$.

###  II-A Intelligent controllers

The loop is closed by an _intelligent proportional controller_ , or _iP_ ,

$u=-\frac{\hat{F}-\dot{y}^{\ast}+K_{P}e}{\alpha}$ (2)

where

  * •

$y^{\star}$ is the reference trajectory,

  * •

$e=y-y^{\star}$ is the tracking error,

  * •

$K_{P}$ is the usual tuning gain.

Combining equations (1) and (2) yields:

$\dot{e}+K_{P}e=0$

where $F$ does not appear anymore. Local exponential stability is ensured if
$Kp>0$:

  * •

The gain $K_{P}$ is thus easily tuned.

  * •

Robustness with respect to different types of disturbances and model
uncertainties is achieved.

###  II-B Estimation of $F$

$F$ is estimated in real-time according to recent algebraic identification
techniques [26, 27, 28].

####  II-B1 First approach

The term $F$ in Equation (1) may be assumed to be “well” approximated by a
piecewise constant function $\hat{F}$. Rewrite then Equation (1) in the
operational domain (see, e.g., [29]):

$\displaystyle sY=\frac{\Phi}{s}+\alpha U+y(0)$ (3)

where $\Phi$ is a constant. We get rid of the initial condition $y(0)$ by
multiplying both sides on the left by $\frac{d}{ds}$:

$\displaystyle Y+s\frac{dY}{ds}=-\frac{\Phi}{s^{2}}+\alpha\frac{dU}{ds}$ (4)

Noise attenuation is achieved by multiplying both sides on the left by
$s^{-2}$. It yields in the time domain the real-time estimate, thanks to the
equivalence between $\frac{d}{ds}$ and the multiplication by $-t$,

${\small\hat{F}(t)=-\frac{6}{\tau^{3}}\int_{t-\tau}^{t}\left[(\tau-2\sigma)y(\sigma)+\alpha\sigma(\tau-\sigma)u(\sigma)\right]d\sigma}$
(5)

####  II-B2 Second approach

Close the loop with the iP (2):

$\hat{F}(t)=\frac{1}{\tau}\left[\int_{t-\tau}^{t}\left(\dot{y}^{\star}-\alpha
u-K_{P}e\right)d\sigma\right]$ (6)

###### Remark II.1

Note the following facts:

  * •

Integrals (5) and (6) are low pass filters.

  * •

$\tau>0$ may be chosen quite small.

  * •

Integrals (5) and (6) may of course be replaced in practice by classic digital
filters.

##  III Problem formulation and computer simulations

###  III-A A simple mathematical model

Computer simulations require obviously some mathematical modeling. We have
selected a simple set of linear differential equations, which are time-
invariant, i.e., with constant coefficients. It was derived in [30] and
successfully used in [23, 31, 32]. The dynamic of the interior temperature,
interior wall surface temperature, and exterior wall core temperature are
given by

$\displaystyle\dot{T}_{1}$
$\displaystyle=\frac{1}{C_{1}}\big{[}(K_{1}+K_{2})(T_{2}-T_{1})+K_{5}(T_{3}-T_{1})+u_{c}$
$\displaystyle+\delta_{2}+\delta_{3}\big{]}$ $\displaystyle\dot{T}_{2}$
$\displaystyle=\frac{1}{C_{2}}\big{[}(K_{1}+K_{2})(T_{1}-T_{2})+\delta_{2}\big{]}$
$\displaystyle\dot{T}_{3}$
$\displaystyle=\frac{1}{C_{3}}\big{[}K_{5}(T_{1}-T_{3})+K_{4}(\delta_{1}-T_{3})\big{]}$
(7) $\displaystyle\dot{x}=A_{c}x+B_{c}u+C_{c}w$ (8)

where the state vector $x$, the control $u$ and the disturbance $w$ are given
respectively as

$\displaystyle x=\begin{bmatrix}T_{1}\\\ T_{2}\\\ T_{3}\end{bmatrix},\quad
u=u_{c},\quad w=\begin{bmatrix}\delta_{1}\\\ \delta_{2}\\\
\delta_{3}\end{bmatrix}$

$T_{1}$: room air temperature [${}^{\circ}C$] $22\leq T_{1}\leq 24$,
$T_{2}$: interior-wall surface temperature [${}^{\circ}C$],
$T_{3}$: exterior-wall core temperature [${}^{\circ}C$],
$u_{c}$: cooling power ($\leq 0$) [$kW$],
$\delta_{1}$: outside air temperature [${}^{\circ}C$],
$\delta_{2}$: solar radiation [$kW/m^{2}$],
$\delta_{3}$: internal heat sources [$kW$],
$C_{1}=9.356\times 10^{5}$ kJ/C, $C_{2}=2.970\times 10^{6}$ kJ/C,
$C_{3}=6.695\times 10^{5}$ kJ/C, $K_{1}=16.48$ kW/C, $K_{2}=108.5$ kW/C,
$K_{3}=5$, $K_{4}=30.5$ kW/C and $K_{5}=23.04$ kW/C.

The above equations are only valid for a summer day. The main control variable
is therefore cooling power $u_{c}$. One of our control objectives is to
regulate the interior temperature of several buildings around 23
${}^{\circ}C$. This is the reference trajectory.

###### Remark III.1

A poor knowledge of the coefficients in Equation (8) is unavoidable in
practice. Lack of space prevents us of confirming via simulations the
robustness of our MFC setting in such a situation.

###  III-B Temperature regulation and PV following

Our aim is to regulate the interior temperature of multiple buildings while at
the same time assuring the whole consumed energy tracks as closely as possible
the PV energy profile. It yields the following control objectives:

  * •

Maintain the interior temperature under an acceptable comfort zone, i.e.,
between $22^{\circ}C$ and $24^{\circ}C$.

  * •

Follow the generated PV profile as closely as possible.

###### Remark III.2

The PV profile is the only energy that the HVAC units are supposed to get. The
buildings are also connected to the grid. It means that the total energy
consumed is allowed to fluctuate around the generated PV signal.

Instead of considering the PV profile as a constraint, as it has been
mentioned in the introduction, i.e.,

$\displaystyle\sum_{i=1}^{N_{b}}u_{i}(t)\leq PV(t)$ (9)

where $N_{b}$ is the number of buildings, a better alternative is to assume a
band around the PV profile, such that the available energy is given according
to

$\displaystyle E_{d}(t):\begin{cases}E_{d}(t)=0\quad\text{if}\quad PV(t)=0\\\
PV(t)-\epsilon\leq E_{d}(t)\leq PV(t)+\epsilon\quad\text{if}\quad
PV(t)>0\end{cases}$

where

  * •

$\epsilon>0$ is a constant, which is chosen according to the grid capabilities
and the comfort resulting from the interior temperature;

  * •

it is assumed that, when $E_{d}(t)=0$, all the buildings are connected to the
grid and no PV tracking is performed.

  * •

The constraint $0\leq u_{c}\leq 3$ on the cooling power must be satisfied when
there is no PV energy.

Therefore the total energy consumption is

$\displaystyle PV(t)-\epsilon\leq\sum_{i=1}^{N_{b}}u_{i}(t)\leq
PV(t)+\epsilon$ (10)

Assume for simplicity’s sake that all the buildings are identical. Then
Inequality (10) yields

$\displaystyle\frac{PV(t)}{N_{b}}-\frac{\epsilon}{N_{b}}\leq
u_{i}(t)\leq\frac{PV(t)}{N_{b}}+\frac{\epsilon}{N_{b}},\quad i=1,2,..N_{b}$
(11)

###  III-C Computer simulations

A $10$ minute sampling is used.333The time lapse depends heavily, of course,
on the nature of the plant we want to regulate. Consider first 13 buildings.
The simulations are presented first when only the interior temperatures
$T{{}_{1}}_{i}$, $i=1,\dots N_{b}$, is regulated. Afterwards a second
objective is added, i.e., the tracking of the PV generated signal, i.e.,
Inequality (10) must be satisfied.

####  III-C1 Control design

Since the interior temperatures $T{{}_{1}}_{i}$ are available for
measurements, Equation (1) yields the ultra-local model

$\displaystyle\dot{T}{{}_{1}}_{i}=F_{i}+\alpha_{i}u{{}_{c}}_{i}$ (12)

From Equation (2) we deduce the following intelligent proportional controller

$\displaystyle
u{{}_{c}}_{i}=-\frac{\hat{F}_{i}-\dot{T_{1}}^{*}_{i}+K{{}_{P}}_{i}e{{}_{p}}_{i}}{\alpha_{i}}$
(13)

where

  * •

$\alpha_{i}=5$, $K{{}_{P}}_{i}=2$,

  * •

${e_{p}}_{i}=T{{}_{1}}_{i}-{T_{1}}^{*}_{i}$ is the tracking error,

  * •

the reference trajectory ${T_{1}}^{*}_{i}$ should be as close as possible to
the reference temperature $23^{\circ}C$,

  * •

$\hat{F}_{i}$ is an estimate of $F$, which is obtained according to Section
II-B,

  * •

for a “good” PV tracking each control $u{{}_{c}}_{i}$ must satisfy Inequality
(11).

####  III-C2 Model-free control with no PV tracking

Figure 1 displays the evolution of the interior temperature during three days.
The temperature is perfectly regulated. MFC is very efficient in rejecting
outside disturbances: one can see their effect in the little bumps appearing
essentially from the $8^{th}$ to the $20^{th}$ hour of the day. Note also that
the initial temperatures in each building are different. The control inputs in
Figure 2 satisfy the required constraints. The initial large values correspond
to the buildings with the highest initial temperature. It is important to
notice in Figure 3 that the total energy used by the HVAC units does not take
into account the variations of the PV profile, hence the large deviation of
the total energy signal represented in blue dashed lines from the actual PV
given in red solid line.

Figure 1: Interior temperature variations in ${}^{\circ}C$ for $N_{b}=13$
Figure 2: Time evolution of the $13$ control inputs  Figure 3: Total energy
$\sum_{i=1}^{N_{b}}{u_{c}}_{i}(t)$ and the PV profile

####  III-C3 Simulation with PV profile

As soon as the photovoltaic source starts collecting energy, Inequality (11)
is taken into account. Figure 4 displays the time evolution of the interior
temperature for 13 buildings while the PV tracking is taken into account. Note
that most of the temperatures are within their prescribed comfort zone. The
slight violation of the lower constraints for one or two buildings is due to
the tight restriction on the PV tracking as shown by Figure 5: we select
$\epsilon=1$, i.e., a rather narrow tolerance margin. Temperatures where
slight overshoots appear are those that are starting from the lowest initial
temperatures. By imposing more energy consumption for the buildings than
necessary for regulating its indoor temperature, one would expect the HVAC
units to run during more time or with higher intensity. It would yield large
variations. This behavior is rather normal. It can be significantly improved
by selecting an appropriate number of buildings. It is important to understand
the difficulty of the actual compromise between the requirement of imposing a
certain level of comfort inside the buildings, and that of closely following
the generated PV signal. Those two requirements are very often contradictory.
Our viewpoint helps in satisfying both control objectives without the need of
any complex optimization procedure, especially if one has to deal with a
nonlinear model. As an initial study we decided to focus on a limited number
of HVAC units to assess our results. However scaling the problem to larger
building units is straightforward. It does not necessitate any change in the
problem formulation. More importantly it will not require any additional
significant computational power.

As already mentioned, it is possible to improve the results in Figure 4. Add
for that purpose one building, i.e., $N_{b}=14$. Figure 7 shows a clear-cut
improvement with respect to the interior temperature. Figures 8 and 9 do not
display on the other hand any significant change with respect to the case of
13 buildings. One needs to be cautious when dealing with an unfeasible
scenario where the maximum energy consumed by all buildings is very small when
compared to the energy coming from the PV source, or when the PV energy is not
sufficient to account for all the HVAC units. Thus, a better temperature
regulation emphasizes that we are, in fact, in the scenario where more energy
was available from the PV for the 13 buildings. Adding one more helped in
distributing the extra energy to the new building, putting less weight on the
remaining ones to follow the PV profile.

The above results are quite good. Remember that enabling the buildings to
absorb the PV sources will help the latter to minimize their negative effect
on the distribution grid and consequently to provide responsive loads that
would ultimately help reducing voltage variations and two-way power flow.

Figure 4: Interior temperature variations in ${}^{\circ}C$ for $N_{b}=13$

##  IV Conclusion

In order to offer a better integration of the photovoltaic energy sources into
the grid, a model-free control approach has been proposed. Computer
simulations are showing that the responsive loads (HVAC units) may be
regulated in such a way that the tracking error between demand and PV
generation is minimized. This fact will help reducing voltage variations and
two-way power flow. The two main control objectives mentioned in the
introduction are thus satisfied:

  * •

the interior temperatures inside each unit remains in the band $\pm
1^{\circ}C$ around the reference temperature,

  * •

the tracking error of the PV profile is less than $1$kW.

A forthcoming publication will extend this work to Boolean controls, i.e.,
on/off controls. This step should lead to concrete implementations. Let us
emphasize finally that the hardware implementation of our control strategy is
cheap and easy [33].

Figure 5: Total energy $\sum_{i=1}^{N_{b}}{u_{c}}_{i}(t)$ and the PV profile
($N_{b}=13$)

## ACKNOWLEDGMENT

This material is based upon work supported by the U.S. Department of Energy,
Office of Energy Efficiency & Renewable Energy, the SunShot National
Laboratory Multiyear Partnership (SuNLaMP) program under Contract No. DE-
AC05-00OR22725.

Figure 6: Time evolution of the $13$ control inputs

## References

  * [1] Y.-S. Kim, I.-Y. Chung, S.-I. Moon.  An analysis of variable-speed wind turbine power-control methods with fluctuating wind speed.  _Energies_ , 6, 3323-3338, 2013.
  * [2] D. Jovcic, N. Pahalawaththa, M. Zavahir. Analytical modelling of HVDC-HVAC systems. _IEEE Trans. Power Deliv._ , 14, 506-511, 1999.
  * [3] R. Montgomery, R. McDowall. _Fundamentals of HVAC Control Systems_. Elsevier, 2009.
  * [4] R. Teodorescu, M. Liserre, P. Rodriguez.  _Grid Converters for Photovoltaic and Wind Power Systems_.  Wiley, 2011.
  * [5] M.E. Raoufat, A. Khayatian. A. Mojallal.  Performance recovery of voltage source converters with application to grid-connected fuel cell DGs. _IEEE Trans. Smart Grid_ , 2016.
DOI: 10.1109/TSG.2016.2580945

  * [6] D. Zhou, Q. Hu, C. J. Tomlin. Model comparison of a data-driven and a physical model for simulating HVAC systems. _Amer. Contr. Conf._ , Seattle, 2017.
  * [7] M. Fliess, C. Join. Model-free control. _Int. J. Contr._ , 86, 2228-2252, 2013.
  * [8] S.D. Howell, P. V. Johnson, P. W. Duck. A rapid PDE-based optimization methodology for temperature control and other mixed stochastic and deterministic systems. _Energ. Buildings_ , 43, 1523-1530, 2011.
  * [9] K. J. Åström, T. Hägglund. _Advanced PID Control_ , Instrument Soc. Amer., 2006.
  * [10] H. Abouaïssa, M. Fliess, C. Join. On ramp metering: Towards a better understanding of ALINEA via model-free control. _Int. J. Contr._ , 90, 1018-1026, 2017.
  * [11] O. Bara, M. Fliess, C. Join, J. Day, S. M. Djouadi. Model-free immune therapy: A control approach to acute inflammation.  _Europ. Contr. Conf._ , Aalborg, 2016.
  * [12] H. Wang, X. Ye, Y. Tian, G. Zheng, N. Christov. Model-free based terminal SMC of quadrotor attitude and position. _IEEE Trans. Aerosp. Electron. Syst._ , 52, 2519-2528, 2016.
  * [13] J. Pardelhas, M. Silva, L. Mendonça, L. Baptista. Speed control of an experimental pneumatic engine. P. Garrido, F. Soares, A.P. Moreira (Eds): _Controlo 2016_ , Lect. Notes Electr. Engin. 402, Springer, 2016.
  * [14] S. Maalej, A. Kruszewski, L. Belkoura. Robust control for continuous LPV system with restricted-model-based control. _Circ. Syst. Sign. Process._ , 36, 2499-2520, 2017.
  * [15] T. MohammadRidha, M. Ait-Ahmed, L. Chailloux, M. Krempf, I. Guilhem, J. Y. Poirier, C. H. Moog. Model free iPID control for glycemia regulation of type-1 diabetes. _IEEE Trans. Biomed. Engin._ , 2017.
DOI: 10.1109/TBME.2017.2698036

  * [16] L. Menhour, B. d’Andréa-Novel, M. Fliess, D. Gruyer, H. Mounier. An efficient model-free setting for longitudinal and lateral vehicle control. Validation through the interconnected pro-SiVIC/RTMaps prototyping platform. _IEEE Trans. Intel. Transport. Syst._ , 2017.
DOI: 10.1109/TITS.2017.2699283

  * [17] J. -H. Jeong, D.-H. Lee, Mi. Kim, W.-H. Park, G.-S. Byun. The study of the electromagnetic robot with a four-wheel drive and applied I-PID system. _J. Electr. Eng. Technol._ , 12, 2017.
  * [18] F. Lafont, J. -F. Balmat, N. Pessel, M. Fliess.  A model-free control strategy for an experimental greenhouse with an application to fault accommodation.  _Comput. Electron. Agricult._ , 110, 139-149, 2015.
  * [19] H. Abouaïssa, O. Alhaj Hasan, C. Join, M. Fliess, D. Defer. Energy saving for building heating via a simple and efficient model-free control design: First steps with computer simulations. _Submitted_.
  * [20] H. Hao, Y. Lin, A.-S. Kowli, P. Barooah, S. Meyn.  Ancillary service to the grid through control of fans in commercial building HVAC systems.  _IEEE Trans. Smart Grid_ , 5, 2066-2074, 2014
  * [21] M. Balandat, F. Oldewurtel, M. Chen, C. Tomlin.  Contract design for frequency regulation by aggregations of commercial buildings.  _52nd Allerton Conf. Commun. Contr. Comput._ , Monticello, 2014.
  * [22] C. Join, J. Bernier, S. Mottelet, M. Fliess, S. Rechdaoui-Guérin, S. Azimi, V. Rocher. A simple and efficient feedback control strategy for wastewater denitrification. _20th World IFAC Congr._ , Toulouse, 2017.
  * [23] C. D. Nelson. _Optimal Control of Energy Efficient Buildings_. Master Thesis, Univ. Tennessee, Knoxville, 2014.
  * [24] L. Zheng, L. Cai. A distributed demand response control strategy using Lyapunov optimization. _IEEE Trans. Smart Grid_ , 5, 2075-2083, 2014.
  * [25] L. A. Hurtado, P. H. Nguyen, W.L. Kling. Smart grid and smart building inter-operation using agent-based particle swarm optimization. _Sustain. Energ. Grids Netw._ , 2, 32-40, 2015.
  * [26] M. Fliess, H. Sira-Ramírez.  An algebraic framework for linear identification.  _ESAIM Contr. Optimiz. Calc. Variat._ , 9, 151-168, 2003.
  * [27] M. Fliess, H. Sira-Ramírez.  Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques. H. Garnier & L. Wang (Eds): _Identification of Continuous-time Models from Sampled Data_ , Springer, pp. 362-391, 2008.
  * [28] H. Sira-Ramírez, C. García-Rodríguez, J. Cortès-Romero, A. Luviano-Juárez. _Algebraic Identification and Estimation Methods in Feedback Control Systems_. Wiley, 2014.
  * [29] K. Yosida. _Operational Calculus_ (translated from the Japanese). Springer, 1984.
  * [30] M. Gwerder, J. Tödtli.  Predictive control for integrated room automation.  _8th REHVA World Congr. Bldg Techno. – CLIMA 2005_ , Lausanne, 2005.
  * [31] F. Oldewurtel, A. Parisio, C. N. Jones, M. Morari, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, K. Wirth.  Energy efficient building climate control using stochastic model predictive control and weather predictions.  _Amer. Contr. Conf._ , Baltimore, 2010.
  * [32] R. Gondhalekar, F. Oldewurtel, C. N. Jones. Least-restrictive robust periodic model predictive control applied to room temperature regulation. _Automatica_ , 49, 2760-2766, 2013.
  * [33] C. Join, F. Chaxel, M. Fliess. “Intelligent” controllers on cheap and small programmable devices. _2nd Int. Conf. Contr. Fault-Tolerant Syst._ , Nice, 2013.

Figure 7: Interior temperature variations in ${}^{\circ}C$ for $N_{b}=14$
Figure 8: Total energy $\sum_{i=1}^{N_{b}}{u_{c}}_{i}(t)$ and the PV profile
($N_{b}=14$) Figure 9: Time evolution of the $14$ control inputs


# A short proof of Ledoit-Péché’s RIE formula for covariance matrices

Florent Benaych-Georges  CFM, 23 rue de l’Université, 75007 Paris, France
<EMAIL_ADDRESS>
###### Abstract.

This is a short proof of Ledoit-Péché’s RIE formula for covariance matrices.
The proof is based on the Stein formula, which gives a very simple way to
derive the result. One of the advantages of this approach is that it shows
that the only really needed hypothesis, for the machinery to work, is that the
mean of the eigenvalues of the true covariance matrix and the largest of them
have the same order.

Notation: For $M$ a real matrix, $M^{\prime}$ denotes the transpose of $M$,
$\|M\|_{\text{F}}$ denotes the Frobenius norm of $M$ defined at (2) and
$\|M\|$ denotes the operator norm of $M$ (with respect to the canonical
Euclidian norms). For $Z$ a random variable, $\operatorname{\mathbb{E}}Z$
denotes the expectation of $Z$.

##  1\. Model

Let $n$ and let $X\in\mathbb{R}^{n}$ be a centered Gaussian (column) random
vector with covariance $\Sigma$. Let

$\mathbf{X}=\begin{pmatrix}X(1)&\cdots&X(T)\end{pmatrix}\in\mathbb{R}^{n\times
T}$

be a collection of independent copies of $X$. The problem is to estimate the
true covariance $\Sigma$ out of the _empirical estimator_

$\mathbf{E}=\frac{1}{T}\mathbf{X}\mathbf{X}^{\prime}.$

We are specifically looking for a Rotationally Invariant Estimator111Such
estimators arise naturally in the Bayesian framework where $\Sigma$ has been
chosen at random with a distribution of which we only know that it is
invariant under the action of the orthogonal group by conjugation., i.e. an
estimator of $\Sigma$ which has the same eigenvectors as $\mathbf{E}$ (thus
differs from it only by the eigenvalues). We want this estimator to be
_optimal_ , i.e. to be solution of

$\operatorname{argmin}_{\operatorname{estimators}}\|\operatorname{Estimator}-\Sigma\|_{\text{F}}$
(1)

among the estimators whose eigenvectors are those of the empirical estimator
$\mathbf{E}$. Here, $\|\cdot\|_{\text{F}}$ denotes the _Frobenius norm_ , i.e.
the standard Euclidean norm on matrices:

$\|M\|_{\text{F}}:=\sqrt{\operatorname{Tr}MM^{\prime}}.$ (2)

Let $\lambda_{1},\ldots,\lambda_{n}$ denote the eigenvalues of $\mathbf{E}$,
with associated eigenvectors $\mathbf{u}_{1},\ldots,\mathbf{u}_{n}$. It is
obvious, by invariance of the Frobenius norm under the action of the
orthogonal group by conjugation, that the solution of (1) is the real
symmetric matrix with the eigenvectors $\mathbf{u}_{1},\ldots,\mathbf{u}_{n}$
and respective associated eigenvalues
$\mathbf{u}_{k}^{\prime}\Sigma\mathbf{u}_{k}$ ($1\leq k\leq n$).

It follows that the optimal RIE is the one given by the formula

$\widehat{\mathbf{E}}:=\sum_{k=1}^{n}\widehat{\lambda}_{k}\mathbf{u}_{k}\mathbf{u}_{k}^{\prime}\quad\text{for
}\widehat{\lambda}_{k}:=\mathbf{u}_{k}^{\prime}\Sigma\mathbf{u}_{k}.$ (3)

Introducing the measures

$\sum_{k=1}^{n}\delta_{\lambda_{k}}$

and

$\sum_{k=1}^{n}\mathbf{u}_{k}^{\prime}\Sigma\mathbf{u}_{k}\delta_{\lambda_{k}},$

we can express the $\widehat{\lambda}_{k}$ of (3) as Radon-Nikodym
derivatives, and, using (10), we have that for any $\varepsilon>0$ such that
$[\lambda_{k}-\varepsilon,\lambda_{k}+\varepsilon]\cap\\{\lambda_{1},\ldots,\lambda_{n}\\}=\\{\lambda_{k}\\}$,

$\widehat{\lambda}_{k}=\lim_{\eta\to
0}\frac{\int_{\lambda_{k}-\varepsilon}^{\lambda_{k}+\varepsilon}\mathfrak{Im}L(x+\mathrm{i}\eta)\mathrm{d}x}{\int_{\lambda_{k}-\varepsilon}^{\lambda_{k}+\varepsilon}\mathfrak{Im}G(x+\mathrm{i}\eta)\mathrm{d}x},$
(4)

where $L$ and $G$ are the meromorphic functions on the complex plane defined
by

$L(z):=\frac{1}{T}\operatorname{Tr}\mathbf{G}\Sigma,\qquad
G(z):=\frac{1}{T}\operatorname{Tr}\mathbf{G},$

for $\mathbf{G}=\mathbf{G}(z):=(z-\mathbf{E})^{-1}$ the resolvent matrix of
$\mathbf{E}$.

The main results of [5] can be summed-up in the following one (with main
consequence Equation (6) below):

######  Theorem 1.1.

For any fixed $z\in\mathbb{C}\backslash\mathbb{R}$, we have

$L(z)=1-\frac{1}{1-q+zG(z)}+o(1),$ (5)

where $o(1)$ denotes a random variable tending to zero in probability as $n,T$
tend to infinity in such a way that $q:=n/T$ stays bounded as well as the
operator norm of $\Sigma$ and $T/\operatorname{Tr}\Sigma$.

######  Remark 1.2.

Note that as $z$ tends to $\lambda$ in (5),

$\frac{\mathfrak{Im}L(z)}{\mathfrak{Im}G(z)}\sim\frac{-\mathfrak{Im}\frac{1}{1-q+zG(z)}}{\mathfrak{Im}G(z)}=\frac{1}{|1-q+zG(z)|^{2}}\frac{\mathfrak{Im}(zG(z))}{\mathfrak{Im}G(z)}\sim\frac{\lambda}{|1-q+\lambda
G(z)|^{2}},$

which, given $G$ is normalized by $T$ and not by $n$, is the well known RIE
formula one can find e.g. in [4]. Note that in the finite dimensional case
where this formula is meant to be applied, $G(z)$ is singular at any
eigenvalue $\lambda$, hence one has to slightly regularize the denominator,
i.e. use $G(z)$ for $z=\lambda+\mathrm{i}\eta$ with a positive spectral
resolution $\eta=T^{-\alpha}$ (any exponent $\alpha\in(0,1)$, typically
$\alpha=1/2$, should work). To sum up, for such a spectral resolution $\eta$,
the formula of the cleaned eigenvalue $\widehat{\lambda}$ associated to an
eigenvalue $\lambda$ is given by

$\widehat{\lambda}=\frac{\lambda}{|1-q+\lambda
G(\lambda+\mathrm{i}\eta)|^{2}}.$ (6)

######  Remark 1.3.

One of the advantages of this approach is that it shows that the only really
needed hypothesis, for the machinery to work, is that the mean
$\operatorname{Tr}\Sigma/n$ of the eigenvalues of the true covariance matrix
and the largest of them (i.e. the operator norm $\|\Sigma\|$ of $\Sigma$) have
the same order.

Let us now show how concentration of measure and the Stein formula for
Gaussian random vectors allow to recover the above result very easily. In the
case where the random vector $\mathbf{X}$ is not Gaussian, what follows also
works with a bit more of work, using e.g. [3, Lem. 1.13.9] instead of the
Stein formula and the McDiarmid inequality instead of Gaussian measure
concentration. The idea of the proof is the following one: one starts with the
formula $z\mathbf{G}-I_{n}=\mathbf{G}\mathbf{E}$, pass to the expectation and
expand the RHT using Stein formula for Gaussian vectors. The expansion makes
products of traces appear, each of which concentrates around its expectation
by the concentration of measure principle. At the end, we obtain the desired
relation between $G(z)$ and $L(z)$.

###### Proof.

Let

$H:=\frac{1}{T}\operatorname{Tr}\mathbf{G}\mathbf{E}.$

Using Corollary 2.2,

$\displaystyle\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\mathbf{E}$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\frac{1}{T}\sum_{t}\operatorname{Tr}\mathbf{G}X(t)X(t)^{\prime}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\sum_{t}\operatorname{\mathbb{E}}X(t)^{\prime}\mathbf{G}X(t)$
$\displaystyle=$
$\displaystyle\frac{1}{T}\sum_{t}\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma+\frac{1}{T}\sum_{t}\sum_{k=1}^{n}\operatorname{\mathbb{E}}\mathbf{e}_{k}^{\prime}\Sigma\left(\frac{\partial}{\partial
X(t)_{k}}\mathbf{G}\right)X(t)$

Note that for any $t$, at $X(1),\ldots,\widetilde{X(t)},\ldots X(T)$ fixed,
the differential, at $X(t)\in\mathbb{R}^{n}$, of the function
$X(t)\mapsto\mathbf{G}$ is the function

$x\in\mathbb{R}^{n}\mapsto\frac{1}{T}\mathbf{G}\left(X(t)x^{\prime}+xX(t)^{\prime}\right)\mathbf{G},$

so that

$\frac{\partial}{\partial
X(t)_{k}}\mathbf{G}=\frac{1}{T}\mathbf{G}\left(X(t)\mathbf{e}_{k}^{\prime}+\mathbf{e}_{k}X(t)^{\prime}\right)\mathbf{G}$

and

$\displaystyle\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\mathbf{E}$
$\displaystyle=$
$\displaystyle\frac{1}{T}\sum_{t}\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma+\frac{1}{T^{2}}\sum_{t}\sum_{k=1}^{n}\operatorname{\mathbb{E}}\mathbf{e}_{k}^{\prime}\Sigma\left(\mathbf{G}\left(X(t)\mathbf{e}_{k}^{\prime}+\mathbf{e}_{k}X(t)^{\prime}\right)\mathbf{G}\right)X(t)$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma+\frac{1}{T^{2}}\sum_{t}\sum_{k=1}^{n}\operatorname{\mathbb{E}}\left(\mathbf{e}_{k}^{\prime}\Sigma\mathbf{G}X(t)\mathbf{e}_{k}^{\prime}\mathbf{G}X(t)+\mathbf{e}_{k}^{\prime}\Sigma\mathbf{G}\mathbf{e}_{k}X(t)^{\prime}\mathbf{G}X(t)\right)$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma+\frac{1}{T^{2}}\sum_{t}\operatorname{\mathbb{E}}\left(X(t)^{\prime}\mathbf{G}\Sigma\mathbf{G}X(t)+\operatorname{Tr}(\Sigma\mathbf{G})X(t)^{\prime}\mathbf{G}X(t)\right)$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma+\frac{1}{T}\operatorname{\mathbb{E}}\operatorname{Tr}\mathbf{G}\Sigma\mathbf{G}\mathbf{E}+\frac{1}{T}\operatorname{\mathbb{E}}\operatorname{Tr}(\Sigma\mathbf{G})\operatorname{Tr}(\mathbf{G}\mathbf{E})$

Dividing by $T$ and using Proposition 2.3 with [2, Lem. B.2], we find

$H=L+LH+o(1).$ (7)

Now, we are close to the conclusion but have to divide by $1+H$. Writing
$\mathbf{G}\times(z-\mathbf{E})=I_{n}$, we have

$z\mathbf{G}-I_{n}=\mathbf{G}\mathbf{E},$

so that

$H=zG-q=\frac{q}{n}\sum_{i=1}^{n}\frac{\lambda_{i}}{z-\lambda_{i}}.$

Note that for $x=\mathfrak{Re}z$, $\eta=\mathfrak{Im}z$, we have

$|\mathfrak{Im}H|=\frac{q|\eta|}{n}\sum_{i=1}^{n}\frac{\lambda_{i}}{(x-\lambda_{i})^{2}+\eta^{2}}\geq\frac{|\eta|\operatorname{Tr}\mathbf{E}/T}{2x^{2}+2\|\mathbf{E}\|^{2}+\eta^{2}}.$
(8)

By (8) and Lemma 1.4, there is $c>0$ such that with probability tending to
$1$, $|1+H|>c$.

Hence by (7),

$L=\frac{H}{1+H}+o(1)=1-\frac{1}{1+H}+o(1).$

∎

######  Lemma 1.4.

As $n,T$ tend to infinity in such a way that $q=n/T$ stays bounded and
$\Sigma$ stays bounded in operator norm, there is a constant $C>0$ (depending
on the bounds on $q$ and on $\|\Sigma\|$) such that with probability tending
to $1$, the operator norm $\|\mathbf{E}\|$ of $\mathbf{E}$ and its trace
$\operatorname{Tr}\mathbf{E}$ satisfy

$\|\mathbf{E}\|\leq q\|\Sigma\|C\qquad\textrm{ and
}\qquad\operatorname{Tr}\mathbf{E}\geq\operatorname{Tr}\Sigma/C$

###### Proof.

Let $\mathbf{Y}\in\mathbb{R}^{n\times T}$ be a matrix with independent
standard Gaussian entries. Then $\mathbf{X}$ has the same distribution as
$\Sigma^{1/2}\mathbf{Y}$ and $T\|\mathbf{E}\|$ has the same distribution as
$\|\Sigma\mathbf{Y}\mathbf{Y}^{\prime}\|$. Thus to prove the part about the
operator norm, it suffices to prove that there is a universal constant $C_{0}$
such that with probability tending to $1$,
$\|\mathbf{Y}\mathbf{Y}^{\prime}\|\leq(n+T)C_{0}$. This is a very well known
fact, that can be obtained for example as follows: given the eigenvalues of
$\mathbf{Y}\mathbf{Y}^{\prime}$ are, up to the null eigenvalue, the squares of
the eigenvalues of $\begin{pmatrix}0&\mathbf{Y}\\\
\mathbf{Y}^{\prime}&0\end{pmatrix}$, it follows from [3, Th. 1.13.17]. The
part about the trace follows from the fact that $\operatorname{Tr}\mathbf{E}$
has the same law as
$T^{-1}\operatorname{Tr}\Sigma\mathbf{Y}\mathbf{Y}^{\prime}$ which has
expectation $\operatorname{Tr}\Sigma$ and variance $\leq
2T^{-1}\operatorname{Tr}\Sigma^{2}$. ∎

##  2\. Appendix

###  2.1. Stieltjes transform inversion

Any signed measure $\mu$ on $\mathbb{R}$ can be recovered out of its Stieltjes
transform

$g_{\mu}(z):=\int\frac{\mathrm{d}\mu(t)}{z-t},\quad
z\in\mathbb{C}\backslash\mathbb{R}$ (9)

by the formula

$\mu=-\frac{1}{\pi}\lim_{\eta\to
0^{+}}(\mathfrak{Im}g_{\mu}(x+\mathrm{i}\eta)\mathrm{d}x),$ (10)

where the limit holds in the weak topology (see e.g. [1, Th. 2.4.3] and use
the decomposition of any signed measure as a difference of finite positive
measures).

###  2.2. Stein formula for real Gaussian random vectors

######  Proposition 2.1.

Let $X=(X_{1},\ldots,X_{d})$ be a real centered Gaussian vector with
covariance $\Sigma$ and $f:\mathbb{R}^{d}\to\mathbb{R}$ be a $\mathcal{C}^{1}$
function with gradient having at most polynomial growth at infinity. Then for
all $i_{0}=1,\ldots,d$,

$\operatorname{\mathbb{E}}X_{i_{0}}f(X_{1},\ldots,X_{d})\;=\;\sum_{k=1}^{d}\Sigma_{i_{0}k}\operatorname{\mathbb{E}}(\partial_{k}f)(X_{1},\ldots,X_{d}).$

######  Corollary 2.2.

With the same notation, considering $X$ as a column vector, for
$F:\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ a matrix-valued function with
gradient having at most polynomial growth at infinity. we have

$\displaystyle\operatorname{\mathbb{E}}X^{\prime}F(X)X$ $\displaystyle=$
$\displaystyle\operatorname{Tr}\Sigma\operatorname{\mathbb{E}}F(X)+\sum_{k=1}^{d}\left(\operatorname{\mathbb{E}}\Sigma(\partial_{k}F)(X)X\right)_{k}.$
(11)

###### Proof.

We have, by Proposition 2.1,

$\displaystyle\operatorname{\mathbb{E}}X^{\prime}F(X)X$ $\displaystyle=$
$\displaystyle\sum_{ij}\operatorname{\mathbb{E}}X_{i}X_{j}F(X)_{ij}$
$\displaystyle=$
$\displaystyle\sum_{ijk}\operatorname{\mathbb{E}}\Sigma_{ik}\frac{\partial}{\partial
X_{k}}X_{j}F(X)_{ij}$ $\displaystyle=$
$\displaystyle\sum_{ijk}\operatorname{\mathbb{E}}\Sigma_{ik}\left(\delta_{j=k}F(X)_{ij}+X_{j}(\partial_{k}F)(X)_{ij}\right)$
$\displaystyle=$
$\displaystyle\operatorname{Tr}\Sigma\operatorname{\mathbb{E}}F(X)+\sum_{ijk}\operatorname{\mathbb{E}}\Sigma_{ik}X_{j}(\partial_{k}F)(X)_{ij}$
$\displaystyle=$
$\displaystyle\operatorname{Tr}\Sigma\operatorname{\mathbb{E}}F(X)+\sum_{k}\left(\operatorname{\mathbb{E}}\Sigma(\partial_{k}F)(X)X\right)_{k}$

∎

###  2.3. Concentration of measure for Gaussian vectors

The following proposition can be found e.g. in [1, Sec. 4.4.1] or [6, Th.
5.2.2].

######  Proposition 2.3.

Let $X=(X_{1},\ldots,X_{d})$ be a standard real Gaussian vector and
$f:\mathbb{R}^{d}\to\mathbb{R}$ be a $\mathcal{C}^{1}$ function with gradient
$\nabla f$. Then we have

$\operatorname{Var}(f(X))\;\leq\;\operatorname{\mathbb{E}}\|\nabla
f(X)\|^{2},$ (12)

where $\|\,\cdot\,\|$ denotes the standard Euclidean norm.

Besides, if $f$ is $k$-Lipschitz, then for any $t>0$, we have

$\mathbb{P}(|f(X)-\operatorname{\mathbb{E}}f(X)|\geq
t)\;\leq\;2\mathrm{e}^{-\frac{t^{2}}{2k^{2}}},$ (13)

i.e. $f(X)-\operatorname{\mathbb{E}}f(X)$ is Sub-Gaussian with Sub-Gaussian
norm $\leq k$, up to a universal constant factor.

## References

  * [1] Anderson, G., Guionnet, A., Zeitouni, O. _An Introduction to Random Matrices_. Cambridge Studies in Advanced Mathematics, 118 (2009).
  * [2] Benaych-Georges, F., Couillet, R. _Spectral analysis of the Gram matrix of mixture models_ , ESAIM Probab. Statist., Vol. 20 (2016), 217–237.
  * [3] Benaych-Georges, F., Knowles, A. _Local semicircle law for Wigner matrices. Advanced topics in random matrices_ , 1–90, Panor. Synthèses, 53, Soc. Math. France, Paris, 2017.
  * [4] Bun, J., Bouchaud, J.-P., Potters, M. _Cleaning large correlation matrices: Tools from Random Matrix Theory_ , Physics Reports Volume 666, Review article, 1–109, 2017.
  * [5] Ledoit, O., Péché, S. _Eigenvectors of some large sample covariance matrix ensembles_ Probability Theory and Related Fields, 2011, 151.1, 233–264
  * [6] Vershynin, R. _High-Dimensional Probability: An Introduction with Applications in Data Science_ , Cambridge Series in Statistical and Probabilistic Mathematics, 2018.


# Experimental study of submerged liquid metal jet in transverse magnetic
field

Ivan A. Belyaev
Joint Institute for High Temperature RAS,
Izhorskaya 13 Bd.2, 125412, Moscow, Russia
<EMAIL_ADDRESS>
&Ivan S. Mironov
Joint Institute for High Temperature RAS,
Izhorskaya 13 Bd.2, 125412, Moscow, Russia
&Nikita A. Luchinkin
Joint Institute for High Temperature RAS,
Izhorskaya 13 Bd.2, 125412, Moscow, Russia
& Yaroslav I. Listratov
Moscow Power Engineering Institute,
Krasnokazarmennaya 14, 111250, Moscow, Russia
& Yuri B. Kolesnikov
Technische Universitat Ilmenau,
PF 100565, 98684 Ilmenau,Germany
&Dmitry Kransnov
Technische Universitat Ilmenau,
PF 100565, 98684 Ilmenau,Germany
& Oleg Zikanov
University of Michigan Dearborn,
Dearborn, MI 48128-1491, USA
Sergey Molokov
Technische Universitat Ilmenau,
PF 100565, 98684 Ilmenau,Germany

(November 05, 2021)

###### Abstract

A liquid metal flow in the form of a submerged round jet entering a square
duct in the presence of a transverse magnetic field is studied experimentally.
A range of high Reynolds and Hartmann numbers is considered. Flow velocity is
measured using electric potential difference probes. A detailed study of the
flow in the duct’s cross-section about seven jet’s diameters downstream of the
inlet reveals the dynamics, which is unsteady and dominated by high-amplitude
fluctuations resulting from instability of the jet. The flow structure and
fluctuation properties are largely determined by the value of the Stuart
number N. At moderate N, the mean velocity profile retains a central jet with
three-dimensional perturbations increasingly suppressed by the magnetic field
as N grows. At higher values of N, the flow becomes quasi-two-dimensional and
acquires the form of an asymmetric macrovortex, with high-amplitude velocity
fluctuations reemerging.

##  1 Introduction

Fusion is considered to be a future source of clean energy. It does not
involve nuclear waste characteristic for fission. The fuel for fusion, a
combination of the isotopes of hydrogen deuterium and tritium, are either
abundant in the sea water or can be generated (bred) as a byproduct of the
reaction. In the reactor, these isotopes are heated up to the temperature of
several million Kelvin until the fusion reaction starts. The burning plasma in
the reactor is contained by a very high magnetic field of 5-10T. The plasma
chamber is surrounded by the blankets, which may employ liquid metals, Li or
PbLi, and which act simultaneously as heat exchangers diverting reaction’s
energy into an external circuit, shields, and breeders of tritium (Bühler,
2007; Abdou et al., 2015). Since liquid metals are electrically conducting,
their flow is strongly affected by the magnetohydrodynamic (MHD) interaction.

Liquid metal blankets contain many elements, such as straight ducts, bends,
manifolds, etc. One of the key elements of the blanket is an expansion of a
circular pipe into a larger rectangular duct (see figure 1). Understanding the
flow dynamics in the expansion is important because it largely contributes to
the flow distribution in a subsequent manifold. If both the ratio of the duct
to pipe dimensions and the fluid velocity are high, the resulting flow is in
effect an initially circular jet entering a rectangular duct under the
influence of a high, transverse magnetic field.

Another important technology, in which a similar flow structure is observed is
the electromagnetic braking in continuous casting of steel. Molten steel
entering the solidification mold through a nozzle forms a submerged jet in a
rectangular domain. As explained, for example by Cukierski and Thomas (2008),
a strong transverse magnetic field is often applied to suppress excessive
turbulent fluctuations around the jet and prevent strong impingement of the
jet into the wall.

The flow is also interesting from the fundamental point of view. It can be
classified as one of several archetypal configurations of liquid metal MHD,
other such configurations being, for example, a duct flow in a uniform or non-
uniform transverse magnetic field, an internal shear layer, an MHD boundary
layer, or magnetoconvection in an enclosure or a duct. As we will see in this
paper, the jet flow demonstrates and helps to understand fundamental
properties of flow transformation caused by an imposed magnetic field. In
particular, the flow remains turbulent even in a high magnetic field as it
evolves from three- to quasi-two-dimensional form.

Instability and transition to turbulence in axisymmetric and planar jet of
viscous fluids without the magnetic field effect is a thoroughly studied area
of fluid mechanics (see, e.g., Michalke (1984); Thomas and Goldschmidt (1986);
Stanley et al. (2002)). Jets are known to be unstable at all, but smallest
values of the Reynolds number (Sato, 1960; Sato and Sakao, 1964). The
instability mechanism is that of the Kelvin-Helmholtz type and leads to
vortices in strongly sheared outer regions of a jet and jet’s rapid distortion
and development of turbulence. The process is strongly affected by inlet
conditions (shape of the inlet, velocity profile, and presence of
perturbations).

A sufficiently strong imposed magnetic field completely changes jet flow
behaviour if fluid is electrically conducting. An example of this is the
experiments with a submerged annular jet in an axial uniform magnetic field
(Kalis and Kolesnikov, 1984; Klyukin and Thess, 1993). Here, the transition to
turbulence also occurs through the Kelvin-Helmholtz type instability, but the
process demonstrates unique MHD features, such as development of a quasi-two-
dimensional state and pronounced hysteresis in the plane of the Reynolds and
Hartmann numbers (Kljukin and Kolesnikov, 1989).

Figure 1: Configuration of the experiment: the test section and the
$x$-distribution of the main, $B_{z}$, component of the imposed magnetic
field. The distances are given in mm, while the magnetic field - in Tesla. An
inlet pipe of length 1 m is necessary to position the experimental section in
the magnet.

Flows of laminar, inviscid liquid metal jets in an unbounded fluid in a
transverse, uniform magnetic field $\mathbf{B}$ were theoretically studied by
Moffatt and Toomre (1967) and Davidson (1995). MHD effects in these flows can
be understood by considering the electrodynamic part of the MHD equations
(Müller and Bühler, 2001), namely Ohm’s law

$\mathbf{j}=\sigma(-\nabla\varphi+\mathbf{u}\times\mathbf{B}),$ (1)

and the charge conservation equation

$\nabla\cdot\mathbf{j=0}.$ (2)

Here $\mathbf{u}=(u_{x},u_{y},u_{z})$ is the fluid velocity,
$\mathbf{j}=(j_{x},j_{y},j_{z})$ is the electric current density, $\varphi$ is
the electric potential, $\sigma$ is the electrical conductivity of the fluid,
and $(x,y,z)$ are the Cartesian co-ordinates. The Lorentz force acting on the
fluid is:

$\mathbf{f}=\mathbf{j}\times\mathbf{B}.$ (3)

Moffatt and Toomre (1967) considered a two-dimensional ($y$-independent) jet
flowing in the $x-$direction in a transverse magnetic field
$\mathbf{B}=B\mathbf{e}_{z}$. In a two-dimensional flow $\varphi=$ constant,
so that the first term in the brackets in (1) vanishes. Then
$\mathbf{j}=\sigma\mathbf{u}\times\mathbf{B}=-\sigma u_{x}B\mathbf{e}_{y}$,
the Lorentz force, $\mathbf{f}=-\sigma u_{x}B^{2}\mathbf{e}_{x}$, reduces to
pure braking at every point within the flow, and the axial component of
momentum is annihilated by the field.

(a) |  |  (b)
---|---|---

Figure 2: Schematic diagram of the jet expansion along the magnetic field
lines: an initially round infinite jet (broken line), electric current lines,
and the directions of the Lorentz force (a), evolution of the jet profiles in
the centre-plane $y=0$ as time increases (b).

If the jet is not two-dimensional, axial momentum is conserved in the zero-
viscosity limit, but is redistributed along the field lines by the Lorentz
force (Davidson, 1995) as shown in figure 2. To understand the process we
consider a straight, infinite, initially circular jet flowing in the
$x-$direction. Now $\varphi$ is a function of both $y$ and $z$. In the
original jet cross-section, currents are induced by the electromotive force
$\mathbf{u}\times\mathbf{B}$, and the resulting Lorentz force brakes the flow.
The currents induce the electric potential difference between the left and
right parts of the jet shown by ${}^{\prime}+^{\prime}$ and
${}^{\prime}-^{\prime}$ in figure 2a. This difference in the electric
potential drives the current back through the fluid from
${}^{\prime}+^{\prime}$ to ${}^{\prime}-^{\prime}$ as shown in the figure. The
Lorentz force accelerates the fluid above and below the original jet. The
process continues with time as more and more layers of fluid become involved
in the motion. The development of the axial velocity profiles with time are
shown schematically in figure 2b. The resulting jet has an anisotropic
structure. We will show below that this evolution is very important for
understanding of the flow pattern. It should be noted that similar anisotropic
structures are found in many MHD flows, including flows around solid bodies
(see, e.g., Ludford, 1960; Hasimoto, 1960), rotation of vortices (Davidson,
1995) and solid bodies (Molokov, 1993), electrically driven flows (Moffatt,
1964), etc.

It is evident that the just outlined scenario of jet transformation does not
form a complete picture. Inertial effects, hydrodynamic instabilities, and
turbulence anticipated at sufficiently high Reynolds numbers are not included.
Not less importantly, the effect of walls surrounding the flow domain and
altering paths of electric currents is expected to be strong, much more so
than in a flow without magnetic field. The effect remains poorly understood.

Surprisingly, very little work has been done to analyze the flow. To our best
knowledge, no experiments have been performed in the last three decades. The
only computational work addressing the conclusions of Davidson (1995)
concerning the transformation of a round jet by a transverse magnetic field
was the large eddy simulations of Kim and Choi (2004). The closest we can find
to a study of MHD jet in a wall-bounded domain is the recent work of Krasnov
et al. (2021) and Belyaev et al. (2021), in which transformation of a planar
jet in a transverse magnetic field is analyzed.

The aim of this study is to understand the dynamics of a round jet as it
enters a square duct. The study is divided into two parts. Here, in Part 1, we
present the experimental results. Theoretical results using Direct Numerical
Simulation will be presented in Part 2 by Krasnov et al. (2022).

##  2 MHD transformation of a submerged jet in a duct

MHD flows of liquid metals in ducts are governed by four dimensionless
parameters: the Reynolds number, Re, the Hartmann number, $Ha$, the Stuart
number, $N=Ha^{2}/\mbox{{Re}}$, also known as the interaction parameter, and
the wall conductance ratio, $c$. They characterise the ratios of inertial to
viscous, electromagnetic to viscous, electromagnetic to inertial forces, and
the ratio of wall to fluid electrical conductances, respectively. Their
definitions are given in Table 1.

Consider a flow in a duct (figure 3a). The pairs of walls transverse and
parallel to the magnetic field are called Hartmann walls and the sidewalls,
respectively. In a sufficiently strong magnetic field Hartmann layers of non-
dimensional thickness $\sim Ha^{-1}$ are formed at the Hartmann walls. The
streamwise velocity has an exponential profile inside these layers as
determined by the balance between the viscous and Lorentz forces. There are
also Shercliff (sidewall) layers of thickness $\sim Ha^{-1/2}$ forming at the
sidewalls, but they seem to play no role in this study.

Consider now a flow created specifically by a round jet entering duct. From
the discussion in the previous section it is clear that it will be transformed
by the action of the electromagnetic forces. The transformation is in space,
with the flow evolving with the downstream distance. Initially, the electric
current lines are similar to those shown in figure 2a. As the jet progresses
through the duct, the gradients of the velocity profile along the magnetic
field become smaller, and the jet develops anisotropy. Further downstream, the
gradients virtually disappear. The velocity profile becomes almost independent
of the $z-$coordinate except for thin Hartmann layers. The jet becomes quasi-
two-dimensional, planar-like at a distance $L$ from the inlet as illustrated
in figure 3. If Hartmann walls are electrically insulating, all the lines of
the electric current will pass through these layers. They are very thin, and
their electrical resistance is inversely proportional to their thickness, i.e.
scales as $Ha$. Because of the increased resistance, the induced currents
become much smaller in magnitude, and the resulting Hartmann damping
decreases. If Hartmann walls are electrically conducting, part of the current
will flow in the walls as well. Thus its magnitude increases, and damping
increases as a result.

In a three-dimensional flow, the longitudinal component of electric current
$j_{x}$ also appears due to the variation of the electric potential along the
flow (Müller and Bühler, 2001). Indeed, in the pipe the average fluid velocity
is higher than in the duct. Thus, $\mathbf{u}\times\mathbf{B}$ is higher as
well, and the induced potential difference in the $(x,y)-$planes is higher.
These currents result in the Lorentz force stretching the jet in the
$(x,y)-$planes towards the sidewalls, but the importance of this effect in
this experimental study is unclear.

(a) |  |  (b)
---|---|---

Figure 3: Two-dimensionalization of a round jet by a transverse magnetic field
in a duct.

##  3 Description of experiment

Experiments are performed using the HELME facility (HELMEF) described by
Belyaev et al. (2017). Mercury circulating in a closed loop is used as a model
liquid. Walls are made of stainless steel and have thickness of 2 mm. The
test-section is placed inside of an electromagnet, which can create the
magnetic field of up to $1.8$ Tesla.

The test section is configured as illustrated in figure 1. Mercury enters
through an inlet into a square duct of half-height $a=28$ mm and length
$L_{x}=500$ mm. The inlet has the form of a round pipe, 90 mm long with 6 mm
inner diameter located at the mid-plane of the duct (figure 1a). The flow is
driven in the $x$-direction by a pump. An imposed magnetic field $\bm{B}$ has
the dominant component along the $z$-direction. The uniform part of the
magnetic field covers approximately the first $250$ mm of the duct length.

In order to avoid possible disturbances caused by gas bubbles in the test
section, the loop is vacuumed prior to filling it with mercury. Furthermore,
experiments are repeated, with consistent results, with the test section
oriented both vertically and horizontally.

A range of flow rates from 0.75 l/min to 5.4 l/min corresponding to the duct
Reynolds numbers $600\leq\mbox{{Re}}\leq 7000$ is covered in the experiments.
The bottom limit is set by the impossibility of accurate measurements of flow
velocity occurring at lower flow rates due to the significant impact of
electromagnetic noise. The top limit corresponds to the beginning of
cavitation in the jet area. The mean flow velocity $U_{0}$ used to define the
duct Reynolds number is based on the flow rate measured using an
electromagnetic flow meter. Flow stabilizer with gas cavity is used to avoid
possible flow fluctuations produced by the pump.

Further parameters of the experiment are given in Table 1. The values of the
physical properties of mercury, such as the kinematic viscosity $\nu$ and the
electrical conductivity $\sigma$ at room temperature (about 18 oC) are taken
from Kirillov (2007). The other physical parameters used in the table are the
mean flow velocities $U_{o}$ in the duct and $U_{p}$ in the inlet pipe, and
the thickness $d_{w}$ and electrical conductivity $\sigma_{w}$ of the walls.

Parameter | Quantity | Value
---|---|---
Duct inner section |  $2a$ x $2a$, mm x mm | $56\times 56$
Pipe radius |  $r$, mm | $3$
Characteristic length |  $a$, mm | 28
Magnetic induction |  $B$, T | $0-1.7$
Length of uniform MF |  $L_{B}$, m | $\approx 0.9$
Duct Reynolds number | $\mbox{{Re}}\equiv U_{o}a/{\nu}$ | $600-7000$
Duct Hartmann number | $\mbox{{Ha}}\equiv Ba\sqrt{\sigma_{l}/{\rho\nu}}$ | $0-1300$
Duct Stuart number | $\mbox{{N}}\equiv\mbox{{Ha}}^{2}/\mbox{{Re}}$ | $0-2100$
Pipe Reynolds number | $\mbox{{Re}}_{p}\equiv U_{p}r/{\nu}=24\mbox{{Re}}$ | $14000-82000$
Pipe Hartmann number | $\mbox{{Ha}}_{p}\equiv Br\sqrt{\sigma_{l}/{\rho\nu}}=0.1\mbox{{Ha}}$ | $0-140$
Pipe Stuart number | $\mbox{{N}}_{p}\equiv\mbox{{Ha}}_{p}^{2}/\mbox{{Re}}_{p}$ | $0-1$
Wall conductance ratio | $c\equiv\sigma_{w}d_{w}/\sigma a$ | $\lesssim 0.01$
Table 1:  Physical and non-dimensional parameters of the experiment.

A comment is in order concerning the value of the wall conductance ratio
reported in Table 1. While correct from the formal point of view, the formula
used to define $c$ may provide results quite different from the effective
value observed in the experiment. During the experimental program, it is
impossible to directly assess the electrical conductivity of the wall. No
special methods are used to achieve good contact between mercury and the
stainless steel wall other than vacuuming water vapor. It is well known that
oxides and other deposits worsen the electrical contact with the wall and
decrease the relative electrical conductivity by one to two orders of
magnitude. Therefore, instead of applying the formal definition of $c$ we
estimate the range of the possible effective values of the wall conductance
ratio $c$ in the experiments as between 0.001 and 0.01. The influence of the
wall conductance ratio on the flow is explored via numerical simulations and
further discussed by Krasnov et al. (2022).

Flow velocity is measured using the electric potential difference method (see,
e. g., Mistrangelo and Bühler, 2010). A 4-electrode potential sensor
illustrated schematically in figure 4c can measure two components of velocity
at a point. With such a sensor, it is possible to simultaneously measure the
electric potential $\varphi$ at four electrodes, and then to calculate the
local velocity components using the projection of Ohm‘s law onto the $x-$ and
$y-$ axes. Neglecting electric currents, as their magnitude is small in ducts
with poorly conducting walls, and replacing the partial derivatives with their
difference approximations gives:

$u_{x}\approx\frac{\varphi_{a}-\varphi_{b}}{c_{1}B_{0}\Delta
l_{y}}\approx\frac{\varphi_{c}-\varphi_{d}}{c_{1}B_{0}\Delta l_{y}},$ (4)

$u_{y}\approx\frac{\varphi_{c}-\varphi_{a}}{c_{1}B_{0}\Delta
l_{x}}\approx\frac{\varphi_{b}-\varphi_{d}}{c_{1}B_{0}\Delta l_{x}},$ (5)

where $\varphi$ is the electric potential at the electrodes, $\Delta
l_{y}=2$mm, $\Delta l_{x}=2-4$mm are the distances between the corresponding
electrodes, and $c_{1}$ is the coefficient depending on the conditions of
closing the induced currents. For the Hartmann number values used in the
experiment, $c_{1}\approx 1$.

Sensors of two principal designs were used in the majority of the reported
experiments. The sensor of type 1 (see figure 4a) was used during the first
stage of the experimental program. Its design was primarily motivated by the
desire to avoid damage of the sensor by a strong jet flow. The sensor of type
2 (see figure 4b) was built and utilized at later stages of the program in
order to validate the results obtained with the sensor of type 1 and to
minimize the potential influence of the sensor on measurement data obtained in
situations with reversed flow. Yet another, modified sensor of design similar
to type 2 was utilized to investigate the jet development far downstream (in
the cross-section located at the distance $5.2a$ from jet inlet).

Consistently good agreement between measurements obtained using different
sensors was found in all the experiments. An illustration of that is provided
in figure 5b. The observed small deviations are attributed to uncertainty of
sensor’s positioning. Most of the results reported in this paper are based on
the data obtained with the sensor of type 2, except for the diagram in the
last figure, where the results of the experiments performed using sensors of
all three types are summarised.

|  (a) |
---|---|---

|  (b) |
---|---|---

(c)
---

Figure 4:  Sensors of different construction: a - Type 1; b - Type 2 c -
Schematic representation of the principle of measurements. 1 – tin contact
layer; 2 – copper electrode; 3 – insulating layer of epoxy.

The electrode size of 0.3 mm was chosen to avoid possible errors caused by
sensor damage during measurements near walls. Sensors were installed on
scanning probes extended in the direction opposite to the flow. The probes
allowed us to either record the velocity signal at one location over time or
map the mean velocity in a cross-section of the duct by taking time-averaged
measurements point by point (see figure 5a for the measurement points utilized
for the mapping). Simultaneous measurement of velocities at more than one
point was impossible.

(a) |  |  (b)
---|---|---

Figure 5: (a) Measurement points used to map distributions of time-averaged
velocity in duct’s cross-section. The inlet shape of the jet is shown by a
circle. (b) Profiles of time-averaged streamwise velocity measured by sensors
of type 1 and type 2 at $x=44\textrm{mm}=1.6a$, $z=0$ for
Re=800$,\mbox{{Re}}_{p}=19\cdot 10^{3}$, $Ha=200$. Time averaging is over
about 600 s.

(b) |  |  (c)
---|---|---

Figure 6: Sensor calibration using flow in a uniform duct of rectangular
cross-section 56x16 mm. The range of mean velocities $U_{0}$ and the magnetic
field strengths correspond to the values of the Reynolds and Hartmann numbers
based on $U_{0}$ and the half of the longer side of the duct between 9000 and
55000 and between 80 and 1150, respectively. (a) - Duct configuration and
points at which velocity profile is measured; (b) - Velocity values at the
duct center versus magnetic field; (c) - Velocity profiles at different flow
rates and Ha = 1100. Visible distortion at $-6<y-4$ mm is caused by
deformation of the probe from wall contact.

Time-averaged streamwise velocity values, $<u_{x}>$ and $<u_{y}>$, and
standard deviations $\sigma_{x}$ and $\sigma_{y}$ are used in further
analysis. Standard deviations are calculated as:

$\sigma_{x}=\sqrt{\sum_{j=0}^{n-1}\frac{\left({u_{x}}_{j}-U_{x}\right)^{2}}{n-1}},\quad\sigma_{y}=\sqrt{\sum_{j=0}^{n-1}\frac{\left({u_{y}}_{j}-U_{y}\right)^{2}}{n-1}},$
(6)

where $n$ is the number of measured points in the velocity signal.

The inherent limitation of the electric potential sensors is their inability
to measure velocity in boundary layers. This is due to the facts that the
thickness of MHD boundary layers is significantly smaller than the distance
between the electrodes, and that positioning the sensor very close to the wall
is impossible because of the danger of damaging the sensor and increased
uncertainty of measurements due to possible wall-sensor contact.

In the rest of this section we discuss the factors potentially negatively
affecting the accuracy of electric potential sensors. We also list measures
taken during the experiments to identify and minimize the effect of these
factors. It must be stressed that a consistent moderate shift in measured
absolute values of velocity is not detrimental for achieving the goals of our
study. As we will see below, the flow transformation caused by the magnetic
field is described in terms of relative changes of velocity profiles or in
terms of centered statistical moments of velocity fluctuations. These
characteristics are insensitive to a shift of absolute velocity values,
provided the shift is the same in all the experiments.

One potentially negative factor is the thermoelectrical effects associated
with the contact between dissimilar metals in a non-uniform temperature field.
To minimize it, the mercury loop was water-cooled, so that temperature
remained constant and uniform at 291 degrees K. Variations of temperature in
the test section were within 0.4 K.

Another factor is the error (zero shift) caused by electromagnetic noise in
measurement equipment. To avoid this effect, each experimental session started
with measurements performed at zero imposed magnetic field. The registered
small deviation from the expected zero value of velocity was then used for
correction of measurements in production runs performed at non-zero imposed
magnetic field. Furthermore, two acquisition systems were used, one based on
an NI cDAQ 9213 measuring module with 100Hz antialiasing analog filtering, and
another on NI 4071 digital multimeters with multi-speed digitization and
digital processing of 1kHz signals. After the aforementioned correction, the
two systems showed consistent results.

The biggest potential factor negatively affecting measurement accuracy is the
disturbance of the flow by the probe. The disturbance may happen directly,
with the probe modifying the velocity as a solid obstacle, or indirectly via
MHD effects caused by the flow around the probe. One already mentioned step
undertaken in the experiments to reduce the probe’s influence was to always
insert it in the direction opposite to the mean flow. To evaluate the
remaining effects, a special calibration experiment illustrated in figure 6
was carried out. The geometry of the test section was replaced by that of a
long straight duct of rectangular cross-section $56\times 16$ mm. A uniform
magnetic field was imposed along the $z$-axis. Time-averaged streamwise
velocity $\langle u_{x}\rangle$ was measured along the $y$-axis in the middle
of the duct at the distance of 600 mm downstream from the entrance into the
magnetic field. It is anticipated in such a system that a moderate magnetic
field would produce a flow with experimentally detectable sidewall and the
centerline value of $\langle u_{x}\rangle$ exceeding the mean velocity
$U_{0}$. At stronger magnetic fields, the sidewall boundary layers would be
undetectable thin, and the measured velocity would show a nearly flat profile
along the $y$-axis, with $\langle u_{x}\rangle\sim U_{0}$. The results of
experiments shown in figure 6 agree with the expectations. An exception is
flows with weak magnetic field (approximately Ha below 180), in which case
accuracy of the sensor is expected to be low. Figure 6 illustrates the
agreement by the variation of the ratio $\langle u_{x}\rangle/U_{0}$ with the
magnetic field strength (in figure 6b) and by full measured velocity profiles
at a strong magnetic field (in figure 6c). We conclude that the accuracy of
the sensor at all but the smallest Hartmann numbers considered in our
experiments is confirmed by the test.

The size of the sensor does not allow us to measure velocity distribution
within the boundary layers and, therefore, unambiguously identify the velocity
profile near the walls. It is, however, visible, that the velocity in figure
6c does not show the tendency toward an M-shaped profile typical for a flow in
a duct with thin electrically conducting walls (Müller and Bühler, 2001). This
indicates a very low value of effective conductance rate $c$ consistent with
our estimates in Table 1.

All the results reported in the following sections of the paper satisfy three
criteria: _(i)_ measurements of the time-averaged values of $u_{x}$ are
consistent with $U_{0}$ measured by an electromagnetic flowmeter, _(ii)_
nearly identical results are obtained in experiments performed with different
acquisition systems, and _(iii)_ measurements are performed at magnetic fields
high enough to ensure the stable performance of the sensor demonstrated in
calibration experiments. The experiments, in which these conditions are not
satisfied (this happens, for example in cases with low magnetic field) are
considered insufficiently accurate and not reported.

##  4 Results

Significant unsteadiness of the flow, with high-amplitude fluctuations of
velocity recorded downstream of the jet inlet, is found in all the experiments
reported in this paper. It must be stressed that the unsteadiness cannot be a
result of transition to turbulence within the duct itself. As first identified
by Murgatroyd (1953) and later discussed on the basis of experimental data by
Branover (1978) and from the theoretical perspective by Zikanov et al. (2014),
a transition or laminarization occurs in laboratory flows when
$\mbox{{Re}}/\mbox{{Ha}}$ traverses the range between 200 and 400. All the
experiments discussed below have smaller values of $\mbox{{Re}}/\mbox{{Ha}}$.
The ratio $\mbox{{Re}}_{p}/Ha_{p}=240\mbox{{Re}}/\mbox{{Ha}}$ determining the
possibility of transition within the inlet pipe can be in or above the
transition range, so some turbulent fluctuations of moderate amplitude are
possible at the jet inlet. We conclude that the high-amplitude velocity
fluctuations recorded in the experiments can only be results of the
instability of the jet itself, and that the instability is likely to be
affected by perturbations at the jet’s inlet.

Presence of high-amplitude fluctuations creates a challenge for experimental
characterization of the flow. In particular, time series of long duration (up
to 600 seconds) are required to accurately measure time-averaged velocity at
each point.

Measurements are made in the cross-section $x=44$ mm $=1.6a$, chosen so that
the jet has already developed its unsteady structure resulting from the
instability but has not strongly decayed yet. A distribution of time-averaged
velocity in the duct’s cross-section is found by positioning the sensors along
the lines $z=0$, $z=-a/2$, and $z=a/2$ perpendicular to the magnetic field,
i.e. from one sidewall to another (see figure 5a) and recording a full time
series at each measurement point.

The results for the time-averaged velocity components $\langle u_{x}\rangle$
and $\langle u_{y}\rangle$ are presented in figures 5, 7, 8 and 10. The
tendency towards formation of quasi-two-dimensionality in an initially round
jet is illustrated in figure 7 by velocity profiles at $\mbox{{Re}}=800$. At
$Ha=200$ there is already strong flow at $z=\pm a/2$, but significant
difference between the profile at the center of the duct $z=0$ and the
profiles measured at $z=-a/2$ and $a/2$ is visible (see figure 7a). At
$Ha=1300$, the three profiles are close to each other demonstrating quasi-two-
dimensionality (see figure 7b).

Figure 7 also illustrates existence of two clearly different states of the
flow. The time-averaged velocity profiles are symmetric with respect to the
centerline $y=0$ at low Hartmann numbers, such as $\mbox{{Ha}}=200$ in figure
7a. Strong asymmetry is observed at higher Ha, such as $Ha=1300$ in figure 7b.

(a) |  |  (b)
---|---|---

Figure 7: An illustration of two-dimensionalization of the velocity field
caused by magnetic field. Time-averaged longitudinal velocity profiles at
$x=44$ mm $=1.6a$ are shown for three values of $z$ in the flow at
$\mbox{{Re}}=800$, $\mbox{{Re}}_{p}=19\cdot 10^{3}$. $\mbox{{Ha}}=200$ in (a).
$\mbox{{Ha}}=1300$ in (b).

The two typical behaviours observed at moderate and high values of Ha are
further illustrated in figure 8 on the example of flows at $\mbox{{Re}}=800$,
$\mbox{{Ha}}=200$ and $\mbox{{Re}}=800$, $\mbox{{Ha}}=1300$. Time-averaged
values and waveforms of longitudinal $u_{x}$ and transverse $u_{y}$ velocity
components measured at points distributed along $y$ at $z=0$ are provided. The
waveforms show instantaneous simultaneously measured local values of $u_{x}$
and $u_{y}$, from which time-averaged values computed at the same $y$ are
subtracted. One should keep in mind that waveforms corresponding to different
measurement points are recorded during different periods of time, so any
apparent correlation between them is purely coincidental.

Figure 8:  Longitudinal $u_{x}$ (blue curves) and transverse $u_{y}$ (red
curves) velocity components measured at points located along the line $z=0$ at
$x=1.6a=44mm$ in flows with $\mbox{{Re}}=800,\mbox{{Ha}}=200$ ((a) and (b))
and $\mbox{{Re}}=800,\mbox{{Ha}}=1300$ ((c) and (d)). Profiles of $\langle
u_{x}\rangle$ and $\langle u_{y}\rangle$ obtained by averaging over 610 s are
shown in (a) and (c). Velocity fluctuations around the time-averaged values
are shown in (b) and (d). In the latter, each pair of red and blue curves
shows a velocity signal measured at a $y$-location indicated by a respective
point in the left-hand side graph. The velocity scale in (b) and (d) is such
that the half-distance between the black lines corresponds to $60U_{0}$.

Figures 8a,b illustrate a typical flow observed at moderate values of Ha. The
profile of $\langle u_{x}\rangle$ is symmetric. It has a peak at $y=0$ and
shows positive values of $\langle u_{x}\rangle$ at about $-10<y<10$ mm. This
area coincides with the area of high-amplitude fluctuations in the waveforms.
We conclude that the interval $-10<y<10$ mm corresponds to the zone swept by
an oscillating central jet. Outside this interval, i.e. at $y<-10$ mm or
$y>10$ mm, $\langle u_{x}\rangle<0$ and the amplitude of velocity fluctuations
is significantly lower. These areas are identified as those of a weakly
turbulent reversed flow. The picture of an expanding central jet is further
confirmed by the profile of $\langle u_{y}\rangle$ shown in figure 8a. After a
correction for a zero shift we observe a flow in the positive $y$-direction at
$y>0$ and in the negative $y$-direction at $y>0$.

The waveforms in figure 8b show strong fluctuations in the central area,
especially at $y=0$. The maximum amplitude exceeds $60U_{0}$ at some moments.
The dominant frequencies correspond to the time period of about 6-7 s. There
are also oscillations of much higher frequency and lower amplitude
corresponding to three-dimensional turbulence. Another interesting feature is
the correlation between fluctuations of $u_{x}$ and $u_{y}$ observed in the
side regions of the jet, but not in its center. Positive peaks of $u_{x}$
correlate with positive or negative peaks of $u_{y}$ at, respectively, $y>0$
or $y<0$. Such events are identifiable as excursions of the jet into the areas
of positive or negative $y$.

An entirely different structure of the flow is presented by the measurements
made at $\mbox{{Re}}=800$, $\mbox{{Ha}}=1300$ (see figures 8c,d). The
asymmetry, which we have already presented in figure 7, demonstrates that the
jet is shifted toward one side of the duct. It must be stressed that the
apparent shift is a result of time-averaging. The signal of $u_{x}$ at $y=a/2$
shown in figure 9 illustrates the fact that velocity fluctuates and may become
negative at some moments. Once time-averaging over a sufficiently long time
period is taken, however, the jet is found consistently located on one side of
the duct, at positive (as in figures 10b and 8c) or negative $y$. No
preference for one side and no spontaneous switches of the time-averaged
profiles from one side to another have been found in the experiments. A
possibility of such behaviours cannot be excluded, since no special studies
that could detect them (statistical studies with multiple realizations of the
same flow or very long measurement runs) have been completed.

Figure 9: Signal of longitudinal velocity $u_{x}$ measured at $z=0$, $y=a/2$
in the flow with $\mbox{{Re}}=900$, $\mbox{{Ha}}=1300$

The profiles of $\langle u_{x}\rangle$ and $\langle u_{y}\rangle$ in figure 8c
also show a reversed flow in the negative $y$ part of the duct. It is nearly
as strong as the jet itself. The global structure of the asymmetric flow
structure can be characterized as that of a macrovortex.

The velocity waveforms in figure 8d demonstrate significant fluctuations in
the entire measurement region. The typical amplitude and frequency of the
fluctuations are higher in the zone of the jet at positive $y$ than in the
zone of reversed flow. It can also be seen that high-frequency turbulence
oscillations present in the flow at $\mbox{{Ha}}=200$ (see figure 8b) are not
seen at $\mbox{{Ha}}=1300$. This is attributed to quasi-two-dimensionality of
flow structures at such high Ha.

The profiles of mean velocities $\langle u_{x}\rangle$, $\langle u_{y}\rangle$
and standard deviations $\sigma_{x}$, $\sigma_{y}$ (see (6))
nondimensionalized by the mean duct velocity $U_{0}$ are shown in figure 10.
They further illustrate existence of the two types of flow structure found in
the experiments: a central jet and a macrovortex. The structure of the
macrovortex becomes more complex at high Re as illustrated by multiple extrema
of the profile of $\langle u_{x}\rangle$ in figure 10d.

(a) |  |  (b)
---|---|---

(c) |  |  (d)
---|---|---

Figure 10: Time-averaged characteristics of flows at (a) $\mbox{{Ha}}=200$,
$\mbox{{Re}}=800$; (b) $\mbox{{Ha}}=1300$, $\mbox{{Re}}=800$; (c)
$\mbox{{Ha}}=550$, $\mbox{{Re}}=5000$; and (d) $\mbox{{Ha}}=1300$,
$\mbox{{Re}}=3000$. Profiles measured at $x=1.6a$, $z=0$ are shown. 1 -
$\langle u_{x}\rangle/U_{0}$; 2 - $\langle u_{y}\rangle/U_{0}$; 3 -
$\sigma_{x}/U_{0}$; 4 - $\sigma_{y}/U_{0}$.

(a) |  |  (b)
---|---|---

(c) |  |  (d)
---|---|---

Figure 11: Velocity measurements at the center of the duct $y=z=0$ in the
cross-section $x=1.6a$. Time-averaged velocity components (a,c) and standard
deviations (b,d) based on signals of length 600 s are shown. $\mbox{{Re}}=800$
(a,c) and $\mbox{{Re}}=3000$ (b,d). 1– $\langle u_{x}\rangle$ and
$\sigma_{x}$. 2 – $\langle u_{y}\rangle$ and $\sigma_{y}$. The horizontal
dash-dot line in (a) and (b) shows the value equal to a half of the mean inlet
velocity of the jet. The vertical solid lines indicate $\mbox{{Ha}}_{cr}$
separating central jet and macrovortex regimes (see text).

It is possible to differentiate between flow states with a central jet and a
macrovortex using magnitudes of the time-averaged velocity components and
their standard deviations at the central point $y=z=0$. Results of such an
analysis carried out for two values of Re and the entire explored range of Ha
are presented in figure 11. To indicate the magnitude of velocity with respect
to the mean inlet velocity of the jet $U_{p}=110U_{0}$, a horizontal line
corresponding to $U_{p}/2$ is drawn in figure 11a,b.

The data for $\langle u_{x}\rangle$, $\sigma_{x}$, and $\sigma_{y}$ (although
not for $\langle u_{y}\rangle$, which is nearly zero for all Ha) clearly show
existence of two intervals of Ha with distinct behaviours. The critical
Hartmann numbers separating the intervals (shown as vertical lines in figure
11) are approximately $\mbox{{Ha}}_{cr}=380$ at $\mbox{{Re}}=600$ and
$\mbox{{Ha}}_{cr}=780$ at $\mbox{{Re}}=3000$. At
$\mbox{{Ha}}<\mbox{{Ha}}_{cr}$, $\langle u_{x}\rangle$ is large and gradually
decreases from its maximum $\langle u_{x}\rangle\approx 25U_{0}$ at
$\mbox{{Ha}}=200$ (as explained in section 3, this is the lowest Ha at which
accurate measurements are possible) to nearly zero. The intensity of velocity
fluctuations quantified by $\sigma_{x}$ and $\sigma_{y}$ also decreases and
reaches its minimum at $\mbox{{Ha}}=\mbox{{Ha}}_{cr}$.

The interval $\mbox{{Ha}}>\mbox{{Ha}}_{cr}$ is characterized by small values
of $\langle u_{x}\rangle/U_{0}$ and by $\sigma_{x}$ and $\sigma_{y}$ growing
with Ha to large values (about $12U_{0}$ to $14U_{0}$). Some decrease of the
standard deviation at very high Ha is observed at $\mbox{{Re}}=600$, but not,
in the explored range of Ha, at $\mbox{{Re}}=1300$.

A detailed discussion of the flow transformation causing the radical change of
behavior is provided in section 5. Here we only state the evident conclusion
that the flow takes the form of a central jet at
$\mbox{{Ha}}<\mbox{{Ha}}_{cr}$ and of a macrovortex at
$\mbox{{Ha}}>\mbox{{Ha}}_{cr}$.

The impact of the magnetic field on the properties of velocity signal measured
at the duct’s center is illustrated in figure 12. Properties of $u_{x}$ at
$\mbox{{Re}}=800$ and $\mbox{{Ha}}=200$ (central jet), $\mbox{{Ha}}=1000$
(macrovortex), and $\mbox{{Ha}}=350$ (close to $\mbox{{Ha}}_{cr}$) are shown.
We see the typical features of flow transformation observed in other MHD flows
in experiments and numerical simulations, e.g., by Kolesnikov and Tsinober
(1972); Alemany et al. (1979); Burattini et al. (2010); Verma (2017); Zikanov
et al. (2019). The transition into a quasi-two-dimensional form at growing Ha
is associated with suppression of high-frequency fluctuations, growth of
energy of low-frequency fluctuations, and steeper energy spectra with the
slope approaching $\sim f^{-3}$ (see figures 12a,b). A feature specific to the
jet flow and already visible in figure 11c is illustrated in figure 12c. The
typical amplitude of velocity oscillations varies non-uniformly with Ha. It
decrease from about $10U_{0}$ to about $5U_{0}$ when Ha changes from 200 to
350 and increases back to about $10U_{0}$ when Ha increases to 1000.

(a) |  |  (b)
---|---|---

|  (c) |
---|---|---

Figure 12: Properties of the signal of longitudinal velocity $u_{x}$ measured
at the center $y=z=0$ of the duct in the cross-section $x=1.6a$ at
$\mbox{{Re}}=800$ and $\mbox{{Ha}}=200$, 350, and 1000. (a), signals in a 15 s
window. (b), energy power spectra based on 600 s of observations. The apparent
peaks at $f=50$ Hz are caused by interference from the power supplies. The
slopes $\sim f^{-5/3}$ and $\sim f^{-3}$ are shown for guidance. (c), curves
of the standard deviation calculated in a moving 15 s window.

(a) (b)

Figure 13:  Nondimensionalized standard deviation $\sigma_{x}$ (a) and mean
value $\langle u_{x}\rangle$ (b) of longitudinal velocity measured at the duct
center $y=z=0$ versus the Stuart number N. Data for the cross-sections
$x=1.6a$ (I and II) and $x=5.2a$ (II-a) are shown. Measurements series I and
II are completed using, respectively, sensors of type 1 and 2 shown in figure
4a,b. Measurements of series I are completed using sensor of type 3 similar to
the sensor of type 2 in its design (see section 3.

The results of all the reported measurements of longitudinal velocity at the
duct center are summarized in figure 13. Standard deviation and mean value of
the longitudinal velocity scaled with $U_{0}$ are shown as functions of N. We
see that the two trends presented by figure 11c are universal. The mean
centerline velocity decreases with the strength of the magnetic field effect
and becomes small when the field is strong. The amplitude of the fluctuations
has two peaks, one at small and one at large N, and a minimum at intermediate
values of N.

The line $\langle u_{x}\rangle\sim\mbox{{N}}^{-0.5}$ in figure 13b corresponds
to the theoretically determined MHD decay of velocity in a submerged uniform
jet in an infinite domain (Davidson, 1995). It is not surprising that the
experimental data do not reproduce this slope closely. The real jet is
spatially evolving rather than uniform and occurs in a bounded domain of the
duct.

We see in figure 13a that all the data collected in the cross-section $x=1.6a$
at various values of Re and using sensors of two different types collapse into
one curve. The amplitude of the fluctuations is determined by N as a single
parameter. The situation is more complex for the mean velocity. One may argue
that the decay of $\langle u_{x}\rangle$ at $\mbox{{N}}<170$ is still largely
determined by the value of N. The small values of $\langle u_{x}\rangle$ at
$\mbox{{N}}>170$ do not follow any clear trend. This can be related to
insufficient averaging time in the experiments. The signal in figure 9
indicates flow dynamics that may require much longer than 600 s to measure
true average values. Another possible explanation is the effect of Hartmann
boundary layers, which are not uniquely identified by N.

Four distinct zones characterized by different behaviours of the fluctuation
amplitude and mean velocity can be identified in figure 13. In zone 1,
approximately at $\mbox{{N}}<50$, the fluctuation amplitude increases and the
mean velocity decreases with N. The normalized standard deviation reaches the
peak of $\sigma_{x}\approx 15U_{0}$ at $\mbox{{N}}=50$. Zone 2 at
$50<\mbox{{N}}<170$ is characterized by decrease of $\sigma_{x}$ with N at the
rate about $\sigma_{x}\sim\mbox{{N}}^{-0.5}$. The mean velocity decreases with
the rate higher than in zone 1. The smallest values of the fluctuation
amplitude are observed at $\mbox{{N}}=170$. In zone 3 defined approximately as
$170<\mbox{{N}}<1200$, the fluctuation amplitude grows again, initially as
$\sigma_{x}\sim\mbox{{N}}^{0.5}$ and then at a slower rate. We can also
identify zone 4 at $\mbox{{N}}>1200$, in which the amplitude decreases. The
mean centerline velocity remains low in zones 3 and 4.

It should be stressed that the just described universal behaviour is found in
the cross-section $x=1.6a$, but not in the cross-section $x=5.2a$ further
downstream. We hypothesize that the effect of Re on flow decay becomes
significant at that location, so the flow state cannot be uniquely determined
by N alone.

##  5 Discussion and concluding remarks

In this section we attempt an explanation of the experimental results.
Concluding remarks and suggestions for future work are provided at the end.

One important and evident conclusion of the experiments is that the classical
picture of the transformation of a round jet in a transverse magnetic field
illustrated in figures 2 and 3 is incomplete. The theoretical arguments, on
which the picture is based, fail to predict emergence of high-amplitude
velocity fluctuations and modification of the mean velocity profile found at
high Stuart numbers.

Another outcome of the experiment, namely that some flow properties are
accurately described by a function of N as a single parameter is not
surprising. The jet instability in the presence of a strong magnetic field is
anticipated to be mainly determined by the ratio between the Lorentz and
inertial forces, which is estimated by N.

The key mechanisms determining the flow evolution are the Kelvin-Helmholtz
instability of the jet, transition of the flow into quasi-two-dimensional
state, the Joule damping of velocity fluctuations, and the Hartmann friction.
Viscous dissipation is less significant in flows with such high Ha and Re. It
is ignored in the following discussion.

The jet instability is strong because of high Re and additional perturbations
introduced near the inlet by the Lorentz force in the form of transverse
velocity components associated with the transformation into anisotropic form
(in that way, the instability can be enhanced rather than suppressed by the
magnetic field). Presence of high-amplitude fluctuations at the distance of 44
mm, i.e. about seven jet diameters from the inlet is, therefore, not
surprising.

The Joule dissipation of velocity fluctuations is strong if the fluctuations
are significantly three-dimensional, but vanishes when the flow approaches
quasi-two-dimensional form. The damping effect due to friction in Hartmann
layers appears in quasi-two-dimensional flows, but it is, generally, weaker.
These arguments have been demonstrated earlier to be a valid explanation of
the so-called residual fluctuations (see, e.g., Kolesnikov and Tsinober, 1972;
Sukoriansky et al., 1986; Zikanov et al., 2019) – high-amplitude low-frequency
oscillations of velocity that appear in shear flows at high values of N. One
may take a more general view and consider this phenomenon as a manifestation
of the general effect, found prominently in flows with thermal convection
(Zikanov et al., 2021), of emergence of high-amplitude unsteady structures
when conventional turbulence is suppressed by the magnetic field and the flow
becomes quasi-two-dimensional.

We will now attempt to reconstruct, based on incomplete data, the flow
transformation represented by the results shown in figures 7–13. A sketch
illustrating the discussion is provided in figure 14.

Figure 14: Sketch of the typical flow structures at various values of Stuart
number.

In zones 1 and 2, approximately at $\mbox{{N}}<170$, the mean flow has a
central jet (see figure 14b). The centerline velocity decreases with
increasing N. The typical amplitude of centerline velocity fluctuations
increases with N in zone 1 and decreases with N in zone 2. The flow is
anisotropic, but still not quasi-two-dimensional (see figure 7a). We
hypothesize that the dominant effect of the magnetic field in zone 1 is the
effect consistent with the both observed trends, namely destabilization of the
jet by the Lorentz force. Similarly consistent dominant effect in zone 2 is
the suppression of three-dimensional velocity fluctuations by the Joule
dissipation. In other words, zone 1 corresponds to the situation when the
unstable three-dimensional perturbations are still growing in the measurement
cross-section. Zone 2 corresponds to the situation when the magnetic field is
stronger, and the growth has already been replaced by dissipation-caused decay
by the time fluid arrived at the measurement site.

The sharp minimum of the fluctuation amplitude at the boundary
$\mbox{{N}}=170$ between zones 2 and 3 indicates that the flow undergoes a
radical transformation. We hypothesize that zones 3 and 4 correspond to quasi-
two-dimensional flow states in the measurement cross-section. The central jet
is destroyed and replaced by a quasi-two-dimensional macrovortex (see figure
14c). The mean velocity at the centerline is nearly zero. The amplitude of
velocity fluctuations increases with N in zone 3 and decreases with N in zone
4. It must be stressed that the nature of velocity fluctuations is completely
different from the nature of fluctuations in zones 1 and 2. As illustrated in
figure 12, they are characterized by a steep slope of the energy power
spectrum and dominance of low-frequency oscillations. The only plausible
explanation of emergence of such fluctuations at high N is that they are
caused by large-scale unsteady quasi-two dimensional structures associated
with shear-flow instability of a macrovortex. The increase of the fluctuation
amplitude with N in zone 3 can be explained by a combination of two effects
promoting growth of such structures: stronger quasi-two-dimensionality and the
shift of the jet-to-vortex transition upstream. The decrease of the
fluctuation amplitude in zone 4 is attributed to the effect of growing
friction in the Hartmann boundary layers.

The experimental results are presented in this paper for a broad range of
parameters. This allows us to analyze the effect of the magnetic field and
identify distinct modes of flow’s behavior. At the same time, the experimental
data are inevitably limited in the sense that they do not produce a complete
picture of the flow. Many aspects of the explanation proposed in this section
must, therefore, be viewed as hypothetical until a confirmation in the form of
numerical simulation data is available. Such simulations are a necessary and
logical next step of the analysis. Their goals are (i) to produce detailed
descriptions of various states of the flow and (ii) to provide data for a
quantitative analysis based on detailed comparison between the experimental
and numerical data.

Funding. Financial support was provided by the award 20-69-46067 from the
Russian Science Foundation.

Declaration of interests. The authors report no conflict of interest.

## References

  * Bühler [2007] L. Bühler.  Liquid metal magnetohydrodynamics for fusion blankets.  In _Magnetohydrodynamics_ , pages 171–194. Springer, 2007.
  * Abdou et al. [2015] M. Abdou, N. B. Morley, S. Smolentsev, A. Ying, S. Malang, A. Rowcliffe, and M. Ulrickson.  Blanket/first wall challenges and required R&D on the pathway to DEMO.  _Fusion Eng. Des._ , 100:2 – 43, 2015.
  * Cukierski and Thomas [2008] K. Cukierski and B. G. Thomas.  Flow control with local electromagnetic braking in continuous casting of steel slabs.  _Metallurgical and Materials Transactions B_ , 39(1):94–107, 2008.
  * Michalke [1984] A. Michalke.  Survey on jet instability theory.  _Progress in Aerospace Sciences_ , 21:159–199, 1984.
  * Thomas and Goldschmidt [1986] F.O. Thomas and V.W. Goldschmidt.  Structural characteristics of a developing turbulent planar jet.  _Journal of Fluid Mechanics_ , 163:227–256, 1986.
  * Stanley et al. [2002] S.A. Stanley, S. Sarkar, and J.P. Mellado.  A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation.  _Journal of Fluid Mechanics_ , 450:377–407, 2002.
  * Sato [1960] H. Sato.  The stability and transition of a two-dimensional jet.  _Journal of Fluid Mechanics_ , 7(1):53–80, 1960\.
  * Sato and Sakao [1964] H. Sato and F. Sakao.  An experimental investigation of the instability of a two-dimensional jet at low reynolds numbers.  _Journal of Fluid Mechanics_ , 20(2):337–352, 1964.
  * Kalis and Kolesnikov [1984] Kh. E. Kalis and Yu. B. Kolesnikov.  Stability of a rotating transverse axisymmetric flow in an axial magnetic field.  _Magnetohydrodynamics_ , 19(3):291–294, 1984\.
  * Klyukin and Thess [1993] A. A. Klyukin and A. Thess.  The investigation of the turbulence origination process in an electrically driven MHD flow.  _Magnetohydrodynamics_ , 29(4):341–344, 1993\.
  * Kljukin and Kolesnikov [1989] A.A. Kljukin and Yu.B. Kolesnikov.  MHD instabilities and turbulence in liquid metal metal shear flows.  In J. Lielpeteris and R. Moreau, editors, _Liquid Metal Magnetohydrodynamics_ , pages 449–454. Kluwer Academic Publishers, 1989.
  * Moffatt and Toomre [1967] H.K. Moffatt and J. Toomre.  The annihilation of a two-dimensional jet by a transverse magnetic field.  _Journal of Fluid Mechanics_ , 30(1):65–82, 1967\.
  * Davidson [1995] P.A. Davidson.  Magnetic damping of jets and vortices.  _Journal of Fluid Mechanics_ , 299:153–186, 1995.
  * Müller and Bühler [2001] U. Müller and L. Bühler.  _Magnetofluiddynamics in channels and containers_.  Springer, Berlin, Heidelberg, 2001.
  * Ludford [1960] G.S.S. Ludford.  Inviscid flow past a body at low magnetic reynolds number.  _Reviews of Modern Physics_ , 32(4):1000, 1960\.
  * Hasimoto [1960] H. Hasimoto.  Steady longitudinal motion of a cylinder in a conducting fluid.  _Journal of Fluid Mechanics_ , 8(1):61–81, 1960\.
  * Molokov [1993] S. Molokov.  Single-component magnetohydrodynamic flows in a strong uniform magnetic field. 2. rotation of an axisymmetric body.  _Magnetohydrodynamics_ , 29(2):175–180, 1993\.
  * Moffatt [1964] H.K. Moffatt.  Electrically driven steady flows in magnetohydrodynamics.  pages 945–953, 1964.
  * Kim and Choi [2004] J. Kim and H. Choi.  Large eddy simulation of magnetic damping of jet.  In _Int. Congress of Th. And Appl. Mech_ , 2004.
  * Krasnov et al. [2021] D. Krasnov, Ya. Listratov, Yu. Kolesnikov, I. Belyaev, N. Pyatnitskaya, E. Sviridov, and O. Zikanov.  Transformation of a submerged flat jet under strong transverse magnetic field.  _EPL (Europhysics Letters)_ , 134(2):24003, 2021\.
  * Belyaev et al. [2021] I.A. Belyaev, N. Yu. Pyatnitskaya, N.A. Luchinkin, D Krasnov, Yu B Kolesnikov, Ya. I. Listratov, I.S. Mironov, O. Zikanov, and E.V. Sviridov.  Flat liquid metal jet affected by a transverse magnetic field.  _Magnetohydrodynamics_ , 57(2):211–222, 2021\.
  * Krasnov et al. [2022] D. Krasnov, Yu. Kolesnikov, Ya. Listratov, O. Zikanov, S. Molokov, and I. Belyaev.  Direct numerical simulations of a submerged liquid metal jet in a transverse magnetic field.  _Journal of Fluid Mechanics_ , 2022.  In preparation.
  * Belyaev et al. [2017] I. A. Belyaev, V. G. Sviridov, V. M. Batenin, D. A. Biryukov, I. S. Nikitina, S. P. Manchkha, N. Yu. Pyatnitskaya, N. G. Razuvanov, and E. V. Sviridov.  Test facility for investigation of heat transfer of promising coolants for the nuclear power industry.  _Thermal Engineering_ , 64(11):841–848, Nov 2017\.  ISSN 1555-6301.  doi: 10.1134/S0040601517110027.  URL https://doi.org/10.1134/S0040601517110027.
  * Kirillov [2007] P.L. Kirillov.  _Thermophysical properties of materials for nuclear engineering_.  IzdAT, 2007.
  * Mistrangelo and Bühler [2010] C. Mistrangelo and L. Bühler.  Perturbing effects of electric potential probes on MHD duct flows.  _Experiments in fluids_ , 48(1):157–165, 2010\.
  * Murgatroyd [1953] W. Murgatroyd.  Cxlii. experiments on magneto-hydrodynamic channel flow.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 44(359):1348–1354, 1953.
  * Branover [1978] H. Branover.  Magnetohydrodynamic flow in ducts.  _New York_ , 1978.
  * Zikanov et al. [2014] O. Zikanov, D. Krasnov, T. Boeck, A. Thess, and M. Rossi.  Laminar-turbulent transition in magnetohydrodynamic duct, pipe, and channel flows.  _Applied Mechanics Reviews_ , 66(3), 2014.
  * Kolesnikov and Tsinober [1972] Yu.B. Kolesnikov and A.B. Tsinober.  Two-dimensional turbulent flow behind a circular cylinder.  _Magnetohydrodynamics_ , 8(3):300–307, 1972.
  * Alemany et al. [1979] A. Alemany, R. Moreau, P. L. Sulem, and U. Frisch.  Influence of an external magnetic field on homogeneous MHD turbulence.  _Journal of Fluid Mechanics_ , 18:277–313, 1979.
  * Burattini et al. [2010] P. Burattini, O. Zikanov, and B. Knaepen.  Decay of magnetohydrodynamic turbulence at low magnetic reynolds number.  _Journal of Fluid Mechanics_ , 657:502–538, 2010.
  * Verma [2017] M. K. Verma.  Anisotropy in quasi-static magnetohydrodynamic turbulence.  _Reports on Progress in Physics_ , 80(8):087001, 2017.
  * Zikanov et al. [2019] O. Zikanov, D. Krasnov, T. Boeck, and S. Sukoriansky.  Decay of turbulence in a liquid metal duct flow with transverse magnetic field.  _Journal of Fluid Mechanics_ , 867:661–690, 2019.
  * Sukoriansky et al. [1986] S. Sukoriansky, I. Zilberman, and H. Branover.  Experimental studies of turbulence in mercury flows with transverse magnetic fields.  _Experiments in fluids_ , 4(1):11–16, 1986.
  * Zikanov et al. [2021] O. Zikanov, I. Belyaev, Y. Listratov, P. Frick, N. Razuvanov, and V. Sviridov.  Mixed convection in pipe and duct flows with strong magnetic fields.  _Applied Mechanics Reviews_ , 73(1):010801, 2021\.