Ngoc Minh

Trang Vu

Laurent Lef Èvre
email: laurent.lefevre@lcis.grenoble-inp.fr

Bernhard Maschke
email: maschke@lagep.univ-lyon1.fr

A structured control model for the thermo-magneto-hydrodynamics of plasmas in tokamaks

Introduction

A tokamak is a facility constructed with the shape of a torus (or doughnut) in which a plasma is magnetically confined and heated in order to produce nuclear fusion reactions (see Fig. 1 for a schematic view and the classical Wesson's monograph [START_REF] Wesson | Tokamaks[END_REF] for a large comprehensive reference textbook). The ITER (International Thermonuclear Experimental Reactor ) project is currently building the world's largest experimental tokamak nuclear fusion reactor. However many challenges remain to prove the scientific feasibility of this goal and then to move towards a fully functional plant. Some of these are control problems such as [START_REF] Pironti | Fusion, tokamaks and plasma control[END_REF][START_REF] Walker | Emerging applications in tokamak plasma control[END_REF][START_REF] Ariola | Magnetic Control of Tokamak Plasmas[END_REF] the vertical stabilization of the plasma, controlling the overall shape of the plasma, controlling the magneto-hydrodynamic (MHD) instabilities or controlling the current, temperature and pressure density profiles.

There exist many kinds of models for the plasma dynamics in tokamaks including the 1D resistive diffusion equation for the poloidal flux [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF], Hamiltonian formulation for the microscopic kinetic equations [START_REF] Morrison | Noncanonical Hamiltonian density formulation of hydrodynamics and ideal magnetohydrodynamics[END_REF][START_REF] Marsden | The Hamiltonian structure of the Maxwell-Vlasov equations[END_REF] or macroscopic fluid-like description derived from these microscopic equations [START_REF] Tassi | Hamiltonian formulation and analysis of a collisionless fluid reconnection model[END_REF].

Up to now, the so-called 1D resistive diffusion model for the poloidal flux (see [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF]Chap.6] for model derivation and details) has been extensively used as a "control model" for many control designs related to the plasma total current or current profile in the tokamak. In this model the plasma is assumed to have reached the hydrodynamic equilibrium profile and displacement currents are 1: Schematic view of a tokamak with the electrical solenoids: the magnetic field generated by the three magnets makes the plasma gas ions following helicoidal trajectories along the torus neglected in the Maxwell equations. The model then reduces to a 1D diffusionlike parabolic equation for the radial profile of the plasma poloidal flux (or current density). Readers could refer to [START_REF] Witrant | A control-oriented model of the current profile on tokamak plasma[END_REF] for investigations on this model for control purposes or to [START_REF] Ouarit | Validation of plasma current profile model predictive control in tokamaks via simulations[END_REF] for application to model-based predictive control. A similar model has been used to solve the current profile optimal tracking problem [START_REF] Ou | Optimal tracking control of current profile in tokamaks[END_REF] or to design robust controller for the poloidal magnetic flux profile [START_REF] Ou | Robust control design for the poloidal magnetic flux profile evolution in the presence of model uncertainties[END_REF]. Feedback control using Lyapunov approach in [START_REF] Argomedo | A strict control lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF], or sliding mode in [START_REF] Gaye | Sliding mode stabilization of the current profile in tokamak plasmas[END_REF] are also proposed.

However, this 1D resistive diffusion model makes very restrictive assumptions on the plasma resistivity (which in fact depends strongly and non-linearly of the temperature) and neglects some nonlinear MHD couplings which produces among others the bootstrap current2 . The next objective for plasma high confinement requires now a better understanding and exploitation (for control purposes) of TMHD interdomain couplings. Moreover, a study of the complete TMHD model, including mass and entropy balance equations is necessary to consider the burn control problem (i.e. control of the plasma with the fusion reaction) which comes next after plasma confinement problems (see for instance preliminary works of [START_REF] Schuster | Nonlinear lyapunov-based burn control in fusion reactor[END_REF][START_REF] Schuster | Burn control in fusion reactors via nonlinear stabilization techniques[END_REF][START_REF] Schuster | Control of a nonlinear pde system arising from non-burning tokamak plasma transport dynamics[END_REF][START_REF] Boyer | Adaptive nonlinear burn control in tokamak fusion reactors[END_REF][START_REF] Boyer | Nonlinear burn control in tokamak fusion reactors via output feedback[END_REF]).

On the other hand, both kinetic and fluid models, which include MHD and thermal couplings, have already been considered in the literature. Most of these models already possess a Hamiltonian structure which is considered essential by plasma physicists since the fundamental laws governing charged particle dynamics are Hamiltonian. Therefore the preservation of the Hamiltonian structure provides some confidence that the truncations that are used to derive the fluid model have not introduced unphysical phenomena. The presence of the Hamiltonian structure has the additional benefit of providing important tools for calculations such as the MHD energy principle, solvability conditions for the equilibrium equations, Casimir invariants, etc. Readers could refer to [START_REF] Morrison | Hamiltonian and action principle formulations of plasma physics[END_REF] for a general argumentation about Hamiltonian models for plasma physics, to [START_REF] Morrison | Noncanonical Hamiltonian density formulation of hydrodynamics and ideal magnetohydrodynamics[END_REF][START_REF] Marsden | The Hamiltonian structure of the Maxwell-Vlasov equations[END_REF] for early work on Hamiltonian formulation of 3D kinetic tokamak model or to [START_REF] Tassi | Hamiltonian formulation and analysis of a collisionless fluid reconnection model[END_REF] for recent developments on the Hamiltonian formulation of a MHD fluid model. Although previous models using Hamiltonian formalisms already exist, to the best of our knowledge none of them make use of Dirac structure, nor port-variables, and they only apply to closed systems. Therefore they are not convenient for control or observation purposes. The purpose of this paper is to propose a structured 3D 3 port-Hamiltonian model for the plasma dynamics in tokamaks which is suitable for control (eventually after some symplectic reduction or discretization steps) but still captures most of the physical properties of the real plant. More precisely, we aim at representing explicitly the balance equations in the electromagnetic (EM) and "material" domains (see figure 2), their (multi-domains) closure equations and interconnection (Stokes-Dirac) structures, as well as all the possible distributed or boundary actuations (i.e. boundary energy flows). In the tokamak example, distributed and boundary actuations may occur both in the electromagnetic or "material" domains (e.g. through heat or mass injection). To build this port-Hamiltonian model, we will close the Maxwell field equations written in covariant form (see for instance [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF]) with material balance equations and closure equations (including multi-domain couplings) which will be derived from the kinetic theory of gases and the Boltzmann equation (see [START_REF] Braginskii | Reviews of Plasma Physics[END_REF]). This will require also the derivation of the irreversible entropy source term (from the Gibbs-Duhem relation, following the "port-based" approach in [START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF]Chap.3,p.154] or [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]) which is contributed by all other domains (EM, mechanic, thermal and hydraulic).

This paper is organized as follows. Section 2 recalls the port-Hamiltonian formulation for distributed parameters systems with boundary energy flows. In section 3, the electromagnetic field equations are written using this port-Hamiltonian formulation. In section 4 the material domain balance equations for mass, momentum and energy are derived from the Boltzmann equation (kinetic 3 The chosen approach needs to identify the velocity vector field and the momentum 1-form through the canonical Hodge star operator (see section 2.1.4 hereafter). Therefore the model is expressed in a three dimensional euclidean space. Note however that it may be reduced to a 2D or a 1D model by using either the axisymmetry or the quasi-static equilibrium assumptions and curvilinear coordinates [START_REF] Vu | Structure preserving reduction for thermo-magneto plasma control model[END_REF][START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF]. Another possible approach is to consider the velocity directly as a 1-form such as in [START_REF] Polner | A hamiltonian vorticity-dilatation formulation of the compressible euler equations[END_REF] for the compressible Euler equations theory). The irreversible entropy production is also deduced from these three balance equations by using the Gibbs-Duhem equation. Section 5 is dedicated to closure equations which are the magneto-hydrodynamic coupling (through Lorentz forces); canonical interdomain couplings for the mechanical, hydraulical and thermal domains and a multi-domain resistivity field (R -field) which accounts for the entropy irreversible production. Finally, the obtained set of balance and closure equations for the plasma dynamics in a tokamak is summarized in section 6. A graphical (Bond Graph) representation of the model is proposed to underline its internal multidomain structure. The paper ends with a brief conclusion and perspectives for future works.

Port-based modelling for distributed parameters systems of conservation laws with boundary energy flows

In this section, the port-Hamiltonian formulation for systems of conservation laws is first recalled. Port-based modelling techniques and languages have been extensively used these last decades to model, simulate and control a wide variety of lumped parameters physical systems [START_REF] Karnopp | System Dynamics: Modeling and Simulation of Mechatronic Systems[END_REF][START_REF] Van Der Schaft | L 2 -Gain and Passivity Techniques in Nonlinear Control[END_REF][START_REF] Stramigioli | Modeling and IPC Control of Interactive Mechanical Systems: a Coordinate-free Approach[END_REF][START_REF] Kugi | Non-linear Control Based on Physical Models: Electrical, Hydraulic and Mechanical Systems[END_REF][START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF]Chap.3]. More recently, an intrinsic formulation of port-Hamiltonian models for distributed parameter systems (described by partial differential state space equations) with boundary energy flow have been proposed [START_REF] Maschke | Canonical interdomain coupling in distributed parameter systems: an extension of the symplectic gyrator[END_REF][START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF].

It is known that Stokes-Dirac structures allow to extend the Hamiltonian description of distributed-parameters systems to include variable boundary conditions, leading to open distributed-parameter port-Hamiltonian systems with boundary ports [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF].

This modelling approach has been applied successfully to many hyperbolic systems as varied as transmission lines models [START_REF] Golo | Hamiltonian discretization of boundary control systems[END_REF], beam equations [START_REF] Macchelli | Modeling and control of the Timoshenko beam. the Distributed Port Hamiltonian approach[END_REF] or shallow water equations [START_REF] Hamroun | Port-based modelling for open channel irrigation systems[END_REF][START_REF] Hamroun | A port-controlled hamiltonian approach to geometric reduction of distributed parameters systems -application to the shallow water equations[END_REF]. The same canonical Stokes-Dirac structure may be used to build port-Hamiltonian formulations of parabolic models issued from (linear) irreversible thermodynamics such as in the case of heat and mass transport phenomena in adsorption columns [START_REF] Baaiu | Port-based modelling of mass transfer phenomena[END_REF][START_REF] Baaiu | Structure-preserving infinite dimensional model reduction: Application to adsorption processes[END_REF] or ionic diffusion in fuel cell [START_REF] Franco | Multi-scale bond graph model of the electrochemical dynamics in a fuel cell[END_REF] or Ionic Polymer-Metal Composites [START_REF] Nishida | Multi-scale distributed parameter modeling of ionic polymer-metal composite soft actuator[END_REF][START_REF] Nishida | Multiscale distributed porthamiltonian representation of ionic polymer-metal composite (ipmc)[END_REF]. A more surprising parabolic plasma resistive diffusion model in the next section will show the interest of the approach in multi-physics (e.g. Thermo-MagnetoHydroDynamics) models and more involved geometries (e.g. toroidal).

Co-variant form and exterior calculus notations

Let us first recall some exterior calculus definitions and notations used in the sequel [16] [50, chap.4].

Differential forms

Definition 2.1. On an n-dimensional spatial domain Ω ξ, a k-form (k ≤ n),
is the field ω k ξ of k-alternated forms on T ξ Ω, that is k-linear (i.e. linear with respect to each of its arguments) forms on the tangent vectors: ω k ξ : (T ξ Ω) k → R that satisfy for any permutation π in {1, .., k} and for any k-tuple (ξ 1 , . . . , ξ k ) ∈ (T ξ Ω)

k : ω k ξ ξ π(1) , . . . , ξ π(k) = σ (π) ω k ξ (ξ 1 , . . . , ξ k ) where σ (π) denotes the signature of the permutation π. The set of k-forms is denoted by Λ k (Ω) .

Differential k-forms are endowed with a product, called exterior product, denoted by ∧, which we use to compute the power and define the passivity properties, and a derivation which is called exterior derivation, denoted by d. These operators are defined as follows.

Exterior product

Consider a k-form ω k ∈ Λ k (Ω) and a l-form ω l ∈ Λ l (Ω), their exterior product (also called wedge product or Grassman product) is the (k + l)-form ω k ∧ ω l :

∧ : Λ k (Ω) × Λ l (Ω) → Λ k+l (Ω)
such that for any (k + l)-tuple of vector fields (ξ 1 , . . . , ξ k+l ):

ω k ∧ ω l (ξ 1 , . . . , ξ k+l ) = π∈S σ (π) ω k ξ π(1) , . . . , ξ π(k) ω l ξ π(k+1) , . . . , ξ π(k+l) ∈ Λ k+l (Ω)
where S denotes the sets of shuffle permutation that is the permutation π satisfying π (1) < . . . < π (k) and π (k + 1) < . . . < π (k + l).

Exterior derivation

The exterior derivation (also called co-boundary map) of a k-form ω k ∈ Λ k (Ω), denoted dω k , is a derivation of degree 1 which maps Λ k (Ω) into Λ k+1 (Ω)

d : Λ k (Ω) → Λ k+1 (Ω)
and satisfies:

• linearity d(ω + η) = dω + dη • anti-derivation d ω k ∧ ω l = (-1) k dω k ∧ ω l + ω k ∧ dω l ω k ∈ Λ k (Ω) , ω l ∈ Λ l (Ω)
• if ω 0 is a smooth function, dω 0 is the differential of the function

• d ∧ d = 0
• the derivation is local: for any open set U ⊂ Ω , if the restrictions to U of two k-forms coincide then also the restrictions of their exterior derivatives.

Another operator associated with an inner product on the exterior differential forms is the Hodge star operator.

Hodge star operator

The Hodge star (cf. [11, V.A.,p.295]) is essential in the definition of the Hamiltonian (or energy) functions of many physical systems which admits "quadratic" energy. Assume that the vector space Λ k (Ω) of k-forms on Ω is endowed with an inner product denoted by ( , ), then the Hodge star is a linear map that transforms a k-form ω k into a (n -k)-form, also called the dual or the pseudo (n -k)-form ω k :

: Λ k (Ω) → Λ n-k (Ω) Thus: α, ω k = ˆΩ α ∧ ω k ∀α ∈ Ω k (Ω)
In this paper, the Hodge star operator definition makes use of the canonical inner product. It also depends on the orientation defined for Ω (the sign is changed when reversing orientation). In the sequel, this Hodge star operator will be used to define constitutive closure equations which usually depend on the considered material such as Lorentz forces, Joule dissipation, viscosity, etc. We will also define the macroscopic kinetic energy of the gas as one of the constitutive equations since it must be considered in the plasma example only as an average value for kinetic energies of several species of particles with specific velocities and masses. For other examples (such as anisotropic material domain), non trivial metric must be used to define properly the corresponding Hodge star operator.

The contraction of a k-form with some vector field is defined hereafter. It may be seen as an anti-derivation formulae.

Contraction

The contraction of a k-form ω k ∈ Λ k (Ω) by some vector field v ∈ T ξ Ω , denoted by i v ω k , is a derivation of degree -1 which maps:

i v : Λ k (Ω) → Λ k-1 (Ω)
such that for any (k -1)-tuple of vector fields (ξ 1 , . . . , ξ k-1 ):

i ξ ω k (ξ 1 , . . . , ξ k-1 ) = ω k (ξ, ξ 1 , . . . , ξ k-1 )
Finally we shall recall the Stokes' theorem which is fundamental for the definition of the Dirac structure underlying systems of conservation laws with energy flux at the boundary.

2.1.6. Stokes's theorem Theorem 2.2. Consider a spatial domain Ω ∈ R n being an k-dimensional smooth manifold with smooth (k -1)-dimensional boundary ∂Ω. Then for any (k -1)-form ω k-1 with compact support in R n , one has: ˆΩ dω k-1 = ˆ∂Ω ω k-1 (1)
Note that from here, the spatial domain Ω is considered to be a 3D domain Ω ⊂ R 3 .

Systems of conservation laws

We consider here the case of systems of two conservation laws in canonical interaction.

Definition 2.3. Consider the two conserved quantities as being two 2-forms: q ∈ Λ 2 (Ω) and p ∈ Λ 2 (Ω). Consider also the system of conservation laws, with flux variables β q ∈ Λ 1 (Ω) and β p ∈ Λ 1 (Ω) for each conserved quantity, defined by the Hamiltonian density function H : Λ 2 (Ω) × Λ 1 (Ω) → Λ 3 (Ω) resulting in the total Hamiltonian H := ´Ω H (q, p) ∈ R. The system of two canonically interacting conservation laws is then defined by:

∂ ∂t q p + d β q β p = 0 and β q β p = ε 0 1 1 0 δH δq δH δp (2)
where ε ∈ {-1, +1} depends on the fluxes sign convention on the physical domain.

This system of two conservation laws may be also written as follows:

∂ ∂t q p = ε 0 d d 0 δH δq δH δp (3)
that is as an infinite-dimensional Hamiltonian system defined with respect to the matrix differential operator :

J = ε 0 d d 0 (4)
and generated by the Hamiltonian function H [START_REF] Olver | Applications of Lie Groups to Differential Equations[END_REF].

In order to generate a Hamiltonian systems, the matrix differential operator J defined in (4) should satisfy the properties of a Hamiltonian operator, that is it should be skew-symmetric and satisfy the Jacobi identities. A short calculus shows that the skew-symmetry holds only for functions with compact support in the spatial domain Ω or satisfying Dirichlet or Neumann homogeneous boundary conditions.

Stokes-Dirac structures and systems of balance equations with boundary energy flows

Interested in observation and control problems, one must consider more general (dynamic) boundary conditions where some energy is exchanged through the boundary of the spatial domain. Therefore the matrix differential operator J must be extended to a Dirac structure, called Stokes-Dirac structure [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF][START_REF] Maschke | Advanced Topics in Control Systems Theory[END_REF][START_REF] Macchelli | Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach, chapter Infinite-dimensional Port-Hamiltonian Systems[END_REF]. Dirac structures [START_REF] Courant | Dirac manifolds[END_REF][START_REF] Dalsmo | On representations and integrability of mathematical structures in energy-conserving physical systems[END_REF] are a geometric perspective to skewsymmetric tensors, actually corresponding to their graph, which generalize the tensors associated with Poisson brackets or pre-symplectic forms. They have been introduced in classical mechanics to represent systems with constraints. Then they have been used to include input-output port variables in finite dimensional port-Hamiltonian models [START_REF] Van Der Schaft | The Hamiltonian formulation of energy conserving physical systems with external ports[END_REF].

Definition 2.4. Let F and E be linear spaces with a pairing e|f ∈ L linear space, e ∈ E and f ∈ F. By symmetrizing the pairing we obtain a symmetric bilinear form . on F × E defined as:

(f 1 , e 1 ) , (f 2 , e 2 ) := e 1 |f 2 + e 2 |f 1 (5)
A Dirac structure is a linear subspace D ⊂ F × E such that D = D ⊥ , with ⊥ denoting the orthogonal complement with respect to the bilinear form . .

Making use of the Stokes theorem, this definition of Dirac structure has been used to represent canonical systems of conservation laws with boundary energy flows as described in the proposition 2.5 hereafter, defining a so-called infinite dimensional Stokes-Dirac interconnection structure [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF].

Proposition 2.5. Consider the product spaces of k-forms:

F = Λ 2 (Ω) × Λ 2 (Ω) × Λ 1 (∂Ω) (f p , f q , f b ) E = Λ 1 (Ω) × Λ 1 (Ω) × Λ 1 (∂Ω) (e p , e q , e b ) (6)
Consider the linear subspace D of the bond space B = F × E:

D = {(f p , f q , f b , e p , e q , e b ) ∈ F × E| f p f q = ε 0 d d 0 e p e q , f b e b = ε 0 0 -1 e p|∂Ω e q|∂Ω (7)
where ε ∈ {-1, +1} and | ∂Ω denotes restriction to the boundary

∂Ω = [0, L].
Then D is a Dirac structure with respect to the non degenerated bilinear form between F and E:

(e p , e q , e b ) | (f p , f q , f b ) = ˆΩ [e p ∧ f p + e q ∧ f q ] + e b , f b ∂Ω (8)
The interconnection structure D is called a Stokes-Dirac structure

As a consequence of proposition 2.5 one may define a Hamiltonian system with respect to this Stokes-Dirac structure as follows.

Definition 2.6. The boundary port-Hamiltonian system of two conservation laws with state space Λ 2 (Ω) × Λ 2 (Ω) (q, p) and boundary port variables spaces

Λ 1 (∂Ω) × Λ 1 (∂Ω) (f b , e b )
, is the Hamiltonian system defined with respect to the Stokes-Dirac structure D given in proposition 2.5 and generated by the Hamiltonian functional H (q, p), as follows:

- ∂p ∂t , - ∂q ∂t , f b , δH δp , δH δq , e b ∈ D
The choice of boundary conditions has obviously to be added to the definition of a boundary port-Hamiltonian system in order to define a Cauchy problem. In fact, in the linear case, a boundary port-Hamiltonian system defines a class of well-posed systems. For any solution, the isotropy condition of the Dirac structure implies the balance equation on the Hamiltonian:

dH dt = e b , f b ∂ (9)

Co-variant formulation for the dynamics in the electromagnetic domain

We detail hereafter the port-Hamiltonian co-variant formulation for electrodynamics [49] [50, chap.4]. The model is inspired from [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF][START_REF] Nishida | Port-based modeling of magnetohydrodynamics equations for tokamaks[END_REF] and adapted to include additional Lorentz forces. The main used equations of electromagnetism are summarized in the Table 1 below, both in vector calculus and covariant (differential forms) notations.

The electric and magnetic field intensities are naturally identified in [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF] with the 1-forms: E, H ∈ Λ 1 (Ω); the magnetic, electric flux and free current Name Vector calculus formulation Co-variant formulation Gauss's law ∇.D = ρ e dD = ρ e Gauss's law for magnetism densities are identified with the 2-forms: B, D, J ∈ Λ 2 (Ω); and the charge density is written as a 3-form: ρ e ∈ Λ 3 (Ω). Note that the plasma velocity v is defined as a vector field, which is needed in most of the contraction notations.

∇.B = 0 dB = 0 Maxwell-Faraday ∇ × E = -∂B ∂t dE = -∂B ∂t Maxwell-Ampere ∇ × H = J + ∂D ∂t dH = J + ∂D ∂t Constitutive equations D = E B = µH D = E B = µ H
The electrodynamic equations may be formulated, using the co-variant formulation [START_REF] Jamiolkowski | Classical Electrodynamics[END_REF], as a port-Hamiltonian system [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF]. The Hamiltonian functional is total electromagnetic energy H (D, B) in the domain Ω which is the integration of the energy density

H = 1 2 [E ∧ D + H ∧ B]: H (D, B) = ˆΩ H = 1 2 ˆΩ[ 1 ε D ∧ D + 1 µ B ∧ B] (10)
The electric and magnetic fields intensities are therefore intensive variables defined as variational derivatives of H (D, B), that is respectively the two 1-forms:

E, H ∈ Λ 1 (Ω) : δH δD = E (t, x) δH δB = H (t, x) (11)
From the power continuity and Stokes theorem, one may derive the following Stokes-Dirac structure (see [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF] for the definition and details) on the space of flows

F =Λ 2 (Ω) × Λ 2 (Ω) × Λ 2 (Ω) × Λ 1 (∂Ω) and efforts E = Λ 1 (Ω) × Λ 1 (Ω) × Λ 2 (Ω) × Λ 1 (∂Ω): D EM =            (f el , f mg , f d , f ∂Ω , e el , e mg , e d , e ∂Ω ) ∈ F × E | f el f mg = 0 -d d 0 e el e mg + 1 0 f d , e d = 1 0 e el e mg and f ∂Ω e ∂Ω = e el | ∂Ω e mg | ∂Ω            (12)
where ω 1 ∂Ω denotes the restriction of the 1-forms at the boundary of the domain. Then the electromagnetic field equations may be implicitly defined as a port-Hamiltonian systems:

- ∂D ∂t , - ∂B ∂t , f d , f ∂Ω , E, H, e d , e ∂Ω ∈ D EM (13)
with boundary port variables:

f ∂Ω = E| ∂Ω e ∂Ω = H| ∂Ω
whose wedge product forms the Poynting vector at the boundary of the spatial domain ∂Ω. The distributed port variables in the domain Ω are:

f d = J e d = E
where J = J Ω + J ni is the total current density including the Ohmic and the external non-inductive current (created by the RF waves). This distributed port is terminated with the Ohm's resistive equation η J Ω (with η the Hodge star related to the plasma resistivity) and connected to the material domain via the Lorentz force equation. The Lorentz force E L may be written using the contraction:

E L = -i v B ( 14
)
and we can derive the relation:

E = η J Ω -i v B (15)
The implicit formulation (13) reveals to be extremely useful for the systems and control theory [START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF]. Among others, it satisfies the following energy balance equation:

dH dt = -´∂Ω e ∂Ω ∧ f ∂Ω -´Ω e d ∧ f d
In other words, with the appropriate choice of input and output variables, the system passivity is easily proved and the energy evolution is governed by the external energy supplied through the boundary (and by the Ohmic dissipation rate).

A commonly used control model for plasma current density profiles in tokamaks is the so-called 1D resistive diffusion equation [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF]. This 1D model may be derived from the previous 3D model by using the usual axisymmetry and quasi-static assumptions (the Maxwell-Ampere and Maxwell-Faraday equations may be written in toric magnetic coordinates in the form of a system where all variables depends only on the magnetic surface radius ρ). In [START_REF] Vu | Geometric discretization for a plasma control model[END_REF], a symplectic spectral discretization is proposed to validate this 1D control model.

The influences of the dynamics in material domain on certain parameters in the electromagnetic domain, like the resistivity η and the bootstrap noninductive bootstrap current J bt , is ignored in this resistive diffusion equation and in the many control studies which make use of it. In the sequel, we determine first the material balance equations and then the coupling between the electromagnetic and material domains which allows to include the dependance of the resistivity (and of the diffusivity) with the temperature and induction field profiles. The bootstrap current still remains an unmodelled source term and its explicit formulation in the model (as an inter-domain coupling) is among the prospects of this work.

Balance and closure equations in the material domain

In this section, the material balance equations for mass, momentum, energy and entropy are derived from the Boltzmann equation (in kinetic theory). These balance equations are computed for a moving material domain. The connection between the classical macroscopic transport equation and the port-based formulation is made by using the material derivative in co-variant form and the Gibbs-Duhem relation for the internal energy (cf. [START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF]Chap. 3,p. 154]). The result is a port-based fluid-like model expressed in terms of pairs of power conjugated variables.

Kinetic theory and macroscopic transport equations

Using a kinetic description for the plasma dynamics, the material balance equations (mass, momentum, and energy balances) may be obtained from the integrals of the Boltzmann equation . These computations are usually performed in chosen coordinate systems for a 3D Euclidean space (see e.g. [10, p. 205]). Therefore, in this subsection, we recall the classical results obtained using the Nabla operator to define the grad, curl and div operators respectively as (∇, ∇×, ∇.), before stating the covariant formulation in subsection 4.2.

Boltzmann equation

The Boltzmann equation is a state equation which describes the time evolution of the distribution functions f a (t, r, v) for particles with position r and velocity v in the coordinate system β:

∂f a ∂t + ∂ ∂r β (v β f a ) + ∂ ∂v β F aβ m a f a = C a (16)
Here the subscript a stands for the different considered species, such as electrons a = e, ions a = i, or neutral particles a = n in our case. The external force:

F a = e a E + e a c [vB] (17)
n a (t, r) = ´R3 f a (t, r, v) dv v a (t, r) = 1 n a ´R3 vf (t, r, v) dv = v a T a (t, r) = 1 n a ´R3 m a 3 (v -v a ) 2 f (t, r, v) dv = m a 3 (v -v a ) 2 (18)
which are respectively the a-particles density, average velocity and temperature.

In [START_REF] Gaye | Sliding mode stabilization of the current profile in tokamak plasmas[END_REF] the integrals' domain has been chosen to be R 3 because no restrictions are considered on the admissible velocities. The bracket . denotes the average value of a variable, it should be distinguished with the paring ., . using in covariant form. The third equation is obtained by considering the thermal equilibrium (when v a = 0) m v 2 /2 = 3T /2. In this work, there is no fusion reaction considered. Hence balance equations for one species in the plasma may be usually derived from those for the other species (neglecting the eventual interaction terms). Macroscopic variables related to only one representative species will be considered in this work. Therefore, there will be no more use of the species subscripts a for electrons, ions, and neutral particles. In other words, T, P, v, ... will be representatives of average temperature, pressure, velocity, ... for all species (which usually correspond to electronic values in the plasma control community).

Transport equations (balance equations)

Three macroscopic balance equations may be derived by multiplying the Boltzmann equation ( 16) with 1, mv, and mv 2 2 respectively, and integrating the results over the velocity domain (the detailed derivation, inspired from cf. [10, p.208] may be found in the appendix). Defining the usual material derivative by

d dt = ∂ ∂t + v.∇ (19)
these calculations lead to:

• the equation of continuity (particle transport):

dn dt = -n∇.v (20)
which is obtained assuming the density conservation (no external source term) and the assumption that the third term in ( 16) vanishes rapidly as v → ∞.

• the equation of motion (momentum transport):

n dv dt = -∇P -∇ × τ -en E + 1 c [vB] (21)
where P denotes the scalar pressure, τ the stress tensor, and electron charge e = 1.6 10 -19 C.

• the (total) energy transport equation:

d dt = -∇. [nT v + τ v + q] + en (Ev) ( 22
)
where q is defined as the "heat" flux density, T is the average temperature and the total energy is defined as:

= mn 2 v 2 + 3 2 nT ( 23
)
The internal energy (or heat balance) equation may be obtained when the total energy ( 23) is replaced in equation ( 22), using the continuity and momentum equations:

3 2 n dT dt = -P ∇.v -∇.q -τ (∇ × v) (24)
Defining the entropy in the perfect gas case as s = ln T 3/2 /n = ln P 3/2 /n 5/2 , the heat balance equation may be written:

T n ds dt = -∇.q -τ (∇ × v) (25)
The balance equations ( 20), ( 21) , ( 24) and ( 25) ) will now be written in co-variant form for a representation in the port-Hamiltonian formulation.

Co-variant formulation of the transport equations

The properties and dynamics of the plasma considered as a fluid will be described in this section. The mass, momentum, energy and entropy balance equations will be derived in the co-variant form. We shall follow closely the general frame presented in [START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF]Chap.3,p.154] but express the results in terms of differential forms.

The material properties are defined by the total material energy H M which is the integral over the material domain M of total energy density . It is the sum of the kinetic and internal energy of the system which depends on the specific quantity (per unit mass) of the momentum density p ∈ Λ 2 (M)4 , the entropy s ∈ Λ 3 (M) and the volume v ∈ Λ 3 (M).

This subsection gives a definition of the material derivative in terms of differential form which will be used in the following to write the balance equations on a moving material domain.

Time derivative "following the motion"

The material properties must be described for a fixed amount of mass which (due to some velocity field) occupies a moving material domain. We denote φ t the flow associated with the velocity field v (i.e. v = d/dt (φ t z)). Any moving material domain with a fixed mass may be defined by V (t) = φ t (Ω) where Ω denotes the same material domain at time t = 0. The total variation of some conserved material quantity α k (where α k is a k-form, defined on a n-dimensional spatial domain ) becomes [16, Chap.1, sec.4]:

d dt ˆV (t) α = ˆV (t) ∂α ∂t + L v α (26)
where L v α denotes the Lie derivative of the k-form α. Using Cartan's formula:

L v = i v • d + d • i v ( 27
)
one may define the material derivative Dα Dt in term of the variation of some conserved quantity α such that:

ˆV (t) Dα Dt = d dt ˆV (t) α (28)
Therefore:

Dα Dt := ∂α ∂t + i v dα + d(i v α) (29)
For a volume form α n in an n-dimensional spatial domain (dα n = 0) and ( 29) becomes:

D Dt α = ∂α ∂t + di v α (30)
which is equivalent to the usual ∂ t + v • ∇ vectorial form.

Balance and closure equations

A balance equation for the material quantity α in a fixed frame generally takes the form:

∂α ∂t = -df α + σ α (31)
where f α defines the flux of α, and σ α represents the external α source (by unit of mass, since α is a density). Therefore, when considering a moving material domain, the balance equation becomes (using ( 30) and ( 31)):

D Dt α =-df R α + σ α ( 32
)
where f R α = f α -i v α is now the relative flux of α and D Dt the material derivative. We can now express in co-variant form all the balance equations obtained from the kinetic theory in a moving material domain.

Mass balance (continuity equation)

.

∂ρ ∂t = -df ρ + σ ρ ( 33
)
or using ( 20) and ( 32):

D Dt ρ = -df R ρ + σ ρ ( 34
)
with the mass density ρ = mn, the relative particle flux df R ρ = d (f ρ -i v ρ), and the material source term σ ρ per unit volume.

Remark 4.1. The same continuity equation may be written as a massic volume balance equation. Let v := 1/ρ denotes the specific massic volume and σ v as the massic source term, then from equation ( 34) we obtain

D Dt v = d (i v v) + σ v ( 35
)
with the relative volume flux df R v ≡ d (-i v v). Momentum balance. 1 v Dp Dt = -df R p + σ p ( 36
)
where p is the momentum per unit mass and σ p is the source term caused by external body force per unit mass (the electromagnetic pressure in this work).Note that one could also use a more usual identification of both the velocity v and the current J as vector fields, leading to an alternative expression for the electromagnetic pressure:

σ p = -i ( J) B ( 37
)
Energy balance.

D Dt = -df R + σ ( 38
)
where the source term σ represents the inductive (ohmic) heat source term which is v (E ∧ J).

Closure equations and entropy production

In this subsection, we consider the closure equations for the momentum and energy balance equations. Then, since the entropy is not a conserved quantity, the irreversible entropy production term must be computed.

Momentum

We first consider the momentum closure equation. According to [START_REF] Hamroun | A port-controlled hamiltonian approach to geometric reduction of distributed parameters systems -application to the shallow water equations[END_REF], the tensor f R p , representing the momentum flux per unit of surface area, is the sum of a pressure P and viscous force terms τ :

f R p = (P I + τ ) ( 39
)
Considering that P is a scalar (in an isotropic domain), while τ ∈ Λ 1 is the non-diagonal part of stress tensor, or the so-called shear stress tensor, we have:

df R p = dP + dτ ( 40
)
D Dt s = 1 T Du Dt + P ∧ Dv Dt ⇔ -df R s + σ s = 1 T -df R u + σ u + P ∧ d (i v v) ⇔ -d f q T + σ s = 1 T (-df q -(P ∧ d (i v v) + τ ∧ d (i v v))) + + 1 T (v ( η J ∧ J tot ) + P ∧ d (i v v)) (48)
where the heat flux f q is related to the entropy and internal energy fluxes:

   f R u = f q f R s = 1 T f q (49)
Thus the irreversible entropy source term σ s is obtained with the help of the internal energy source term in equation ( 45):

σ s =d f q T + 1 T -df q -τ ∧ d (i v v) + g ∧ f R ρ = - 1 T 2 f q ∧ dT - 1 T τ ∧ d (i v v) + v T η J ∧ J (50)
This irreversible entropy source contains the heat conduction term -1 T 2 f q ∧ dT , the viscous dissipation term -1 T τ ∧ d (i v v) and the Joule (ohmic) term v T η J ∧ J. The three of them define the constitutive relations for the R -field (see next section).

Multidomain couplings and resistivity field (R -field)

In this section, all power conservative multi-domain couplings are detailed (and represented using Dirac interconnection structures [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF]) as well as a dissipative R -field multi-domain coupling element.

Magneto-hydrodynamic coupling (MHD coupling)

The MHD coupling is the central multi-domain couplings in the tokamak system. It may be represented using a magnetomotive Dirac structure combined with an Eulerian-Lagrangian coordinate transformation used to relate the volumetric and massic domains.

Magnetomotive coupling

Let us now prove that the two relations ( 14) and ( 37) may be associated with a Dirac structure. We will make use of the flow and effort variables identifications (f J , e J ) = (J, E L ) and (f v , e v ) = (σ p , p). Let

F 1 = Λ 2 (Ω) × Λ 1 (Ω) (f J , f v ) and E 1 = Λ 1 (Ω) × Λ 2 (Ω) (e J ,
e v ) be respectively the spaces of conjugated flow and effort variables. For a given magnetic induction field B, we can note that equations ( 37) and ( 15) define two anti-adjoint maps:

Λ 2 (Ω) → Λ 1 (Ω) f J → f v = -( f J ∧ B) and Λ 1 (Ω) → Λ 2 (Ω) e v → ( e v ∧ B)
Proposition 5.1. For any given B ∈ Λ 2 (Ω), the linear subspace D 1 (B)

D1 (B) =      (fJ , fv, eJ , ev) ∈ F1×E 1 | fv = -( fJ ∧ B) eJ = ( ev ∧ B)      (51)
is a Dirac structure in the bond space in F 1 ×E 1 .

Let us first precise that the Dirac structure is defined with respect to the symmetrized pairing of the power pairing on F 1 ×E 1 defined by:

e J e v , f J f v = ˆΩ e v ∧ f v -e J ∧ f J First note that D 1 (B) ⊂ D 1 (B) ⊥ (see subsection 2.3) is equivalent to: e e v , f J f v = 0 , ∀ (f J , f v , e J , e v ) ∈ F 1 ×E 1
which easily follows from [START_REF] Ou | Optimal tracking control of current profile in tokamaks[END_REF]. Now let us prove that D 1 (B)

⊥ ⊂ D 1 (B). For this, assume that (f J , f v , e J , e v ) ∈ F 1 ×E 1 satisfies ∀ (f J , f v , e J , e v ) ∈ D 1 (B): e J e v , f J f v + e J e v , f J f v = 0 Then, equivalently ∀ (f J , e v ) ∈ Λ 2 (Ω) × Λ 1 (Ω): ( e v ∧ B) e v , f J f v + e J e v , f J -( f J ∧ B) = 0
That is:

ˆΩ (e v ∧ f v -( e v ∧ B) ∧ f J ) -(e J ∧ f J -e v ∧ ( f J ∧ B)) = 0 (52) ∀ (f J , e v ) ∈ Λ 2 (Ω) × Λ 1 (Ω). Using the fact that ´Ω α ∧ β = ´Ω α ∧ β that α = α one gets: ´Ω ( e v ∧ B) ∧ f J = ´Ω e v ∧ B ∧ f J = - ´Ω e v ∧ f J ∧ B = -´Ω e v ∧ ( f J ∧ B)
Choosing f J = 0, the preceding equation becomes:

ˆΩ e v ∧ (f v + ( f J ∧ B)) = 0 , ∀e v ∈ Λ 1 (Ω) hence f v = -( f J ∧ B).
In a similar way, when choosing e v = 0, (52) becomes:

ˆΩ (e J ∧ f J + e v ∧ ( f J ∧ B)) = 0 , ∀f J ∈ Λ 2 (Ω)
or with similar calculus:

´Ω e v ∧ ( f J ∧ B) = -´Ω e v ∧ ( B ∧ f J ) = - ´Ω e v ∧ B ∧ f J = -´Ω ( e v ∧ B) ∧ f J = - ´Ω ( e v ∧ B) ∧ f J ˆΩ (e J -( e v ∧ B) ∧ f J ) ∧ f J = 0 , ∀f J ∈ Λ 2 (Ω)
yielding e J = ( e v ∧ B). Hence (f J , f v , e J , e v ) ∈ D 1 (B) which completes the proof.

Eulerian-Lagrangian transformation

We shall define hereafter a Dirac structure associated with the change from Eulerian to Lagrangian coordinates (that is from integration on a fixed domain Ω to a moving domain M = φ t (Ω)). Therefore let us write the power product between the effort and the flow variables in the two coordinates and look for the following equality:

ˆM=φt(Ω) e 2 ∧ f 2 = ˆΩ φ * (e 2 ∧ f 2 ) = ˆΩ φ * (e 2 ) ∧ f 1 (53)
with φ * (e 2 ) = 1/ ve 2 •φ -1 , where v is the specific massic volume of the plasma. This defines the following vector subspace :

D (v) =    (f 1 , f 2 , e 1 , e 2 ) ∈ Λ 1 (Ω) × Λ 1 (M ) × Λ 2 (Ω) × Λ 2 (M ) f 2 = f 1 • φ -1 ; e 1 = φ * (e 2 ) • φ; f 1 ∈ Λ 1 (Ω) , e 2 ∈ Λ 2 (φ (Ω))    (54)
The subspace D (v) (as well as the subspace resulting from the interconnection of D 1 and D (v)) may be proved to be Dirac structure using the power product in the material domain suggested in [START_REF] Walker | Emerging applications in tokamak plasma control[END_REF] and following arguments similar to those in [START_REF] Van Der Schaft | Hamiltonian formulation of distributed parameter systems with boundary energy flow[END_REF]. In the tokamak model the pairs of flow and effort variables connected to the Dirac interconnection structures D 1 and D (v) are respectively (f 1 , e 1 ; f 2 , e 2 ) 1 = (-i J B, p; -i J B, i v v) and (f 1 , e 1 ; f 2 , e 2 ) 2 = (J, η J; J, η Jv).

Interconnection structures in the material domains

The material domains are the mechanical, hydraulical and thermal domains. The pairs of energy variables in these subdomains are interconnected together through Stokes-Dirac structures of the form:

D i =            (f p , f q , f i∂M , e p , e q , e i∂M ) ∈ F i × E i | f p f q = 0 d d 0 e p e q , f i∂M e i∂M = 0 1 1 0 e p e q | ∂M            i = {2, 3, 4} (55)
• The Stokes-Dirac structure D 2 for the hydraulical domain is:

(-dP, d (i v v) , f 2∂M , i v v, -P, e 2∂M ) = (f p , f q , f 2∂M , e p , e q , e 2∂M ) ∈ F 2 × E 2 (56)
where

F 2 = Λ 1 (M ) × Λ 3 (M ) × Λ 0 (M ) E 2 = Λ 2 (M ) × Λ 0 (M ) × Λ 2 (M )
• The Stokes-Dirac structure D 3 for the mechanical domain may be written:

(dτ, d (i v v) , f 3∂M , i v v, τ, e 3∂M ) = (f p , f q , f 3∂M , e p , e q , e 3∂M ) ∈ F 3 × E 3 (57)
where

F 3 = Λ 2 (M ) × Λ 2 (M ) × Λ 1 (M ) E 3 = Λ 1 (M ) × Λ 1 (M ) × Λ 1 (M )
• The Fourrier heat transport in the thermal domain may also be represented by using a Stokes-Dirac structure D 4 of the form (see for instance [START_REF] Baaiu | Port-based modelling of mass transport phenomena[END_REF], here the exterior derivative of T is a spatial gradient):

d f q T , dT, f 4∂M , T, f q T , e 4∂M = (f p , f q , f 3∂M , e p , e q , e 3∂M ) ∈ F 4 × E 4 ( 58
)
where

F 4 = Λ 3 (M ) × Λ 1 (M ) × Λ 2 (M ) E 4 = Λ 0 (M ) × Λ 2 (M ) × Λ 0 (M )

Multi-domains fields elements

In this section we formally introduce a dissipative R -field accounting for the irreversible entropy production (50) and a conservative C -field which represents the first principle in the material domain with the energy density given in [START_REF] Jamiolkowski | Classical Electrodynamics[END_REF]. We follow a classical approach from irreversible (near equilibrium) thermodynamics [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] and apply it to model transport phenomena in toroidal plasmas. Therefore we will assume the existence of a (semi) linear relation between the generalized flows and thermodynamical forces.

Resistive field (R -field)

According to [START_REF]Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF], the irreversible entropy production σ s creation is contributed by different resistivity sources: the heat conduction, the ohmic diffusion and the viscosity, with associated dissipated power:

T σ s = - 1 T f q ∧ dT -τ ∧ d (i v v) + v( η J ∧ J) (59)
According to [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] and [START_REF] Boozer | Onsager symmetry of transport in toroidal plasmas[END_REF] or [START_REF] Garbet | Thermodynamics of neoclassical and turbulent transport[END_REF], there is a transport matrix Γ which yields the "linear" relation between efforts and flows variables:

f = Γde (60)
Onsager's reciprocal relations [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] implies that the transport matrix Γ is symmetric, for some special choice for the flows variables (and consequently of their conjugated efforts). As a first approximation, for the R-field definition, we may consider (as in [START_REF] Boozer | Onsager symmetry of transport in toroidal plasmas[END_REF]), that the transport matrix Γ

  f q τ J   =   χ • • • κ • • • 1 /η   Γ d   T i v v V loop   (61)
to be diagonal and the diagonal terms χ , κ , 1 /η to represent respectively the thermal diffusion coefficient, the viscosity, and the electrical resistivity. V loop is the loop voltage such that the electrical intensity is E Ω = dV loop . This constitutive relation (59) completes our model taking into account the irreversible entropy production (dissipation) in the R-field. The off-diagonal terms in (61) may play significant roles in the plasma dynamics. Nevertheless, the are until now poorly known and the subject of active researches [START_REF] Garbet | Thermodynamics of neoclassical and turbulent transport[END_REF].

Energy storage element (C-field)

The energy storage element in the material domain may be define according to equation (41):

= p ∧ p 2 + u (s, v) (62)
where the power contributions come from the mechanical domain via kinetic

energy (i v v ∧ 1 v Dp Dt
), from the thermal domain with (T ∧ Ds Dt ) and from the hydraulical domain by (P ∧ Dv Dt ) (see the Gibbs-Duhem equation ( 46)). The blue part describes the electromagnetic domain (fixed volume balance equations). The Dirac structure D EM in [START_REF] Courant | Dirac manifolds[END_REF] couples the electric and magnetic domains and stands for the usual Maxwell's equations [START_REF] Dalsmo | On representations and integrability of mathematical structures in energy-conserving physical systems[END_REF] together with the 0-junctions (which represents the magnetic and electric densities balance equations). The 1-junction stands for the electrical field balance equation ( 15) coupling together Lorentz forces ( 14) and Ohm's law. The electromagnetic energy is accumulated in the energy storage elements C :: H (D) and C :: H (B) whose constitutive equations are given by [START_REF] Choquet-Bruhat | Analysis Manifolds and Physics[END_REF]. It is dissipated by the Joule effect via the R -field element with a constitutive equation being part of (61). Two usual controlled variables in the electromagnetic domain are "loop voltage" (or total plasma current) from the boundary pair (V loop , I p ) (generated by the main magnetic coils) and the non-inductive current density from the distributed pair of variables(J ni , η J ni ) (usually generated by antennas).

Integration of the complete model with the port-based approach

The red, green, and black parts in the Bondgraph respectively sketch the models for the thermal, mechanical and hydraulic subdomains in the moving frame material domain. The energy storage in the material domain This modular multi-physics model allows a separate analysis for (some) interesting subdomain(s) with the appropriate simplifying assumptions. For inbetween energy storage subsystems (i.e. C :: H elements), dissipative subsystems (i.e. Rfield elements) , and power preserving interconnection elements (such as Dirac structures D or junctions 0 and 1). Power exchanges between the system energy and the environment are "flowing" through boundary or distributed "external" ports). Each power flow in each domain is identified as a pairing between effort and flow variables denoted e ∧ f . stance, if only the electromagnetic domain is considered and when axisymmetry and quasistatic equilibrium assumptions are made, one obtains the classical 1D resistive diffusion equation commonly used for the control of the poloidal magnetic flux radial profile (see [START_REF] Vu | Structure preserving reduction for thermo-magneto plasma control model[END_REF] for this reduction and equivalence with the resistive diffusion equation). Besides, the stability of Dirac structures under serial, parallel and feedback interconnection allows guarantees the passivity and allows passivity based control approaches for the interconnected model. The balance and closure equations in the moving material domain are derived from the Boltzmann equation, using the classical kinetic theory approach. A covariant port-Hamiltonian formulation of these equations is obtained by introducing material derivatives. The entropy balance and the Gibbs-Duhem equations allowed us to compute the irreversible entropy production (source) terms. These entropy production terms and transport phenomena have been used to define TMHD couplings including a transport R -field (which satisfies Onsager's principle) and Stokes-Dirac interconnection structures.

Now that the global model interconnection structure and constitutive equations are determined, one expects to be able to gradually augment the usual resistive diffusion equation control model (which corresponds only to the electromagnetic balance equations with some restricting assumptions). The first step will be to include a simple thermal model with the heat production from the Joule effect (ohmic resistivity), a simple heat transport diffusion (dispersion) model and the variation of the resistivity with the temperature. This model will allow us to consider the important effect of temperature variations on the resistivity and also to use the non inductive actuators for the heating of the plasma. The second step will be to consider the magneto-hydrodynamic and hydrodynamic part to obtain a control model which includes the bootstrap current effect. Finally, the third step will be the inclusion of the mass balance equation and the supplementary terms and equations related to the fusion reaction.

AppendixA. Microscopic model (transport equations)

Using usual tensor notations, the Boltzmann equation for distribution functions f a (t, r, v) for particles of species a with position r and velocity v in the coordinates β reads: as given in equation [START_REF] Hamroun | Port-based modelling for open channel irrigation systems[END_REF].

AppendixA.2. Equation of motion (momentum transport)

From now, we denote respectively m, n , v the corresponding quantities m a , n, v a for the species a. Multiplying the Boltzmann [START_REF] Frankel | The Geometry of Physics : an Introduction[END_REF] with e a = e e = -e for electrons. Decomposing the velocity field v = v + v , with the mean velocity v and a random velocity v whose v = 0, then:

vv β = vv β + v v β (A.7)
and the second term of (A. where P denotes the scalar pressure and τ the stress tensor, whereas: In first term, the total energy (which includes the kinetic energy) may be written with the help of (A.10):

P = nm v 2 /3 = nT τ aβ = nm v v β -v 2 /
= mn 2 v 2 = mn 2 v 2 + mn 2 v 2 = mn 2 v 2 + 3 2 nT (A.14)
The second term in (A.13) corresponds to the divergence of the energy flux:

v 2 v β = (v + v ) 2 v β + v β = v 2 v β + v 2 v β + 2v v v β + (v ) 2 v β = v 2 v β + (3nT ) v β nm + 2 (τ + nT ) v nm + 2q nm (A.15)
Hence thanks to the definition of the pressure P and the stress tensor τ in (A.10), using the flux density of heat:

q = ˆR3 m 2 v 2 vf dv = mn v 2 2 v
the energy balance equation in (A.13) may be written:

Figure

  Figure 1: Schematic view of a tokamak with the electrical solenoids: the magnetic field generated by the three magnets makes the plasma gas ions following helicoidal trajectories along the torus

Figure 2 :

 2 Figure 2: Model structure for the developed plasma port-Hamiltonian model. The electromagnetic balance equations are computed on a fixed volume integration domain while mass, entropy and momentum balance equations must be computed in a moving material domain. Lorentz forces and the Joule effect are the main couplings between these electromagnetic and "material" domains.

  is caused by electric and magnetic fields, respectively E and B. The collision term is usually of the form C a = b C ab (f a , f b ) where C ab denotes the rate of change for the distribution function of particles a due to elastic and/or inelastic collisions with b. This collision term is neglected in this work. It is a quite usual assumption since the plasma gas in tokamaks is low density and usually considered as a perfect gas (cf.[10, p.206] or[START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF] Chap.6]). Averaging momenta of the distribution function over all possible velocities, the average or macroscopic classical (fluid-like) quantities may be derived:

Figure 3

 3 Figure 3 presents the complete tokamak 3D model, using the Bondgraph 5 (cf.[50, Chap.2, p.53]) notations. It includes all energy domains, as well as their

  is represented by a field element C :: H = ´M p ∧ p 2 + u(s, v) with the internal energy given in the Gibbs-Duhem equation (46). Commonly the ideal gas constitutive equation is used in the material domain since the plasma density is very low. All the balance equations developed in section 4 are realized with the 0-junctions in the three subdomains. The transport equations (61) and the irreversible entropy source terms (50) are realized in the dissipative R-field which transfers all the dissipated energy to the thermal domain. In each material subdomain, there is a Stokes-Dirac structure (see section 5.2) such as the structure D 2 for the hydraulical domain (equation (56)), D 3 for the mechanical domain (equation (57)) and D 4 for the thermal domain (equation (58)). In each of these Stokes-Dirac structures, boundary port variables are available for connections with external sensors/actuators systems or more generally with the external environment. Two other used control actions in the material domain, common and important in tokamak applications, are the distributed heating source S heat in thermal domain and the material injection source σ v in the hydraulical domain. The electromagnetic volumetric balance equations in the blue part are related to the material domain through the Dirac structure D (v) defined in (54) and the Eulerian-Lagrangian transformation T F :: B = D 1 as in (51).

Figure 3 :

 3 Figure 3: Bondgraph of thermodynamic tokamak system including electromagnetic domain (blue), mechanical domain (green), thermal domain (red) and hydraulic domain (black)

1 )=

 1 For the plasma example in the tokamak we may consider the collisionless case (C a = 0) and external forces F a = e a E + e a c[vB]. We perform then the integration of the products of that equation with 1, m a v a , and m a v 2 a 2 respectively, over the velocity domain, in order to deduce the mass, momentum, and energy balances (cf.[10, p.208]).AppendixA.1. Equation of continuity (particle transport equation)When integrating equation (16) over the velocity domain, we get the particle transport equation: a = e or i. In this last expression one assumes that the third term of (16) vanishes rapidly as v → ∞ and the resulting term in the integration may be neglected when compared to the two first terms. Omitting the species subscripts, one gets: ∂n ∂t + ∇. (nv) = 0 (A.3) Using the usual material derivative: -∇. (nv) + v∇.n = -n∇.v -v∇n + v∇.n = -n∇.v (A.5)

  3) is:∇. (mn vv β ) = ∇. (mnvv β ) + ∇. mn v v β = mnv β ∇. (v) + mv∇. (nv β ) + +∇. mn v v β (A.8)One obtains then the pressure tensor in the form:P aβ = ˆR3 mv v β f dv = mn v v β = P δ aβ + τ aβ (A.9)

  ( v β ) = -∇. [( + nT ) v β + τ v + q] + en (Ev) + ∇. ( v β ) = -∇. [nT v β + τ v + q] + en (Ev) (A.17

  )

Table 1 :

 1 The vector notation form and the co-variant formulation of the characteristic equations of electrodynamics. The Hodge star products , µ are associated with the permittivity and permeability tensors respectively.

  by mv, and integrating over the velocity domain, one gets:

		m	ˆR3	v	∂f a ∂t	+	∂ ∂r β	(v β f a ) +	∂ ∂v β	F aβ m a	f a	dv=0
												(A.6)
	⇔	∂ (mnv) ∂t	+ ∇. (mn vv β ) -e a n E a +	1 c	[vB] a =0

  3 δ aβ (A.10) On the other hand, thanks to the continuity equation in (A.3), the first term of equation (A.6) is derived:

	AppendixA.3. Energy transport
	Multiplying the Boltzmann equation (16) by velocity domain, one gets:	mv 2 2	and integration over the
		ˆR3	mv 2 2		∂f a ∂t	+	∂ ∂x β	(v β f a ) +	∂ ∂v β	F aβ m a	f a	dv =0
											(A.13)
	⇔	∂ ∂t	mn 2	v 2 + ∇.	mn 2	v 2 v β -en (Ev)	=0
					∂ (mnv a ) ∂t		= mn	∂v ∂t	+ mv	∂n ∂t	(A.11)
										= mn	∂v ∂t	-mv [∇. (nv β )]
	Then from (A.15-A.11), and (A.4) the momentum transport can be rewritten
	as:		nm	dv dt	= -∇P -∇ × τ + en E a +	1 c	[vB] a	(A.12)

The bootstrap current comes from an magnetohydrodynamical effect and has been considered in the plasma control literature (until now) as a perturbation. However it is now hoped to be one of the major source of inductive heating in the ITER project.

The conjugated variable of the momentum density p ∈ Λ 1 (M) is a 1-form equivalent to the average plasma velocity v in[50, section 3.4]

Bondgraph is a graphical modelling methodology which applies for "first principle" multidomain dynamical systems. It is based on the explicit representation of the power transfers

Internal energy

The total energy density (per unit mass) is:

with u is the internal energy. Its time variation is defined according to:

where the contraction i v v represents to the plasma velocity v ( p) of a plasma volume v. The material derivative of internal energy u is written following the general formulation [START_REF] Nishida | Multiscale distributed porthamiltonian representation of ionic polymer-metal composite (ipmc)[END_REF]:

Using the momentum balance equation [START_REF] Ou | Robust control design for the poloidal magnetic flux profile evolution in the presence of model uncertainties[END_REF] and the internal energy balance equation ( 43), the total energy balance(42) becomes:

with

Irreversible entropy production

The internal energy is also determined by the Gibbs-Duhem equation:

which will allow us to compute the irreversible entropy production source term σ s (without nuclear fusion reaction in the studied case). The entropy production rate is then:

Applying the material derivative formula for s, u and v in (47), we get:

AppendixA.4. Equation of internal energy and entropy equation

The internal energy balance may be derived by using the equations of continuity and motion here above to eliminate the kinetic energy in equation (A. [START_REF] Frankel | The Geometry of Physics : an Introduction[END_REF]). We first remind the kinetic rate term:

due to (A.5) and (A.12), and the fact that v [vB] = 0 . Then (A.17) becomes:


Shlomo Barak
email: shlomo@tagapro.com

Equations of Physics and Space Density

Keywords: Space, Navier, Maxwell, Einstein, General Relativity, Ricci Flow, GDM

We show that the "hidden" variable of the Navier, Maxwell, Einstein and Ricci Flow equations is  -the space lattice density. In our GeometroDynamic Model of the physical reality -the GDM (presented in a series of papers during 2017) -we show that sources are merely highly deformed zones of space (highly contracted or dilated) whereas their fields are continuations of these deformations, reduced by orders of magnitudes, with only a small deviation from the normal (standard) space density. Thus the equations of physics express different aspects of the same reality -the geometrodynamics of space.

Introduction 1

Space as a Lattice 1. [START_REF] Barak | The Geometrodynamics of Space[END_REF] By attributing a cellular structure to space we can explain its expansion, its elasticity and can introduce a cut-off in the wavelength of the vacuum state spectrum of vibrations [START_REF] Barak | The Geometrodynamics of Space[END_REF].

Space density  1.2

Space density  is defined as the number of space cells per unit volume [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF]. Space density in a zone of space without deformations (far away from masses and charges) is denoted  0 .

Let dn be the number of space cells in a given volume dV. Since dn = ρ 0 dV and also In Appendix A of [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF] we prove that the relative volumetric change equals the divergence of the Elastic Displacement Vector u:

dV dV dV'    u  (1)
Thus:

ρ ρ ρ 0    u  ( 2
)
This proof is a cornerstone of the GDM. It is based on the strain tensor u ik , which is defined as:

                   k i l l i k k i ik x x u u x u x u 2 1 u
The Small Deformation Strain Tensor as a Fundamental Metric Tensor 1. [START_REF] Rindler | Relativity Oxford[END_REF] The authors of [START_REF] Palacios | The Small Deformation Strain Tensor as a Fundamental Metric Tensor Journal of High Energy Physics[END_REF], see also [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF], conclude: "……. the small deformation strain tensor could be used as a fundamental metric tensor, instead of the usual fundamental metric tensor." Note, however, that space deformation is a local feature of curvature, whereas manifold curvature can also be a global feature of curvature (like that of a closed spherical surface).

Deformed Spaces versus Bent Manifolds 1.5

We relate to space not as a passive static arena for fields and particles but as an active elastic entity, which in our model is the one and only entity that exists. Physicists have different, sometimes conflicting, ideas about the physical meaning of the mathematical objects of their models. The mathematical objects of General Relativity, as an example, are n-dimensional manifolds in hyper-spaces with more dimensions than n. These are not necessarily the physical objects that General Relativity accounts for, and n-dimensional manifolds can be equivalent to n-dimensional elastic spaces. Rindler [START_REF] Rindler | Relativity Oxford[END_REF] uses this equivalence to visualize bent manifolds, whereas Steane [START_REF] Steane | Relativity made relatively easy[END_REF] considers this equivalence to be a real option for a presentation of reality. Callahan [START_REF] Callahan | The Geometry of Spacetime Springer[END_REF], being very clear about this equivalence, declares: "…in physics we associate curvature with stretching rather than bending". After all, in General Relativity gravitational waves [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF] are space waves and the attribution of elasticity to a 3D space is thus a must.

The deformation of space is the change in size, of its cells [START_REF] Barak | On Bent Manifolds and Deformed Spaces[END_REF]. The terms positive deformation and negative deformation, around a point in space, are used to indicate that space cells grow or shrink, respectively, from this point outwards. Positive deformation is equivalent to positive curving and negative deformation to negative curving. In [START_REF] Barak | On Bent Manifolds and Deformed Spaces[END_REF] we show that the curvature K around a point in space is:

π 4 K             (3)
The GDM and the Equations of Physics 1.6

In our "GeometroDynamic Model of reality", the GDM, the variable in the equations of physics is directly or indirectly related to the density ρ of the elastic space lattice. These equations yield wave equations that are also related to the density ρ of the elastic space.

In the Navier equation [START_REF] Barak | The Geometrodynamics of Space[END_REF], the displacement vector u is related to ρ:

∇•u = (ρ-ρ 0 )/ρ.
In the Maxwell equations, [START_REF] Barak | Electromagnetism as the Geometrodynamics of Space[END_REF] and [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF], the charge density q is defined as q = 1/4π•(ρ-ρ 0 )/ρ, and the fields E and B are related to u in the following way: E = Hu and B = 1/c v×E , where v is the velocity of the charge that creates E and hence B.

In the General Relativity (GR) equation, [START_REF] Barak | On the Essence of Gravitation and Inertia B, Part 2: The Curving of Space by an Elementary Particle[END_REF], the variable is the metric g ij and as shown in Section 1.4 it is related to the strain tensor u ik, which is related to  . The Ricci tensor R ij in the GR equation expresses curvature. But the curvature K is also related to space density:

K = 4π/45•(∇ρ/ρ) 2 .
In the Ricci Flow equation:

∂ t g ij = -2R ij (Section 6)
with the Ricci tensor and the metric, the situation is the same.

Necessarily, all waves are space waves of different kinds.

Note that transversal waves do not require a time-variable space density ρ, whereas longitudinal waves require a time-variable space density ρ.

The Navier equation yields both transversal and longitudinal waves, whereas the Maxwell and Einstein equations seem to yield "only" transversal waves. We relate to this situation in Sections 3 and 5.

The GDM Postulate 1.7

The elastic and vibrating three-dimensional Space Lattice is all there is. Elementary particles are transversal or circulating longitudinal wavepackets of this vibrating space [START_REF] Barak | The Geometrodynamics of Space[END_REF].

Navier Equation -The Vibrating Space 2

Elastic Wavesa Reminder 2.1

The equation of equilibrium in an elastic media, with displacement vector, u, (53.6) in [START_REF]VI Elasticity and Plasticity[END_REF], is the Navier equation:

    0 t m 2 2 2            u u u    (4)
 and  are the elastic Lamé coefficients, and m is the media mass density.

Mass in the GDM is only a practical attribute [START_REF] Palacios | The Small Deformation Strain Tensor as a Fundamental Metric Tensor Journal of High Energy Physics[END_REF]. The relevant fundamental attribute is energy U or energy density . Hence, for space as the elastic media, m represents the energy density of vibrations:

m 2 c   (5)
To solve equation (4), we adopt the Kelvin [START_REF]VI Elasticity and Plasticity[END_REF] method and decompose u as follows:

T L u u u   Where: 0 0 T L     u u   Therefore: f T L        u u
 stands for the scalar potential and f for the vectorial potential. This decomposition is true under the boundary condition u → 0 at ∞. The known equations for u L and u T are obtained by substituting

T L u u u   in (4). For u L :   0 t m 2 2 2 2 L L         u u or: 0 t c 1 2 2 2 2 L L L      u u
which is the vector wave equation for waves which move at a speed c L , where:

2 1 m 2 c L           (6) Since 0 L   u 
this is the contractional /dilational, longitudinal wave.

For u T :

0 t m 2 T 2 T 2       u u or: 0 t c 1 2 T 2 2 T T 2      u u
This is again a vector wave equation for waves with speed c T , where:

2 1 T m c         (7)
Since •u T = 0, u T and c T , correspond to the shear, transverse wave. a negative compressibility (1/k) of the aether was required.

On the Transversal and Longitudinal

In the GDM, material (non-zero rest energy) elementary particles are circulating longitudinal wavepackets [START_REF] Barak | On the Essence of Electric Charge B, Part 2: How Charge Curves Space[END_REF], [START_REF] Barak | On the Essence of Gravitation and Inertia A, Part 1: Inertia and Free Fall of an Elementary Particle[END_REF]. This circulation is complex [START_REF] Barak | On the Essence of Electric Charge B, Part 2: How Charge Curves Space[END_REF] and its basic motion is at the longitudinal velocity c L > c T = c. Hence we require:

λ+2μ > μ or λ+μ > 0 see ( 6) and [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF].

In [START_REF] Barak | On the Essence of Electric Charge B, Part 2: How Charge Curves Space[END_REF] we show that c L = 1.6068 c ( 8) hence (λ+2μ)/ μ = 1.6068 2 = 2.5818 and:

λ = 0.5818 μ (9)

Maxwell Equations 3 Classical Electromagnetism 3.1

Electromagnetic theory can be fully expressed by the four Maxwell's equations [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF] and [START_REF] Barak | Electromagnetism as the Geometrodynamics of Space[END_REF]:

I q 4   E  Coulomb's Law II t c 1      B E  Faraday's Law III B   = 0
No magnetic monopoles

IV t c 1 c 4       E j B 

Ampère's Law

In these equations B is the vector [START_REF] Barak | Electromagnetism as the Geometrodynamics of Space[END_REF]:

  E v B   c 1 Magnetic Field (10)
The magnetic field is created by a charge, with the field E, moving at speed v.

The dimensions of B and E are the same, and in the GDM they are:

[E] = [B] = LT -2 , see [2],
[9], and both express the elastic displacement.

In the GDM the charge density q is defined as q = 1/4π•(ρ-ρ 0 )/ρ, and the fields E and B are related to u in the following way: E = Hu and B = 1/c v×E , where v is the velocity of the charge that creates E. The displacement vector u, though, is related to ρ, as follows:

∇•u = (ρ-ρ 0 )/ρ.
The Maxwell electromagnetic equations are, thus, related to space density ρ only.

On electromagnetic longitudinal waves see [START_REF] Heitler | The Quantum Theory Of Radiation[END_REF].

The Einstein General Relativity (GR) Equation 4The Equation 4.1

The right hand-side of Einstein's field equation [START_REF] Barak | The Geometrodynamics of Space[END_REF] of GR (below) should express curvature exactly as the left-hand side does. In [START_REF] Barak | On the Essence of Gravitation and Inertia B, Part 2: The Curving of Space by an Elementary Particle[END_REF] we show, for the first time, that this is indeed the case.

The need to express curvature by Riemannian geometry and obtain a covariant formulation of physical laws, by using tensors, led Einstein to the equation of General Relativity (GR):

ij 4 ij ij ij T c G 8 g Rg 2 1 R       (11)
R ij is the Ricci contracted Riemannian tensor and R is the Ricci scalar, R It was Einstein's vision that "future physics" will show that the right-hand side of ( 1

): c G 8 4  • T ij is also curvature.
In equation (2) of [START_REF] Steane | Relativity made relatively easy[END_REF], R c is the Gaussian radius of curvature, and 1/R c 2 is the Gaussian curvature at the radial distance r from the center of a black hole with Schwarzschild Radius r S (our notations are different from those in [START_REF] Steane | Relativity made relatively easy[END_REF]):

1/R c 2 = r S /r 3 [START_REF] Sommerfeld | Mechanics of Deformable Bodies[END_REF] Equation ( 12) enables us, for the first time to show, [START_REF] Barak | On the Essence of Gravitation and Inertia B, Part 2: The Curving of Space by an Elementary Particle[END_REF], that 8πG/c 4 •T ij expresses curvature.

Curvature 4.2

The Ricci tensor R ij in the GR equation expresses curvature. But the curvature K is also related to space density [START_REF] Barak | On Bent Manifolds and Deformed Spaces[END_REF]:

K = 4π/45•(∇ρ/ρ) 2

The Metric 4.3

In the GR equation the variable is the metric g ij which, as shown in Section 1.4, is related to the strain tensor u ik, which is itself related to .

Thus the GR equation relates to space density only.

The Extended GR Equation of the GDM 5

The Extended Field Equation of General Relativity for Energy/Momentum and 5.1

Charge/Current Using ( 9), [START_REF] Barak | On the Essence of Gravitation and Inertia B, Part 2: The Curving of Space by an Elementary Particle[END_REF] and [START_REF]VI Elasticity and Plasticity[END_REF] in [START_REF] Barak | On the Essence of Gravitation and Inertia B, Part 2: The Curving of Space by an Elementary Particle[END_REF] the extended Einstein GR field equation becomes:

R  -1/2Rg  = 8πG/c 4 • T m  + 4π√G /s 2 • T q  [START_REF] Barak | On the Essence of Electric Charge B, Part 2: How Charge Curves Space[END_REF] And in the GDM, by inserting the constant, H, it takes the form: R  -1/2Rg  = 8πG/c 4 • T m  + 4π√HG /s 2 • T q

 ( 14
)
Electromagnetism relates to space density only, Section 3, and necessarily this is also the case for the extended GR equation.

Note that the gravitational weak field approximation yields the gravitational Poisson equation.

Similarly the electric weak field approximation should yield the electrical Poisson equation.

  = is the energy density of space. In this paper we ignore the Cosmological Constant term Ʌ ij g .

Wave Velocities 2.2

	Historically, to account for the absence of electromagnetic longitudinal waves, Cauchy (19 th
	century) suggested that	λ		2μ		0	, see [12] P.108. Hence	2μ λ  	, but the bulk modulus is
	k				3 2		and since			

We have already explained that both the metric and the curvature are the expression of space density, and necessarily so is also the Ricci Flow.

Summary 7

Equations of physics express different aspects of the same reality. This reality is the geometrodynamics of space, and the same "hidden" variable of these equations is the density of space.

Laws of physics express the geometrodynamics of the elastic space and there is no need to ask where they come from.

Our papers can be found by inserting shlomo barak in the HAL Search: HAL (https://hal.archives-ouvertes.fr)

Acknowledgements

We would like to thank Mr. Roger M. Kaye for his linguistic contribution and technical assistance.

Note that the weak force is electromagnetic (Salam and Weinberg) and according to the GDM so is the strong force.

On gravitational longitudinal waves see the introduction of [START_REF] Fesik | Polarization states of gravitational waves detected by LIGO-Virgo antennas[END_REF].

The Ricci Flow 6

The Ricci Flow in topology, [START_REF] Topping | Lectures on the Ricci Flow March[END_REF], is governed by the equation:

where g ij is the metric of space and R ij is the Ricci tensor. This equation expresses the idea that the rate of change in the metric in a given zone is dependent on the local curvature. This is kind of an "ironing" process of space or a curved manifold, analogues to heat diffusion. We consider this idea in topology to be also relevant to physics. Curved zones of space (contracted or dilated) or bent manifolds have a built-in tension, hence higher energy, due to the elasticity of space. We assume that in order to ease the tension and lower the energy the rate of change, with time, of the metric ∂ t g ij is dependent on the amount of curvature 2R ij . This is the idea of space being "ironed". Thus two zones of space, both curved positively or negatively will repel each other, whereas if one of them is curved positively and the other negatively they will attract each other. This is the case for the bivalent electrical charges [START_REF] Barak | On the Essence of Electric Charge A, Part 1: Charge as Deformed Space[END_REF]. We contend that this should also be the case for masses since they curve space positively. Thus, in contrast to the current understanding, masses should repel each other. Free Fall appears as an attraction force but, as we have seen, it is not.

Locally Free Fall is stronger than the repulsion, but on the large scale of the universe it is the repulsion that takes over.

We know that the Hubble Flow (space expansion) takes place in the inter-galactic space and not within the galaxies themselves. Thus one can argue that galaxies, with their inner black holes, grip space locally and that their mutual repulsion is the reason for space expansion.