- Email: [email protected]
Port-Hamiltonian modeling and reduction of a burning plasma system
6th IFAC Workshop on Lagrangian and Hamiltonian Methods or 6th IFAC Workshop on Lagrangian and Hamiltonian Methods or Nonlinear Control on 6th IFAC IFAC Workshop Workshop Lagrangian and and Hamiltonian Hamiltonian Methods Methods or or 6th Nonlinear Control on Lagrangian Available online at www.sciencedirect.com Universidad Técnicaon Federico Santa María Nonlinear Control 6th IFAC Workshop Lagrangian and Hamiltonian Methods or Nonlinear Control Universidad Técnica Federico Santa María Valparaíso, Chile, May 1-4, 2018 Universidad Técnica Federico Santa Nonlinear Control Universidad Técnica Federico Santa María María Valparaíso, Chile, May 1-4, 2018 Valparaíso, Chile, Chile, May 1-4, 2018 2018 Universidad Técnica Federico Santa María Valparaíso, May 1-4, Valparaíso, Chile, May 1-4, 2018
ScienceDirect
IFAC PapersOnLine 51-3 (2018) 68–73
Port-Hamiltonian Port-Hamiltonian modeling modeling and and reduction reduction Port-Hamiltonian modeling and reduction of a burning plasma system Port-Hamiltonian modeling and reduction of a burning plasma system of a burning plasma system of a burning plasma ∗∗ system ∗∗∗ Benjamin Vincent ∗,∗∗ ∗,∗∗ Trang Vu ∗∗ Nicolas Hudon ∗∗∗
Benjamin Vincent ∗,∗∗ Trang Vu ∗∗ Nicolas Hudon ∗∗∗ ∗∗ ∗ ∗,∗∗ ∗∗∗ Benjamin Vincent Trang Vu Nicolas e Dochain ∗∗ Denis ∗ BenjaminLaurent VincentLef` Trang Vu ∗∗ Nicolas Hudon Hudon Laurent Lef` evre vre Denis Dochain ∗,∗∗ ∗∗ ∗∗∗ ∗∗ ∗ ∗∗ Denis ∗ BenjaminLaurent Vincent Trang Vu Nicolas Hudon Lef` e vre Dochain Laurent Lef` evre ∗∗ Denis Dochain ∗ ∗ Laurent Lef` evre Denis Dochain ee Catholique a ∗ ICTEAM, Universit´ Catholique de de Louvain, Louvain, Bˆ Bˆ atiment timent Euler, Euler, av. av. ∗ ICTEAM, Universit´ ∗ ICTEAM, Universit´ eB-1348 Catholique de Louvain, Louvain, Bˆ Bˆ atiment timent Euler, av. av. Lemaˆ ıtre 4-6, Louvain-La-Neuve, Belgium. ICTEAM, Universit´ e Catholique de a Euler, Lemaˆ ıtre 4-6, B-1348 Louvain-La-Neuve, Belgium. ∗ ICTEAM, Universit´ Catholique de Louvain, Bˆ atiment Euler, av. Lemaˆ ıtre 4-6, 4-6,eB-1348 B-1348 Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) Lemaˆ ıtre Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) ∗∗Lemaˆ ıtre 4-6, B-1348 Alpes, Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) ∗∗ Universit´ (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) Universit´ee Grenoble Grenoble Alpes, LCIS, LCIS, F-26902, F-26902, France. France. ∗∗ ∗∗ Universit´ (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) e Grenoble Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) Universit´ e Grenoble Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) ∗∗∗ ∗∗ Universit´ Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) Chemical Engineering, Queen’s University, ∗∗∗ Department eofGrenoble (e-mail:laurent.lefevre,[email protected]) Department of Chemical Engineering, Queen’s University, ∗∗∗ ∗∗∗ (e-mail:laurent.lefevre,[email protected]) Department of Chemical Engineering, Queen’s University, University, Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) Department of Chemical Engineering, Queen’s Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) ∗∗∗ Department Chemical Queen’s University, Kingston, ON, 3N6 Canada. (e-mail: Kingston, ON, K7L K7L of 3N6 Canada.Engineering, (e-mail: [email protected]) [email protected]) Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) Abstract: Abstract: In In this this contribution, contribution, we we develop develop aa structured structured port-Hamiltonian port-Hamiltonian model model for for aa class class Abstract: In this we structured port-Hamiltonian model class of parameter systems exploit obtained for Abstract: In distributed this contribution, contribution, we develop develop structured port-Hamiltonian model for for class of irreversible irreversible distributed parameter systems aaand and exploit the the obtained formulation formulation for aamodel model Abstract: In distributed this contribution, we develop aand structured port-Hamiltonian modelis for class of irreversible distributed parameter systems and exploit the obtained formulation for amodel model reduction from dimension 3 to dimension 1 in space. The proposed methodology motivated of irreversible parameter systems exploit the obtained formulation for reduction from dimension 3 to dimension 1 in space. The proposed methodology is motivated of irreversible parameter systems and exploit the obtained formulation for model reduction fromdistributed dimension 3 to to dimension dimension 1 in in space. The proposed methodology is motivated by the control of burning plasma profiles in Tokamak reactors. The burning plasma is viewed reduction from dimension 3 1 space. The proposed methodology is motivated by the control burning plasma profiles in Tokamak reactors. The methodology burning plasma is viewed reduction from of dimension 3built to dimension 1 in space. The proposed ismomentum, motivated by the control of burning plasma profiles in Tokamak reactors. The burning plasma is as system on Maxwell equations and total species, by a themulti-physics control of burning plasma profiles in Tokamak reactors. Themass, burning plasma is viewed viewed as a multi-physics system built on Maxwell equations and total mass, species, momentum, by the control of burning plasma profiles in Tokamak reactors. The burning plasma is viewed as a multi-physics system built on Maxwell equations and total mass, species, momentum, energies, and entropy balance equations. Moreover, the system presents nonlinear couplings, as a multi-physics system builtequations. on Maxwell equations totalpresents mass, species, Moreover, the and system nonlinearmomentum, couplings, energies, and entropy balance as a multi-physics system built on Maxwell equations totalpresents mass, species, energies, and entropy balance equations. Moreover, the and system presents nonlinearmomentum, couplings, especially through transport coefficients, and evolves over time scales. energies, and entropy balance equations. Moreover, the system nonlinear couplings, especially through transport and its its dynamic dynamic evolves over multiple multiple time scales. energies, and entropy balancecoefficients, equations. Moreover, the system presents nonlinear couplings, especially through transport coefficients, and Joule its dynamic dynamic evolves overforces, multiple time scales. The main couplings considered here are the effect, the Lorentz and the fusion especially through transport coefficients, and its evolves over multiple time scales. The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion especially through transport coefficients, and its dynamic evolves overrely multiple time scales. The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction on the use of the The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction rely on the use of the The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction rely on the use of the Gibbs relation, Onsager linear transport theory, Stokes–Dirac structures, and energy preserving reaction kinetics. The port-based modeling formulation and structures, reduction rely on the preserving use of the Gibbs relation, Onsager linear transport theory, Stokes–Dirac and energy reaction kinetics. The port-based modeling formulation and structures, reduction rely the preserving use of the Gibbs relation, relation, Onsager linear transport transport theory, Stokes–Dirac structures, and on energy preserving geometric reduction. Gibbs Onsager linear theory, Stokes–Dirac and energy geometric reduction. Gibbs relation, Onsager linear transport theory, Stokes–Dirac structures, and energy preserving geometric reduction. geometric reduction. © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. geometric reduction. Keywords: Keywords: Port-Hamiltonian Port-Hamiltonian systems, systems, Tokamak Tokamak reactors, reactors, Modeling, Modeling, Model Model reduction, reduction, Keywords: Port-Hamiltonian systems, Tokamak reactors, Modeling, Model Irreversible thermodynamics Keywords: Port-Hamiltonian systems, Tokamak reactors, Modeling, Model reduction, reduction, Irreversible thermodynamics Keywords: Irreversible thermodynamics IrreversiblePort-Hamiltonian thermodynamics systems, Tokamak reactors, Modeling, Model reduction, Irreversible thermodynamics 1. INTRODUCTION INTRODUCTION model 1. model has has been been proposed proposed (Vu (Vu et et al., al., 2016), 2016), and and then then 1. model has has been proposed proposed (Vu et al., 2016), and then then reduced. It should be noted that the model considered in 1. INTRODUCTION INTRODUCTION model been (Vu et al., 2016), and reduced. It should be noted that the model considered in 1.promising INTRODUCTION model has been did proposed (Vu al., 2016), and then reduced. It 2016) should be not noted that et the model considered in (Vu et al., consider the material field, i.e., reduced. It should be noted that the model considered in Nuclear fusion is a source of renewable energy et al., 2016) did not consider the material field, i.e., Nuclear fusion is a promising source of renewable energy (Vu reduced. It should be balance noted that the model considered in (Vu et al., 2016) did not consider the material field, i.e., the mass and species equations, and by extension Nuclear fusion is a promising source of renewable energy (Vu et al., 2016) did not consider the material field, i.e., able to solve present and future energy production issues. Nuclear fusion is a promising source of renewable energy mass and species balance equations, and by field, extension able to solve present and future energy production issues. the (Vu et al., 2016) did not consider the material i.e., mass and species balance equations, and by extension Nuclear fusion is a promising source of renewable energy the fusion reaction, were not considered. able to solve present and future energy production issues. specieswere balance equations, and by extension Tokamak reactors areand onefuture of the theenergy investigated technology able to solve present production issues. the mass fusionand reaction, not considered. Tokamak reactors are one of investigated technology specieswere balance equations, and by extension the mass fusionand reaction, were not considered. considered. able to solve present production issues. Tokamak reactors are one of the investigated technology fusion reaction, not to achieve fusion reaction; however, many technical chal- the Tokamak reactors areand onefuture of theenergy investigated technology to achieve fusion reaction; however, many technical chalIn paper, we aa structured the fusion reaction, were not considered.model Tokamak reactors are one before of the energy investigated technology In this this paper, we propose propose structured model for for burning burning to achieve fusion reaction; however, many technical challenges remain to be solved can be produced to achieve fusion reaction; however, many technical challenges remain to be solved before energy can be produced plasma In this paper, we propose aa structured model for burning in Tokamak reactors, i.e., aa structured model takIn this paper, we propose structured model for burning to achieve fusion reaction; however, many technical chalplasma in Tokamak reactors, i.e., structured model taklenges remain to be solved before energy can be produced using this class of reactors. One of the remaining problem lenges remain to be solved before energy can be produced this paper, wethe propose a reaction. structured model for burning using this class of reactors. One of the remaining problem In plasma inaccount Tokamak reactors, i.e., aa structured structured model taking into fusion This reaction occurs plasma in Tokamak reactors, i.e., model taklenges remain to be solved before energy can be produced ing into account the fusion reaction. This reaction occurs using this class of reactors. One of the remaining problem concerns the stabilization of One plasma around desiredproblem profiles plasma using thisthe class of reactors. of the remaining inaccount Tokamak reactors, i.e., a structured model takconcerns stabilization of plasma around desired profiles ingainto into account the fusion reaction. Thiscannot reaction occurs at microscopic level of description and be added ing the fusion reaction. This reaction occurs using thisthe class of reactors. of remaining ainto microscopic levelfusion of description and cannot be occurs added concerns stabilization of plasma around desired profiles such that the fusion reaction is the sustained. To problem achieve concerns the stabilization of One plasma around desired profiles at ing account the reaction. This reaction such that the fusion reaction is sustained. To achieve at a microscopic level of description and cannot be added directly to principle model at a a microscopic of description and cannot added concerns stabilization of plasma aroundcurrent desired profiles at directly to aa first firstlevel principle model developed developed at be a macromacrosuch that the fusion reaction is To achieve this goal, control of plasma plasma temperature, density, such goal, thatthe the fusion reaction is sustained. sustained. To density, achieve at a microscopic level of description and cannot be added this control of temperature, current directly to aaoffirst principle model developed at a macroscopic level description. A formulation for the burning directly to first principle model developed at a macrosuch that the fusion reaction is sustained. To achieve scopic level of description. A formulation for the burning this control of density, and particle density profilestemperature, is required. required. current Actuators and directly this goal, goal, control of plasma plasma temperature, current density, to aof principle developed macroand particle density profiles is Actuators and scopic level level offirst description. Amodel formulation for at theaburning burning plasma model from a microscopic to a macroscopic level scopic description. A formulation for the this of plasma temperature, model a microscopic to a macroscopic level of of and particle density profiles is Actuators sensors arecontrol designed to interact with the the current plasma density, andand to plasma and goal, particle density profiles is required. required. Actuators and scopic of from description. A formulation for thewhere burning sensors are designed to interact with plasma and to plasmalevel model from a microscopic microscopic toetaa al., macroscopic levelthe of description was derived in (Vincent 2017), plasma model from a to macroscopic level of and particle density profiles is required. Actuators and description was derived in (Vincent et al., 2017), where the sensors are designed to interact with the plasma and to measure the plasma; however, as advanced control design sensors are to interact the plasma to plasma model from a microscopic toet aisal., macroscopic level of measure thedesigned plasma; however, as with advanced control and design description was derived in (Vincent 2017), where the total energy associated to the plasma characterized, and description was derived in (Vincent et al., 2017), where the total energy associated to the plasma is characterized, and sensors to interact the plasma and to description was derived in (Vincent et al., 2017), where the measure the plasma; as advanced control design methodologies rely on onhowever, first principle models, this contribucontribumeasureare thedesigned plasma; however, as with advanced control design methodologies rely first principle models, this total energy associated to the plasma is characterized, and mass, species, momentum, energies, and entropy balance total energy associated to the plasma and is characterized, and measure the plasma; asofadvanced control design mass, species, momentum, energies, entropy balance methodologies rely on first this contribution considers the development structured controlled methodologies rely development onhowever, first principle principle models, thiscontrolled contributotal associated tousing the plasma is characterized, and tion considers the of aa models, structured mass,energy species, momentum, energies, and entropy balance equations are developed aa port-based methodology mass, species, momentum, energies, and entropy balance methodologies rely on first principle models, this contribuequations are developed using port-based methodology tion considers the development of a structured controlled model for the the the Tokamak, viewedof as as distributed multi- mass, tion considers development a structured controlled species, momentum, energies, and entropy balance model for Tokamak, viewed aa distributed multiequations are developed using a port-based port-based methodology (Duindam et al., 2009). The port-Hamiltonian formulation equations are developed using a methodology tion considers development of as a structured controlled et al., 2009). The port-Hamiltonian formulation model for the Tokamak, viewed aa distributed multiphysics system, with nonlinear couplings and evolving evolving over (Duindam model for the the Tokamak, viewed as distributed multiequations are developed using astructures port-based methodology physics system, with nonlinear couplings and over (Duindam al., 2009). The port-Hamiltonian is using Stokes–Dirac systems (Duindam et al., 2009). The port-Hamiltonian formulation is derived derived et using Stokes–Dirac structures for for formulation systems of of model for thescales Tokamak, viewed as a distributed multiphysics system, with nonlinear couplings and evolving over multiple time (Wesson and Campbell, 2004). physics system, with nonlinear couplings and evolving over (Duindam et al., 2009). The port-Hamiltonian formulation multiple time scales (Wesson and Campbell, 2004). is derived using Stokes–Dirac structures for systems of balance equations following (van der Schaft and Maschke, is derived using Stokes–Dirac structures for systems of physics system, with nonlinear couplings and evolving over balance equations following (van der Schaft and Maschke, multiple time (Wesson and Campbell, 2004). multiple time scales scales (Wesson and Campbell, 2004). ap- is derived using Stokes–Dirac for of balance equations following (vanstructures der Schaft andsystems Maschke, 2002). As the system is subjected to symmetries (vertical balance equations following (van der Schaft and Maschke, From previous works on the tokamak, a structured multiple time scales Campbell, 2004). ap- 2002). As the system is subjected to symmetries (vertical From previous works(Wesson on the and tokamak, a structured balance equations following (van der and reduction Maschke, 2002). As As the system is subjected subjected to Schaft symmetries (vertical revolution axis of a torus), we consider the model From previous works on the tokamak, a structured ap2002). the system is to symmetries (vertical proach for modeling and control of plasma profiles proved From previous works on the tokamak, a structured apaxissystem of a torus), we consider the model reduction proach for modeling and control of plasma profiles proved revolution Asfor the is subjected to symmetries (vertical revolution aa torus), we consider the reduction From previous works oncontrol the tokamak, a structured approblem theof developed port-Hamiltonian model of the the proach for and of profiles revolution axis ofdeveloped torus), we consider the model model reduction its effectiveness. To achieve achieve practical implementation of 2002). proach for modeling modeling and control of plasma plasma profiles proved proved problem foraxis the port-Hamiltonian model of its effectiveness. To aa practical implementation of revolution axis of a torus), we consider the model reduction problem for the developed port-Hamiltonian model of proach for modeling and control of plasma profiles proved Tokamak. Assuming that the momentum field is at equilibits effectiveness. To achieve a practical implementation of problem for the developed port-Hamiltonian model of the the passivity-based control of internal internal profiles in the the TCV TCV exexits effectiveness.control To achieve a practical implementation of Tokamak. Assuming that the momentum field is at equilibpassivity-based of profiles in for the dimensional developed port-Hamiltonian model of one the Tokamak. Assuming that the the momentum field is at at equilibequilibits effectiveness. To achieve a practical implementation of problem rium, the three model can be reduced to passivity-based control of internal profiles in Tokamak. Assuming that momentum field is perimental Tokamak (Vu et al., 2017), port-Hamiltonian passivity-based control ofet internal profiles in the the TCV TCV exexrium, the three dimensional model can be reduced to one perimental Tokamak (Vu al., 2017), aa port-Hamiltonian Tokamak. Assuming that the momentum field is at equilibrium, the three dimensional model can be reduced to one passivity-based control ofet internal profiles in the TCV ex- rium, dimension. The reduction achieved aa geometric perimental (Vu al., aa port-Hamiltonian the three can using be reduced to one perimental Tokamak Tokamak (Vu et al., 2017), 2017), port-Hamiltonian dimension. The dimensional reduction is ismodel achieved using geometric rium, the three dimensional model canvariable be reduced to one 1 dimension. The reduction is achieved using aa properties geometric perimental Tokamak (Vu et al., 2017), a port-Hamiltonian procedure that preserves the energy dimension. The reduction is achieved using geometric research was supported by the Belgian Network DYSCO, 1 This procedure that preserves the energy variable properties This research was supported by the Belgian Network DYSCO, 1 dimension. The reduction is achieved using a geometric procedure that preserves the energy variable properties 1 and the system structure (Vu et al., 2014). Inputs and funded by the Inter-university Attraction Poles Program, initiated This research was supported by the Belgian Network DYSCO, procedure that preserves the energy variable properties This by research was supportedAttraction by the Belgian Network DYSCO, and the system structure the (Vuenergy et al., variable 2014). Inputs and funded the Inter-university Poles Program, initiated 1 procedure that preserves properties and the system structure (Vu et al., 2014). Inputs and by the by Belgian Science Policy Attraction Office. first Program, author a FRIA funded by the Inter-university Inter-university Attraction Poles Program, initiated This research was supported by theThe Belgian Networkis DYSCO, outputs are characterized from a control perspective, such and the system structure (Vu et al., 2014). Inputs and funded the Poles initiated by the Belgian Science Policy Office. The first author is a FRIA outputs are characterized from a control perspective, such and the are system structurefrom (Vu aaetcontrol al., 2014). Inputs such and fellow by the Belgian Science Office. first author a funded by the Inter-university outputs characterized perspective, by the(F.R.S.-FNRS). Belgian Science Policy Policy Attraction Office. The ThePoles first Program, author is is initiated a FRIA FRIA outputs are characterized from control perspective, such fellow (F.R.S.-FNRS). fellow (F.R.S.-FNRS). by the(F.R.S.-FNRS). Belgian Science Policy Office. The first author is a FRIA outputs are characterized from a control perspective, such fellow
fellow (F.R.S.-FNRS). 2405-8963 © 2018 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 68 Copyright © under 2018 IFAC IFAC 68 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 68 Copyright © 2018 IFAC 68 10.1016/j.ifacol.2018.06.017 Copyright © 2018 IFAC 68
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
that control design could be achieved using passivity-based control.
69
Drive (ICRH/ICCD), the Lower Hybrid Heating and Current Drive (LHCD), and the Electron Cyclotron Heating and Current Drive (ECRH/ECCD). The pros and cons of non-inductive current drive actuators are reviewed in Wesson and Campbell (2004). The fuel supply is actuated by various methods: The gas injector releases a mixture of ions with a fixed fraction of deuterium/tritium, the pellet injector ejects deuterium and tritium at the plasma core, and the Neutral Beam Injectors is also a source of current and momentum.
The paper is structured as follows. In Section 2, the problem of nuclear fusion in a Tokamak is described, and the available actuators and sensors structure is presented. In Section 3, the first principle model is presented, and the entropy balance equation is derived using the local version of the Gibbs equation, leading to a three dimensional structured model. In Section 4, the model reduction to dimension one is presented based on a geometric considerations. Future works are discussed in Section 5.
Sensors and diagnostics Operating conditions in a Tokamak are extreme, hence measurements are sporadic and discontinuous (Wesson and Campbell, 2004). Available reliable measurements include the total plasma current, the temperature at the plasma boundary (or at the center, depending on the Tokamak), and the electron density (based on reflectometry). Despite limited available measurements, estimation of the plasma state (magnetic, temperature profiles, etc.) and estimation of various parameters (transport coefficients, etc.) via real time diagnostics can be achieved using simulation codes, for example RAPTOR, as described in (Felici et al., 2011).
2. TOKAMAK INFRASTRUCTURE AND CONTROL PROBLEM Tokamak reactors are torus-shaped plasma confinement reactors designed with the aim to sustain nuclear fusion at steady-state operation. 2.1 Fusion reaction and Lawson criteria The fusion reaction is the result of the collision of two hydrogen isotopes that generates a heavier element and energy. The reaction does not conserve mass, and the mass deficit δm is converted into average kinetic energy distributed among the reaction products in an inverse proportion to their weight. For the case of the deuteriumtritium reaction, a neutron quickly escapes the plasma while an helium element (the alpha-particle) transfers most of its energy to the plasma. The thermonuclear reaction is triggered at high level of energy density and the ignition state, where the plasma sustains itself, is reached when the inequality, i.e., the Lawson criteria, nT τe > 3.1021 m−3 keV s, (1) holds. in the above inequality, n, T , and τe denote the total plasma density, the average temperature, and the confinement time, respectively.
3. 3D BURNING PLASMA MODEL We first recall the model developed in (Vincent et al., 2017) and re-write it as a structured port-Hamiltonian model. We use notation from differential geometry (Marsden et al., 2001). For a vector field X ∈ Γ(Rn ) and a differential k-form α ∈ Λk , and using the Euclidean metric tensor, we denote the exterior derivative by dα : Λk (Rn ) ⇒ Λk+1 (Rn ), the Hodge star operator by ∗α : Λk (Rn ) ⇒ Λn−k (Rn ), the interior product iX α : Λk × Γ(Rn ) ⇒ Λk−1 (Rn ), and the flat operator X to transform a vector into a differential one-form.
Instrumentation is required to achieve the ignition point. Therefore, Tokamak reactors are equipped with sensors and actuators for real-time control of the plasma state (Wesson and Campbell, 2004).
The burning plasma is defined on a fixed domain Ω ∈ R3 with boundary ∂Ω ∈ R2 . The total energy is defined by the integral of the sum of the total electric, magnetic, kinetic, and internal energy over the spatial domain, i.e., 1 1 H= D∧E+ B∧H 2 Ω 2 Ω (2) 1 1 + ρ(∗v ) ∧ v + ρ u s, ∗ , ωk , 2 Ω ρ Ω
Actuators We can identify three classes of actuators in Tokamak reactors: The inductive and non-inductive current drive actuators and the fuel injectors. The plasma is magnetically confined by three types of coils: inner and outer poloidal, and toroidal coils, which are positioned around the plasma. The plasma is a resistive fluid and is heated by Ohmic current induced by the central coil, where the plasma itself acts as the secondary transformer. Since the plasma resistivity, denoted η, decreases as the temper−3/2 ature rises with a rate of Te (Wesson and Campbell, 2004), the Ohmic heating decreases gradually and auxiliary heating sources are required to maintain the plasma temperature profile around a desired one. The plasma is excited by electromagnetic waves at different frequency ranges with various antennae disposed on the inner vessel. Non-inductive heating and current drive actuators most commonly used are the Ion Cyclotron Heating and Current
Plasmas are governed in the electric and magnetic fields by Maxwell’s equations:
2.2 Tokamak infrastructure
where B and D ∈ Λ2 (Ω) are 2-forms denoting the magnetic and electric fluxes, respectively. The electric and magnetic fields are given by 1-forms E ∈ Λ1 (Ω) and H ∈ Λ1 (Ω), respectively. The plasma density is defined as a 0-form: ρ ∈ Λ0 (Ω), and the barycentric velocity vector is denoted by v, such that v ∈ Λ1 Ω. The specific internal energy is a 3-form u ∈ Λ3 (Ω) and is function of the specific internal entropy, the volume, and the mass fraction that are 3-forms: s, ∗ ρ1 , and ωk ∈ Λ3 (Ω), respectively. It should be noted that since the density ρ is a 0-form, the specific volume is a 3-form defined through the Hodge star operator ∗. In this system, we consider N species denoted by index k. The state of the plasma in the domain is described by balance equations for electric and magnetic fields, mass, species, momentum, energies, and entropy.
69
IFAC LHMNC 2018 70 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
∂D ∂B = −dE, and = dH − j, (3) ∂t ∂t where j ∈ Λ2 (Ω) denotes the total current density. Furthermore, Maxwell equations (3) are subject to ∇B = 0, and ∇D = 0, (4) where the quasi-neutrality of plasma is applied, i.e., the plasma charge is equal to zero. The total plasma current is the sum of the Ohmic current (Ohmic heating source), the bootstrap current (nonlinearity due to neoclassical transports), and the external current source (manipulated variable).
The internal energy density u, is a function of the entropy, the specific volume, and the mass fraction, following the Gibbs fundamental relation N 1 + µk dωk , (13) du = T ds − pd ∗ ρ k=1
from which the entropy balance equation is derived as N 1 Ds = −d ∧ Jq − ρ + σs . µ k ∧ Jk (14) Dt T k=1
The irreversible entropy source term is expressed as N 1 Jτ Jk ∧ d (µk ) − ∧ dv σs = − 2 Js ∧ dT − T T T
The total mass balance equation (or continuity equation) is given by Dρ = −d(iv ρ) − δm rfus , (5) Dt where the reaction mass deficit is denoted by δm , and the fusion reaction rate by rfus . Without loss of generalities, D· ∂· we use the material derivative: Dt = ∂t + v d· associated to the barycentric velocity v.
k=1 N
µk 1 rfus (ν k Mk + δm ωk ) (15) + ∧ (j ∧ E) − T T k=1 rfus ∗v ∧ v + Efus + δm T s + δm , T 2 and represents the generation of entropy from different phenomena: heat and material diffusion, viscosity, Joule effect, reaction, and mass deficit, respectively.
The species balance equations are given by Dωk ρ = −dJk + ωk rfus δm + ∗rfus ν k Mk , (6) Dt with Jk , ν k , and Mk , the particle flux, the stoichiometric coefficient, and the molar mass of species k, respectively.
The system balance equations are closed using the phenomenological laws for the electric and magnetic domains D = ∗ε E, and B = ∗µ H, (16) where µ and ε denote the vacuum permeability and permittivity, respectively. Total mass, species, momentum, and energy balances equations are closed using Onsager linear transport theory. A diagonal model is considered at this point of the research project, i.e., T Jq 1 v Jτ (17) j = diag(χ , κ , − , γk ) ∗ d , VE η Jk µk where χ(t, z), κ(t, z), η(t, z), and γk (t, z) denote the plasma thermal diffusivity, viscosity, resistivity, and particles transport coefficients. Here we use the relation E = −dVE , where VE ∈ Λ0 (Ω) is the electric potential. Coupled parameters in the transport matrix are not considered in the control model even if couplings exist (Boozer, 1992).
The momentum balance equation is written as Dv = −d ∗ (pI + Jτ ) − ∗(∗j ∧ ∗B) (7) Dt where p ∈ Λ0 (Ω) and Jτ ∈ Λ2 (Ω) denotes the scalar pressure and the viscous flux, respectively. The matrix identity is written as I. The last term in the right hand side of the momentum balance equation denotes the Lorentz forces (or long range forces) applied to the plasma and can be rewritten using the relationship ∗(∗j ∧ ∗B) = ij B. (8) ρ
The total energy balance equation takes the form De = −d (pI + Jτ ) ∧ v + Jq + j ∧ E + Efus ∗ rfus , (9) ρ Dt where e is the total energy per unit mass, and Efus denotes the energy released by the fusion reaction carried by the helium particle (for the deuterium/tritium reaction). The total energy mass density is given by the sum of the kinetic energy and the internal energy (∗v ) ∧ v 1 ρe = ρ + ρu s, ∗ , ωk . (10) 2 ρ
At the plasma boundary ∂Ω ∈ R2 (a torus), electric field and magnetic field are imposed by the confinement coils surrounding the plasma such that D|∂Ω = D∂ and B|∂Ω = B∂ . (18) The plasma mass density ρ is given by the sum of mass densities ρk = ρωk with boundary conditions ρk |∂Ω = 0. (19) The ideal gas law applies to plasmas (Braginskii, 1965), hence the specific (local) entropy can be defined as 3/2 T , (20) s = ln n where the total plasma density n is derived from the total mass density ρ. We can define the temperature at the plasma edge to characterize entropy of the plasma, i.e., T |∂Ω = T∂ .
The total kinetic energy balance is given by D ∗v ∧ v = −d (pI + Jτ ) ∧ v + pI ∧ dv ρ Dt 2 (11) ∗v ∧ v +Jτ ∧ dv − v ∧ (ij B) − δm rfus , 2 and the internal energy balance equation is written as Du = − d (Jq ) + j ∧ E + ∗Efus rfus + v ∧ (ij B) ρ Dt (12) ∗v ∧ v − pI ∧ dv − Jτ ∧ dv + δm rfus . 2
To develop a structured formulation of the model, we follow the work of van der Schaft and Maschke (2002) for
70
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
balance equations. The methodology is based on Stokes theorem to derive a Stokes–Dirac structure where the interconnections are explicitly written. As a result, the derived model takes the form ∂x = J ∇H(x) + g(x)u(t), (21) ∂t where x ∈ L2 (Ω) is the state, J = −J ∗ is a skewsymmetric differential operator, H(x) is the total energy function, g(x) is the input map, and u(t) defined for all t ≥ 0, is the distributed sink/source term. The interconnection operator is derived directly from the Stokes–Dirac interconnection structure. Using the mass (5), species (6), momentum (7), entropy (14) balance equations, the total energy equation (2) and the phenomenological equations (16) and (17), one can write the port-Hamiltonian system as (21).
0 −d T = d 0 Js µk Sk Sheat + σs + T T . + k 0
(28)
4. 1D BURNING PLASMA MODEL
The reduction approach for the burning plasma model relies on symmetry and one assumption on the state of the plasma, explained below.
The non-inductive current jni is the sum of the bootstrap current jbt , a nonlinear term resulting from neoclassical transport (Wesson and Campbell, 2004), and the external non-inductive current drive jext as discussed in Section 2
Axial symmetry: The toric geometry of Tokamak reactors admits an axial symmetry around the vertical revolution axis, denoted φ in Figure 1. However it should be noted that the presence of toroidal coils at discrete locations all around the plasma generates asymmetries (for example ripples effect in the electro-magnetic domain). These can be treated as perturbations in the reduced model developed below.
(23)
Using the diagonal transport model, the resistivity can be explicitly written since we have: 1 ∗ E. (24) j= η(t, z) The species balance equations are given by Dωk 0 −d µk ρ = Dt d 0 Jk Fk (ωk δm + ν k Mk ) ∗ rfus + Sk + , 0
In the proposed model, the fusion reaction term affects the momentum (26), species balance equations (25), and the thermal model (28). The model evolves on multiple time scales and can be reduced to a one dimensional system assuming that the momentum balance equation is at equilibrium (equivalently, that the dynamics of the momentum equation is fast compared to the other dynamics) (Blum, 1989).
Following (van der Schaft and Maschke, 2002), the Maxwell relations can be written as ∂ D 0 d E j + jni = − . (22) −d 0 H 0 ∂t B
jni = jbt + jext .
Ds ρ Dt FT
71
Quasi-static equilibrium: The plasma is assumed to be at a quasi-static equilibrium, in particular it is assumed that the relevant dynamics for control purposes are larger than the Alfv´en time constant (Blum, 1989). As a result, the momentum dynamics (26) is at steady-state and the pressure gradient balances the magnetic force, i.e., the Lorentz force, j × B = ∇p ⇔ ij B = d(pI), (29) is such that the constant magnetic field lines generates nested torus, called iso-surfaces. The iso-surfaces are also iso-baric, iso-thermal, and has iso-poloidal flux. Each magnetic surface is indexed in the new coordinate Φ z= , (30) πB0 where Φ and B0 denote the toroidal magnetic flux and the toroidal magnetic field at the plasma center.
(25)
where Sk represent the input source of specie k and Fk is the thermodynamic force associated to the species, i.e. Fk = dµk . The momentum and the continuity balance equation are written as Dv i B − v δ r ρ j m fus 0 −d v . δm rfus + DDt 1 = d 0 pI + Jτ ∗ ρ ρ Dt ρ (26) The resistive part in the momentum balance equation can be explicitly written using the expression of viscosity term Jτ in the transport matrix (17).
Here we use modified magnetic toric coordinates where each magnetic surface is indexed by coordinates z. As a result, the radial coordinate is function of both the magnetic surface index z and the angle θ, i.e., r = r(z, θ). The metric associated to the transformation between the toric coordinates (r,√θ, φ) and the magnetic coordinates √ (z, θ, φ) is given by g = gz gθ gφ with
The internal energy and the thermodynamic force FT = dT are linked in a port-Hamiltonian structured model as Du 0 −d T j ∧ E + Efus ∗ rrf + Sheat ρ + , Dt = d 0 Jq 0 FT (27) where Sheat is the distributed external heating source.
√
gz =
√ √
Finally, the entropy transport model is written as 71
gθ =
∂r ∂z
(31) ∂r ∂θ
2
+ r2
gφ = R0 + rcos(θ),
(32) (33)
IFAC LHMNC 2018 72 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
Fourier’s law is given by
X3 φ
√ gθ gφ J q = χ(t, z) √ dT, gz and the Fick’s law reduces to √ gθ gφ J q = γk (t, z) √ dµk . gz
δ
B
θ θ
z
z X2
R0
X1
a
0
z
Ip
Ω ∈ R3
∈ R2
Π ∈ R1
φ
where R0 denotes the large radius. The transformation coefficients can be recovered by solving the Grad–Shafranov equation and are considered for simplicity as external input parameters (Blum, 1989). As a result, the plasma is defined in one dimension with coordinate z such that z ∈ Π = [0, a] with a the small plasma radius.
β=
0
2π
√
0
gz gφ ββ dφ,
0
√
gφ αφ dφ
2π
√
gz gθ βφ dθ .
0
(36)
0
S heat (t, z) = Pheat (t)fheat (z), S k (t, z) = PMk (t)fMk (z),
(38) 2π
√
β gdθdφ.
(39)
0
and
2π 0
0
√ gθ dθ B θ = µH θ √ gz gφ
√ gz gθ dθ E φ = η(t, z)(j φ − j niφ ). √ gφ
(51) (52)
Using both reduction assumptions and following the approach proposed in (Vu et al., 2014), the three-dimensional model equations (22)(26)(25)(28) is reduced to a onedimensional model. From the quasi-static equilibrium assumption, it is assumed that the balance equation (26) is at steady state and is no longer used. The remaining balance equations (22)(25)(28) are reduced to a one-dimensional port-Hamiltonian model of the form ∂H ∂x = [J − R] + g(x)u(t), (53) ∂t ∂z
A1
2π
(50)
where fext , fheat (z), and, fMk (z) are continuous and differentiable deposition functions, and Pext (t), Pheat (t) and PMk (t) are the input total current, power and masses, respectively. We assume for simplicity, that the deposition functions are known, and that the manipulated variables are Pext (t), Pheat (t) and PMk (t).
The reduction methodology is power conservative, and the reduced variables result from their integration on the magnetic surfaces. The same analysis is applied to reduce constitutive relations. Hence the reduced Maxwell’s constitutive relations take the form
(46)
j ext = Pext (t)fext (z),
with the reduced variables
2π
Dφ (z)|z=a = Vloop ,
For the resistive diffusion equation, the reduced noninductive current denoted j ni = j bt + j ext is the sum of the reduced bootstrap current and external current. This last term is a distributed input while the plasma edge loop voltage Vloop is a boundary input. Distributed inputs j ext , S heat , and, S k result from linear combinations of all external current, heating and material sources, respectively. The distributions are defined by Gaussian deposition functions along the radial axis z. For control purposes, we assume that the control inputs are decoupled such that
Π
α=α β=
(45)
z=a
Alternatively, for α ∈ Λ0 and β ∈ Λ3 , the expression of the total energy reduces to √ α ∧ β gdzdθdφ = α ∧ β, (37) HII (α, β) = Ω
B θ (z)|z=a = Ip ,
Boundary conditions at z = 0 ensure the plasma smoothness and symmetry around the magnetic axis.
(35)
(44)
Boundary conditions of mixed type applied to the thermal and species diffusion equations, given as ∂ T ∂ρk = = 0, (47) ∂z z=0 ∂z z=0 ρk |z=a = 0 (48) = T0 . (49) T (z)
Π
2π
B θ (z)|z=0 = Dφ (z)|z=0 = 0,
where Vloop is the loop voltage applied to the plasma and Ip is the total plasma current.
Reduction methodology: Following (Vu et al., 2014), we observe that the total energy is defined as the integral of a sum of products between effort and flow variables defined in three dimensions. The reduction methodology uses the system geometry by successive integrations as depicted in Figure 1. As we are working in dimension 3, two cases of products between effort and flow variables are identified: either products of a 1-form with a 2-form (case I), or products of a 0-form with a 3-form (case II). For α ∈ Λ1 and β ∈ Λ2 , the expression of the total energy reduces to √ HI (α, β) = α ∧ β gdzdθdφ = α∧β (34) Ω
(43)
Considering the reduced one-dimensional model with toric coordinate index z ∈ [0, a], boundary conditions are given by Dirichlet conditions for the resistive diffusion equation
θ
Fig. 1. Reduction procedure to toric coordinates
with the reduced variables 2π √ α= gθ αθ dθ,
(42)
(40)
(41)
A2
72
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
73
model using axis-symmetry and steady state equilibrium of the momentum balance equation. This multi-physics portHamiltonian model consists in energy, species, entropy, and electro-magnetic balance equations. Couplings and nonlinearities resulting from the transport coefficients and the fusion reaction are the new features of the proposed model. Based on the obtained infinite dimensional model, control of Tokamak reactors with fusion reaction is to be considered using passivity-based control.
where x represents the reduced plasma state, where B, D, s, ω k denote the plasma magnetic, electric, temperature and density profiles expressed in the spatial domain Π, respectively. Note that this formulation differs from the representation (21) as dissipation is explicitly written using the phenomenological equations. For control purposes, it is more convenient to use the usual time differentiation instead of the material time derivative. A Lagrangian to Eulerian coordinates transformation is carried out when required. The total energy (10) is given by the integral H = Π H, where H is the reduced total energy density function. It consists of the sum of the electro-magnetic and internal energy densities 1 H= B ∧ H + D ∧ E + u(s, ω k ). (54) 2
REFERENCES Blum, J. (1989). Numerical Simulation and Optimal Control in Plasma Physics: With Applications to Tokamaks. Gauthier-Villars, New York. Boozer, A.H. (1992). Onsager symmetry of transport in toroidal plasmas. Physics of Fluids B, 4(9), 2845–2853. Braginskii, S. (1965). Reviews of Plasma Physics, volume 1, chapter Transport processes in a plasma, 205– 311. Consultant Bureau, New York. A. M. A. Leontovich. Duindam, V., Macchelli, A., Stramigioli, S., and Bruyninckx, H. (2009). Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach. Springer, Berlin Heidelberg. Felici, F., Sauter, O., Coda, S., Duval, B.P., Goodman, T.P., Moret, J.M., Paley, J.I., and Team, T. (2011). Real-time physics-model-based simulation of the current density profile in tokamak plasmas. Nuclear Fusion, 51(8), 083052. Marsden, J.E., Ratiu, T., and Abraham, R. (2001). Manifolds, Tensor Analysis, and Applications. SpringerVerlag, New York, 2nd edition. van der Schaft, A.J. and Maschke, B.M. (2002). Hamiltonian formulation of distributed-parameter systems with boundary energy flow. Journal of Geometry and Physics, 42(12), 166–194. Vincent, B., Hudon, N., Lef`evre, L., and Dochain, D. (2017). Burning magneto-hydrodynamics plasmas model: A port-based modelling approach. IFACPapersOnLine, 50, 13580 – 13585. Vu, N.M.T., Lef`evre, L., and Maschke, B. (2016). A structured control model for the thermo-magnetohydrodynamics of plasmas in tokamaks. Mathematical and Computer Modelling of Dynamical Systems, 22(3), 181–206. Vu, N.M.T., Lef`evre, L., Nouailletas, R., and Br´emond, S. (2014). Structure preserving reduction for thermomagneto plasma control model. In Proceedings of the 21st International Symposium on Mathematical Theory of Networks and Systems (MTNS), 716–723. Groningen, Netherlands. Vu, N.M.T., Nouailletas, R., Maljaars, E., Felici, F., and Sauter, O. (2017). Plasma internal profile control using IDA-PBC: Application to TCV. Fusion Engineering and Design. doi:10.1016/j.fusengdes.2017.02.074. In Press. Wesson, J. and Campbell, D.J. (2004). Tokamaks. Clarendon Press, New York.
Maxwell relations applied to the plasma in Tokamak reduced in one dimension gives two decoupled sub-models: a poloidal model (Dφ , B θ , H θ , E φ , and J φ ), and a toroidal model (Dθ , B φ , H φ , E θ , and J θ ). The well known resistive diffusion equation for the poloidal magnetic flux (Wesson and Campbell, 2004), when the electric equilibrium is achieved (Blum, 1989), results from the poloidal sub model, and leads to −1 ∂ Dφ ∂ Eφ 0 1 Rem (η) 0 = − 0 0 Hθ ∂t B θ ∂z − 1 0 (55) j ext − , 0
where Dφ , B θ , E φ and H θ denote the poloidal electric and toroidal magnetic fluxes, and the poloidal electric and toroidal magnetic field, respectively.
−1 (η) is derived from the plasma The dissipative element Rem resistivity η(t, z). Thus the resistivity operator is given by 1 Rem (η) = η(t, z). (56) A2
The external current source jext is the sum of the noninduced current that is an input source, and the bootstrap current a result of neoclassical transport arising in plasmas, leading to the structures √ ∂ω k ∂ 0 −1 µk 0 0 ρ g = − −1 ∂t 1 0 (γ ) 0 R J ∂z k k k 0 r (δ ω + ν k Mk ) + S k + fus m k , 0 (57) and √ ∂s ∂ 0 −1 T gρ ∂t = ∂z 1 0 J s FT (58) S heat + k µk S k . + σ s + T 0 5. CONCLUSION We presented a structured three-dimensional burning plasma model and derived a one-dimensional reduced 73
ScienceDirect
IFAC PapersOnLine 51-3 (2018) 68–73
Port-Hamiltonian Port-Hamiltonian modeling modeling and and reduction reduction Port-Hamiltonian modeling and reduction of a burning plasma system Port-Hamiltonian modeling and reduction of a burning plasma system of a burning plasma system of a burning plasma ∗∗ system ∗∗∗ Benjamin Vincent ∗,∗∗ ∗,∗∗ Trang Vu ∗∗ Nicolas Hudon ∗∗∗
Benjamin Vincent ∗,∗∗ Trang Vu ∗∗ Nicolas Hudon ∗∗∗ ∗∗ ∗ ∗,∗∗ ∗∗∗ Benjamin Vincent Trang Vu Nicolas e Dochain ∗∗ Denis ∗ BenjaminLaurent VincentLef` Trang Vu ∗∗ Nicolas Hudon Hudon Laurent Lef` evre vre Denis Dochain ∗,∗∗ ∗∗ ∗∗∗ ∗∗ ∗ ∗∗ Denis ∗ BenjaminLaurent Vincent Trang Vu Nicolas Hudon Lef` e vre Dochain Laurent Lef` evre ∗∗ Denis Dochain ∗ ∗ Laurent Lef` evre Denis Dochain ee Catholique a ∗ ICTEAM, Universit´ Catholique de de Louvain, Louvain, Bˆ Bˆ atiment timent Euler, Euler, av. av. ∗ ICTEAM, Universit´ ∗ ICTEAM, Universit´ eB-1348 Catholique de Louvain, Louvain, Bˆ Bˆ atiment timent Euler, av. av. Lemaˆ ıtre 4-6, Louvain-La-Neuve, Belgium. ICTEAM, Universit´ e Catholique de a Euler, Lemaˆ ıtre 4-6, B-1348 Louvain-La-Neuve, Belgium. ∗ ICTEAM, Universit´ Catholique de Louvain, Bˆ atiment Euler, av. Lemaˆ ıtre 4-6, 4-6,eB-1348 B-1348 Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) Lemaˆ ıtre Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) ∗∗Lemaˆ ıtre 4-6, B-1348 Alpes, Louvain-La-Neuve, Belgium. (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) ∗∗ Universit´ (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) Universit´ee Grenoble Grenoble Alpes, LCIS, LCIS, F-26902, F-26902, France. France. ∗∗ ∗∗ Universit´ (e-mail:{benjamin.vincent,denis.dochain}@uclouvain.be) e Grenoble Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) Universit´ e Grenoble Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) ∗∗∗ ∗∗ Universit´ Alpes, LCIS, F-26902, France. (e-mail:laurent.lefevre,[email protected]) Chemical Engineering, Queen’s University, ∗∗∗ Department eofGrenoble (e-mail:laurent.lefevre,[email protected]) Department of Chemical Engineering, Queen’s University, ∗∗∗ ∗∗∗ (e-mail:laurent.lefevre,[email protected]) Department of Chemical Engineering, Queen’s University, University, Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) Department of Chemical Engineering, Queen’s Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) ∗∗∗ Department Chemical Queen’s University, Kingston, ON, 3N6 Canada. (e-mail: Kingston, ON, K7L K7L of 3N6 Canada.Engineering, (e-mail: [email protected]) [email protected]) Kingston, ON, K7L 3N6 Canada. (e-mail: [email protected]) Abstract: Abstract: In In this this contribution, contribution, we we develop develop aa structured structured port-Hamiltonian port-Hamiltonian model model for for aa class class Abstract: In this we structured port-Hamiltonian model class of parameter systems exploit obtained for Abstract: In distributed this contribution, contribution, we develop develop structured port-Hamiltonian model for for class of irreversible irreversible distributed parameter systems aaand and exploit the the obtained formulation formulation for aamodel model Abstract: In distributed this contribution, we develop aand structured port-Hamiltonian modelis for class of irreversible distributed parameter systems and exploit the obtained formulation for amodel model reduction from dimension 3 to dimension 1 in space. The proposed methodology motivated of irreversible parameter systems exploit the obtained formulation for reduction from dimension 3 to dimension 1 in space. The proposed methodology is motivated of irreversible parameter systems and exploit the obtained formulation for model reduction fromdistributed dimension 3 to to dimension dimension 1 in in space. The proposed methodology is motivated by the control of burning plasma profiles in Tokamak reactors. The burning plasma is viewed reduction from dimension 3 1 space. The proposed methodology is motivated by the control burning plasma profiles in Tokamak reactors. The methodology burning plasma is viewed reduction from of dimension 3built to dimension 1 in space. The proposed ismomentum, motivated by the control of burning plasma profiles in Tokamak reactors. The burning plasma is as system on Maxwell equations and total species, by a themulti-physics control of burning plasma profiles in Tokamak reactors. Themass, burning plasma is viewed viewed as a multi-physics system built on Maxwell equations and total mass, species, momentum, by the control of burning plasma profiles in Tokamak reactors. The burning plasma is viewed as a multi-physics system built on Maxwell equations and total mass, species, momentum, energies, and entropy balance equations. Moreover, the system presents nonlinear couplings, as a multi-physics system builtequations. on Maxwell equations totalpresents mass, species, Moreover, the and system nonlinearmomentum, couplings, energies, and entropy balance as a multi-physics system built on Maxwell equations totalpresents mass, species, energies, and entropy balance equations. Moreover, the and system presents nonlinearmomentum, couplings, especially through transport coefficients, and evolves over time scales. energies, and entropy balance equations. Moreover, the system nonlinear couplings, especially through transport and its its dynamic dynamic evolves over multiple multiple time scales. energies, and entropy balancecoefficients, equations. Moreover, the system presents nonlinear couplings, especially through transport coefficients, and Joule its dynamic dynamic evolves overforces, multiple time scales. The main couplings considered here are the effect, the Lorentz and the fusion especially through transport coefficients, and its evolves over multiple time scales. The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion especially through transport coefficients, and its dynamic evolves overrely multiple time scales. The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction on the use of the The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction rely on the use of the The main couplings considered here are the Joule effect, the Lorentz forces, and the fusion reaction kinetics. The port-based modeling formulation and reduction rely on the use of the Gibbs relation, Onsager linear transport theory, Stokes–Dirac structures, and energy preserving reaction kinetics. The port-based modeling formulation and structures, reduction rely on the preserving use of the Gibbs relation, Onsager linear transport theory, Stokes–Dirac and energy reaction kinetics. The port-based modeling formulation and structures, reduction rely the preserving use of the Gibbs relation, relation, Onsager linear transport transport theory, Stokes–Dirac structures, and on energy preserving geometric reduction. Gibbs Onsager linear theory, Stokes–Dirac and energy geometric reduction. Gibbs relation, Onsager linear transport theory, Stokes–Dirac structures, and energy preserving geometric reduction. geometric reduction. © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. geometric reduction. Keywords: Keywords: Port-Hamiltonian Port-Hamiltonian systems, systems, Tokamak Tokamak reactors, reactors, Modeling, Modeling, Model Model reduction, reduction, Keywords: Port-Hamiltonian systems, Tokamak reactors, Modeling, Model Irreversible thermodynamics Keywords: Port-Hamiltonian systems, Tokamak reactors, Modeling, Model reduction, reduction, Irreversible thermodynamics Keywords: Irreversible thermodynamics IrreversiblePort-Hamiltonian thermodynamics systems, Tokamak reactors, Modeling, Model reduction, Irreversible thermodynamics 1. INTRODUCTION INTRODUCTION model 1. model has has been been proposed proposed (Vu (Vu et et al., al., 2016), 2016), and and then then 1. model has has been proposed proposed (Vu et al., 2016), and then then reduced. It should be noted that the model considered in 1. INTRODUCTION INTRODUCTION model been (Vu et al., 2016), and reduced. It should be noted that the model considered in 1.promising INTRODUCTION model has been did proposed (Vu al., 2016), and then reduced. It 2016) should be not noted that et the model considered in (Vu et al., consider the material field, i.e., reduced. It should be noted that the model considered in Nuclear fusion is a source of renewable energy et al., 2016) did not consider the material field, i.e., Nuclear fusion is a promising source of renewable energy (Vu reduced. It should be balance noted that the model considered in (Vu et al., 2016) did not consider the material field, i.e., the mass and species equations, and by extension Nuclear fusion is a promising source of renewable energy (Vu et al., 2016) did not consider the material field, i.e., able to solve present and future energy production issues. Nuclear fusion is a promising source of renewable energy mass and species balance equations, and by field, extension able to solve present and future energy production issues. the (Vu et al., 2016) did not consider the material i.e., mass and species balance equations, and by extension Nuclear fusion is a promising source of renewable energy the fusion reaction, were not considered. able to solve present and future energy production issues. specieswere balance equations, and by extension Tokamak reactors areand onefuture of the theenergy investigated technology able to solve present production issues. the mass fusionand reaction, not considered. Tokamak reactors are one of investigated technology specieswere balance equations, and by extension the mass fusionand reaction, were not considered. considered. able to solve present production issues. Tokamak reactors are one of the investigated technology fusion reaction, not to achieve fusion reaction; however, many technical chal- the Tokamak reactors areand onefuture of theenergy investigated technology to achieve fusion reaction; however, many technical chalIn paper, we aa structured the fusion reaction, were not considered.model Tokamak reactors are one before of the energy investigated technology In this this paper, we propose propose structured model for for burning burning to achieve fusion reaction; however, many technical challenges remain to be solved can be produced to achieve fusion reaction; however, many technical challenges remain to be solved before energy can be produced plasma In this paper, we propose aa structured model for burning in Tokamak reactors, i.e., aa structured model takIn this paper, we propose structured model for burning to achieve fusion reaction; however, many technical chalplasma in Tokamak reactors, i.e., structured model taklenges remain to be solved before energy can be produced using this class of reactors. One of the remaining problem lenges remain to be solved before energy can be produced this paper, wethe propose a reaction. structured model for burning using this class of reactors. One of the remaining problem In plasma inaccount Tokamak reactors, i.e., aa structured structured model taking into fusion This reaction occurs plasma in Tokamak reactors, i.e., model taklenges remain to be solved before energy can be produced ing into account the fusion reaction. This reaction occurs using this class of reactors. One of the remaining problem concerns the stabilization of One plasma around desiredproblem profiles plasma using thisthe class of reactors. of the remaining inaccount Tokamak reactors, i.e., a structured model takconcerns stabilization of plasma around desired profiles ingainto into account the fusion reaction. Thiscannot reaction occurs at microscopic level of description and be added ing the fusion reaction. This reaction occurs using thisthe class of reactors. of remaining ainto microscopic levelfusion of description and cannot be occurs added concerns stabilization of plasma around desired profiles such that the fusion reaction is the sustained. To problem achieve concerns the stabilization of One plasma around desired profiles at ing account the reaction. This reaction such that the fusion reaction is sustained. To achieve at a microscopic level of description and cannot be added directly to principle model at a a microscopic of description and cannot added concerns stabilization of plasma aroundcurrent desired profiles at directly to aa first firstlevel principle model developed developed at be a macromacrosuch that the fusion reaction is To achieve this goal, control of plasma plasma temperature, density, such goal, thatthe the fusion reaction is sustained. sustained. To density, achieve at a microscopic level of description and cannot be added this control of temperature, current directly to aaoffirst principle model developed at a macroscopic level description. A formulation for the burning directly to first principle model developed at a macrosuch that the fusion reaction is sustained. To achieve scopic level of description. A formulation for the burning this control of density, and particle density profilestemperature, is required. required. current Actuators and directly this goal, goal, control of plasma plasma temperature, current density, to aof principle developed macroand particle density profiles is Actuators and scopic level level offirst description. Amodel formulation for at theaburning burning plasma model from a microscopic to a macroscopic level scopic description. A formulation for the this of plasma temperature, model a microscopic to a macroscopic level of of and particle density profiles is Actuators sensors arecontrol designed to interact with the the current plasma density, andand to plasma and goal, particle density profiles is required. required. Actuators and scopic of from description. A formulation for thewhere burning sensors are designed to interact with plasma and to plasmalevel model from a microscopic microscopic toetaa al., macroscopic levelthe of description was derived in (Vincent 2017), plasma model from a to macroscopic level of and particle density profiles is required. Actuators and description was derived in (Vincent et al., 2017), where the sensors are designed to interact with the plasma and to measure the plasma; however, as advanced control design sensors are to interact the plasma to plasma model from a microscopic toet aisal., macroscopic level of measure thedesigned plasma; however, as with advanced control and design description was derived in (Vincent 2017), where the total energy associated to the plasma characterized, and description was derived in (Vincent et al., 2017), where the total energy associated to the plasma is characterized, and sensors to interact the plasma and to description was derived in (Vincent et al., 2017), where the measure the plasma; as advanced control design methodologies rely on onhowever, first principle models, this contribucontribumeasureare thedesigned plasma; however, as with advanced control design methodologies rely first principle models, this total energy associated to the plasma is characterized, and mass, species, momentum, energies, and entropy balance total energy associated to the plasma and is characterized, and measure the plasma; asofadvanced control design mass, species, momentum, energies, entropy balance methodologies rely on first this contribution considers the development structured controlled methodologies rely development onhowever, first principle principle models, thiscontrolled contributotal associated tousing the plasma is characterized, and tion considers the of aa models, structured mass,energy species, momentum, energies, and entropy balance equations are developed aa port-based methodology mass, species, momentum, energies, and entropy balance methodologies rely on first principle models, this contribuequations are developed using port-based methodology tion considers the development of a structured controlled model for the the the Tokamak, viewedof as as distributed multi- mass, tion considers development a structured controlled species, momentum, energies, and entropy balance model for Tokamak, viewed aa distributed multiequations are developed using a port-based port-based methodology (Duindam et al., 2009). The port-Hamiltonian formulation equations are developed using a methodology tion considers development of as a structured controlled et al., 2009). The port-Hamiltonian formulation model for the Tokamak, viewed aa distributed multiphysics system, with nonlinear couplings and evolving evolving over (Duindam model for the the Tokamak, viewed as distributed multiequations are developed using astructures port-based methodology physics system, with nonlinear couplings and over (Duindam al., 2009). The port-Hamiltonian is using Stokes–Dirac systems (Duindam et al., 2009). The port-Hamiltonian formulation is derived derived et using Stokes–Dirac structures for for formulation systems of of model for thescales Tokamak, viewed as a distributed multiphysics system, with nonlinear couplings and evolving over multiple time (Wesson and Campbell, 2004). physics system, with nonlinear couplings and evolving over (Duindam et al., 2009). The port-Hamiltonian formulation multiple time scales (Wesson and Campbell, 2004). is derived using Stokes–Dirac structures for systems of balance equations following (van der Schaft and Maschke, is derived using Stokes–Dirac structures for systems of physics system, with nonlinear couplings and evolving over balance equations following (van der Schaft and Maschke, multiple time (Wesson and Campbell, 2004). multiple time scales scales (Wesson and Campbell, 2004). ap- is derived using Stokes–Dirac for of balance equations following (vanstructures der Schaft andsystems Maschke, 2002). As the system is subjected to symmetries (vertical balance equations following (van der Schaft and Maschke, From previous works on the tokamak, a structured multiple time scales Campbell, 2004). ap- 2002). As the system is subjected to symmetries (vertical From previous works(Wesson on the and tokamak, a structured balance equations following (van der and reduction Maschke, 2002). As As the system is subjected subjected to Schaft symmetries (vertical revolution axis of a torus), we consider the model From previous works on the tokamak, a structured ap2002). the system is to symmetries (vertical proach for modeling and control of plasma profiles proved From previous works on the tokamak, a structured apaxissystem of a torus), we consider the model reduction proach for modeling and control of plasma profiles proved revolution Asfor the is subjected to symmetries (vertical revolution aa torus), we consider the reduction From previous works oncontrol the tokamak, a structured approblem theof developed port-Hamiltonian model of the the proach for and of profiles revolution axis ofdeveloped torus), we consider the model model reduction its effectiveness. To achieve achieve practical implementation of 2002). proach for modeling modeling and control of plasma plasma profiles proved proved problem foraxis the port-Hamiltonian model of its effectiveness. To aa practical implementation of revolution axis of a torus), we consider the model reduction problem for the developed port-Hamiltonian model of proach for modeling and control of plasma profiles proved Tokamak. Assuming that the momentum field is at equilibits effectiveness. To achieve a practical implementation of problem for the developed port-Hamiltonian model of the the passivity-based control of internal internal profiles in the the TCV TCV exexits effectiveness.control To achieve a practical implementation of Tokamak. Assuming that the momentum field is at equilibpassivity-based of profiles in for the dimensional developed port-Hamiltonian model of one the Tokamak. Assuming that the the momentum field is at at equilibequilibits effectiveness. To achieve a practical implementation of problem rium, the three model can be reduced to passivity-based control of internal profiles in Tokamak. Assuming that momentum field is perimental Tokamak (Vu et al., 2017), port-Hamiltonian passivity-based control ofet internal profiles in the the TCV TCV exexrium, the three dimensional model can be reduced to one perimental Tokamak (Vu al., 2017), aa port-Hamiltonian Tokamak. Assuming that the momentum field is at equilibrium, the three dimensional model can be reduced to one passivity-based control ofet internal profiles in the TCV ex- rium, dimension. The reduction achieved aa geometric perimental (Vu al., aa port-Hamiltonian the three can using be reduced to one perimental Tokamak Tokamak (Vu et al., 2017), 2017), port-Hamiltonian dimension. The dimensional reduction is ismodel achieved using geometric rium, the three dimensional model canvariable be reduced to one 1 dimension. The reduction is achieved using aa properties geometric perimental Tokamak (Vu et al., 2017), a port-Hamiltonian procedure that preserves the energy dimension. The reduction is achieved using geometric research was supported by the Belgian Network DYSCO, 1 This procedure that preserves the energy variable properties This research was supported by the Belgian Network DYSCO, 1 dimension. The reduction is achieved using a geometric procedure that preserves the energy variable properties 1 and the system structure (Vu et al., 2014). Inputs and funded by the Inter-university Attraction Poles Program, initiated This research was supported by the Belgian Network DYSCO, procedure that preserves the energy variable properties This by research was supportedAttraction by the Belgian Network DYSCO, and the system structure the (Vuenergy et al., variable 2014). Inputs and funded the Inter-university Poles Program, initiated 1 procedure that preserves properties and the system structure (Vu et al., 2014). Inputs and by the by Belgian Science Policy Attraction Office. first Program, author a FRIA funded by the Inter-university Inter-university Attraction Poles Program, initiated This research was supported by theThe Belgian Networkis DYSCO, outputs are characterized from a control perspective, such and the system structure (Vu et al., 2014). Inputs and funded the Poles initiated by the Belgian Science Policy Office. The first author is a FRIA outputs are characterized from a control perspective, such and the are system structurefrom (Vu aaetcontrol al., 2014). Inputs such and fellow by the Belgian Science Office. first author a funded by the Inter-university outputs characterized perspective, by the(F.R.S.-FNRS). Belgian Science Policy Policy Attraction Office. The ThePoles first Program, author is is initiated a FRIA FRIA outputs are characterized from control perspective, such fellow (F.R.S.-FNRS). fellow (F.R.S.-FNRS). by the(F.R.S.-FNRS). Belgian Science Policy Office. The first author is a FRIA outputs are characterized from a control perspective, such fellow
fellow (F.R.S.-FNRS). 2405-8963 © 2018 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 68 Copyright © under 2018 IFAC IFAC 68 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 68 Copyright © 2018 IFAC 68 10.1016/j.ifacol.2018.06.017 Copyright © 2018 IFAC 68
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
that control design could be achieved using passivity-based control.
69
Drive (ICRH/ICCD), the Lower Hybrid Heating and Current Drive (LHCD), and the Electron Cyclotron Heating and Current Drive (ECRH/ECCD). The pros and cons of non-inductive current drive actuators are reviewed in Wesson and Campbell (2004). The fuel supply is actuated by various methods: The gas injector releases a mixture of ions with a fixed fraction of deuterium/tritium, the pellet injector ejects deuterium and tritium at the plasma core, and the Neutral Beam Injectors is also a source of current and momentum.
The paper is structured as follows. In Section 2, the problem of nuclear fusion in a Tokamak is described, and the available actuators and sensors structure is presented. In Section 3, the first principle model is presented, and the entropy balance equation is derived using the local version of the Gibbs equation, leading to a three dimensional structured model. In Section 4, the model reduction to dimension one is presented based on a geometric considerations. Future works are discussed in Section 5.
Sensors and diagnostics Operating conditions in a Tokamak are extreme, hence measurements are sporadic and discontinuous (Wesson and Campbell, 2004). Available reliable measurements include the total plasma current, the temperature at the plasma boundary (or at the center, depending on the Tokamak), and the electron density (based on reflectometry). Despite limited available measurements, estimation of the plasma state (magnetic, temperature profiles, etc.) and estimation of various parameters (transport coefficients, etc.) via real time diagnostics can be achieved using simulation codes, for example RAPTOR, as described in (Felici et al., 2011).
2. TOKAMAK INFRASTRUCTURE AND CONTROL PROBLEM Tokamak reactors are torus-shaped plasma confinement reactors designed with the aim to sustain nuclear fusion at steady-state operation. 2.1 Fusion reaction and Lawson criteria The fusion reaction is the result of the collision of two hydrogen isotopes that generates a heavier element and energy. The reaction does not conserve mass, and the mass deficit δm is converted into average kinetic energy distributed among the reaction products in an inverse proportion to their weight. For the case of the deuteriumtritium reaction, a neutron quickly escapes the plasma while an helium element (the alpha-particle) transfers most of its energy to the plasma. The thermonuclear reaction is triggered at high level of energy density and the ignition state, where the plasma sustains itself, is reached when the inequality, i.e., the Lawson criteria, nT τe > 3.1021 m−3 keV s, (1) holds. in the above inequality, n, T , and τe denote the total plasma density, the average temperature, and the confinement time, respectively.
3. 3D BURNING PLASMA MODEL We first recall the model developed in (Vincent et al., 2017) and re-write it as a structured port-Hamiltonian model. We use notation from differential geometry (Marsden et al., 2001). For a vector field X ∈ Γ(Rn ) and a differential k-form α ∈ Λk , and using the Euclidean metric tensor, we denote the exterior derivative by dα : Λk (Rn ) ⇒ Λk+1 (Rn ), the Hodge star operator by ∗α : Λk (Rn ) ⇒ Λn−k (Rn ), the interior product iX α : Λk × Γ(Rn ) ⇒ Λk−1 (Rn ), and the flat operator X to transform a vector into a differential one-form.
Instrumentation is required to achieve the ignition point. Therefore, Tokamak reactors are equipped with sensors and actuators for real-time control of the plasma state (Wesson and Campbell, 2004).
The burning plasma is defined on a fixed domain Ω ∈ R3 with boundary ∂Ω ∈ R2 . The total energy is defined by the integral of the sum of the total electric, magnetic, kinetic, and internal energy over the spatial domain, i.e., 1 1 H= D∧E+ B∧H 2 Ω 2 Ω (2) 1 1 + ρ(∗v ) ∧ v + ρ u s, ∗ , ωk , 2 Ω ρ Ω
Actuators We can identify three classes of actuators in Tokamak reactors: The inductive and non-inductive current drive actuators and the fuel injectors. The plasma is magnetically confined by three types of coils: inner and outer poloidal, and toroidal coils, which are positioned around the plasma. The plasma is a resistive fluid and is heated by Ohmic current induced by the central coil, where the plasma itself acts as the secondary transformer. Since the plasma resistivity, denoted η, decreases as the temper−3/2 ature rises with a rate of Te (Wesson and Campbell, 2004), the Ohmic heating decreases gradually and auxiliary heating sources are required to maintain the plasma temperature profile around a desired one. The plasma is excited by electromagnetic waves at different frequency ranges with various antennae disposed on the inner vessel. Non-inductive heating and current drive actuators most commonly used are the Ion Cyclotron Heating and Current
Plasmas are governed in the electric and magnetic fields by Maxwell’s equations:
2.2 Tokamak infrastructure
where B and D ∈ Λ2 (Ω) are 2-forms denoting the magnetic and electric fluxes, respectively. The electric and magnetic fields are given by 1-forms E ∈ Λ1 (Ω) and H ∈ Λ1 (Ω), respectively. The plasma density is defined as a 0-form: ρ ∈ Λ0 (Ω), and the barycentric velocity vector is denoted by v, such that v ∈ Λ1 Ω. The specific internal energy is a 3-form u ∈ Λ3 (Ω) and is function of the specific internal entropy, the volume, and the mass fraction that are 3-forms: s, ∗ ρ1 , and ωk ∈ Λ3 (Ω), respectively. It should be noted that since the density ρ is a 0-form, the specific volume is a 3-form defined through the Hodge star operator ∗. In this system, we consider N species denoted by index k. The state of the plasma in the domain is described by balance equations for electric and magnetic fields, mass, species, momentum, energies, and entropy.
69
IFAC LHMNC 2018 70 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
∂D ∂B = −dE, and = dH − j, (3) ∂t ∂t where j ∈ Λ2 (Ω) denotes the total current density. Furthermore, Maxwell equations (3) are subject to ∇B = 0, and ∇D = 0, (4) where the quasi-neutrality of plasma is applied, i.e., the plasma charge is equal to zero. The total plasma current is the sum of the Ohmic current (Ohmic heating source), the bootstrap current (nonlinearity due to neoclassical transports), and the external current source (manipulated variable).
The internal energy density u, is a function of the entropy, the specific volume, and the mass fraction, following the Gibbs fundamental relation N 1 + µk dωk , (13) du = T ds − pd ∗ ρ k=1
from which the entropy balance equation is derived as N 1 Ds = −d ∧ Jq − ρ + σs . µ k ∧ Jk (14) Dt T k=1
The irreversible entropy source term is expressed as N 1 Jτ Jk ∧ d (µk ) − ∧ dv σs = − 2 Js ∧ dT − T T T
The total mass balance equation (or continuity equation) is given by Dρ = −d(iv ρ) − δm rfus , (5) Dt where the reaction mass deficit is denoted by δm , and the fusion reaction rate by rfus . Without loss of generalities, D· ∂· we use the material derivative: Dt = ∂t + v d· associated to the barycentric velocity v.
k=1 N
µk 1 rfus (ν k Mk + δm ωk ) (15) + ∧ (j ∧ E) − T T k=1 rfus ∗v ∧ v + Efus + δm T s + δm , T 2 and represents the generation of entropy from different phenomena: heat and material diffusion, viscosity, Joule effect, reaction, and mass deficit, respectively.
The species balance equations are given by Dωk ρ = −dJk + ωk rfus δm + ∗rfus ν k Mk , (6) Dt with Jk , ν k , and Mk , the particle flux, the stoichiometric coefficient, and the molar mass of species k, respectively.
The system balance equations are closed using the phenomenological laws for the electric and magnetic domains D = ∗ε E, and B = ∗µ H, (16) where µ and ε denote the vacuum permeability and permittivity, respectively. Total mass, species, momentum, and energy balances equations are closed using Onsager linear transport theory. A diagonal model is considered at this point of the research project, i.e., T Jq 1 v Jτ (17) j = diag(χ , κ , − , γk ) ∗ d , VE η Jk µk where χ(t, z), κ(t, z), η(t, z), and γk (t, z) denote the plasma thermal diffusivity, viscosity, resistivity, and particles transport coefficients. Here we use the relation E = −dVE , where VE ∈ Λ0 (Ω) is the electric potential. Coupled parameters in the transport matrix are not considered in the control model even if couplings exist (Boozer, 1992).
The momentum balance equation is written as Dv = −d ∗ (pI + Jτ ) − ∗(∗j ∧ ∗B) (7) Dt where p ∈ Λ0 (Ω) and Jτ ∈ Λ2 (Ω) denotes the scalar pressure and the viscous flux, respectively. The matrix identity is written as I. The last term in the right hand side of the momentum balance equation denotes the Lorentz forces (or long range forces) applied to the plasma and can be rewritten using the relationship ∗(∗j ∧ ∗B) = ij B. (8) ρ
The total energy balance equation takes the form De = −d (pI + Jτ ) ∧ v + Jq + j ∧ E + Efus ∗ rfus , (9) ρ Dt where e is the total energy per unit mass, and Efus denotes the energy released by the fusion reaction carried by the helium particle (for the deuterium/tritium reaction). The total energy mass density is given by the sum of the kinetic energy and the internal energy (∗v ) ∧ v 1 ρe = ρ + ρu s, ∗ , ωk . (10) 2 ρ
At the plasma boundary ∂Ω ∈ R2 (a torus), electric field and magnetic field are imposed by the confinement coils surrounding the plasma such that D|∂Ω = D∂ and B|∂Ω = B∂ . (18) The plasma mass density ρ is given by the sum of mass densities ρk = ρωk with boundary conditions ρk |∂Ω = 0. (19) The ideal gas law applies to plasmas (Braginskii, 1965), hence the specific (local) entropy can be defined as 3/2 T , (20) s = ln n where the total plasma density n is derived from the total mass density ρ. We can define the temperature at the plasma edge to characterize entropy of the plasma, i.e., T |∂Ω = T∂ .
The total kinetic energy balance is given by D ∗v ∧ v = −d (pI + Jτ ) ∧ v + pI ∧ dv ρ Dt 2 (11) ∗v ∧ v +Jτ ∧ dv − v ∧ (ij B) − δm rfus , 2 and the internal energy balance equation is written as Du = − d (Jq ) + j ∧ E + ∗Efus rfus + v ∧ (ij B) ρ Dt (12) ∗v ∧ v − pI ∧ dv − Jτ ∧ dv + δm rfus . 2
To develop a structured formulation of the model, we follow the work of van der Schaft and Maschke (2002) for
70
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
balance equations. The methodology is based on Stokes theorem to derive a Stokes–Dirac structure where the interconnections are explicitly written. As a result, the derived model takes the form ∂x = J ∇H(x) + g(x)u(t), (21) ∂t where x ∈ L2 (Ω) is the state, J = −J ∗ is a skewsymmetric differential operator, H(x) is the total energy function, g(x) is the input map, and u(t) defined for all t ≥ 0, is the distributed sink/source term. The interconnection operator is derived directly from the Stokes–Dirac interconnection structure. Using the mass (5), species (6), momentum (7), entropy (14) balance equations, the total energy equation (2) and the phenomenological equations (16) and (17), one can write the port-Hamiltonian system as (21).
0 −d T = d 0 Js µk Sk Sheat + σs + T T . + k 0
(28)
4. 1D BURNING PLASMA MODEL
The reduction approach for the burning plasma model relies on symmetry and one assumption on the state of the plasma, explained below.
The non-inductive current jni is the sum of the bootstrap current jbt , a nonlinear term resulting from neoclassical transport (Wesson and Campbell, 2004), and the external non-inductive current drive jext as discussed in Section 2
Axial symmetry: The toric geometry of Tokamak reactors admits an axial symmetry around the vertical revolution axis, denoted φ in Figure 1. However it should be noted that the presence of toroidal coils at discrete locations all around the plasma generates asymmetries (for example ripples effect in the electro-magnetic domain). These can be treated as perturbations in the reduced model developed below.
(23)
Using the diagonal transport model, the resistivity can be explicitly written since we have: 1 ∗ E. (24) j= η(t, z) The species balance equations are given by Dωk 0 −d µk ρ = Dt d 0 Jk Fk (ωk δm + ν k Mk ) ∗ rfus + Sk + , 0
In the proposed model, the fusion reaction term affects the momentum (26), species balance equations (25), and the thermal model (28). The model evolves on multiple time scales and can be reduced to a one dimensional system assuming that the momentum balance equation is at equilibrium (equivalently, that the dynamics of the momentum equation is fast compared to the other dynamics) (Blum, 1989).
Following (van der Schaft and Maschke, 2002), the Maxwell relations can be written as ∂ D 0 d E j + jni = − . (22) −d 0 H 0 ∂t B
jni = jbt + jext .
Ds ρ Dt FT
71
Quasi-static equilibrium: The plasma is assumed to be at a quasi-static equilibrium, in particular it is assumed that the relevant dynamics for control purposes are larger than the Alfv´en time constant (Blum, 1989). As a result, the momentum dynamics (26) is at steady-state and the pressure gradient balances the magnetic force, i.e., the Lorentz force, j × B = ∇p ⇔ ij B = d(pI), (29) is such that the constant magnetic field lines generates nested torus, called iso-surfaces. The iso-surfaces are also iso-baric, iso-thermal, and has iso-poloidal flux. Each magnetic surface is indexed in the new coordinate Φ z= , (30) πB0 where Φ and B0 denote the toroidal magnetic flux and the toroidal magnetic field at the plasma center.
(25)
where Sk represent the input source of specie k and Fk is the thermodynamic force associated to the species, i.e. Fk = dµk . The momentum and the continuity balance equation are written as Dv i B − v δ r ρ j m fus 0 −d v . δm rfus + DDt 1 = d 0 pI + Jτ ∗ ρ ρ Dt ρ (26) The resistive part in the momentum balance equation can be explicitly written using the expression of viscosity term Jτ in the transport matrix (17).
Here we use modified magnetic toric coordinates where each magnetic surface is indexed by coordinates z. As a result, the radial coordinate is function of both the magnetic surface index z and the angle θ, i.e., r = r(z, θ). The metric associated to the transformation between the toric coordinates (r,√θ, φ) and the magnetic coordinates √ (z, θ, φ) is given by g = gz gθ gφ with
The internal energy and the thermodynamic force FT = dT are linked in a port-Hamiltonian structured model as Du 0 −d T j ∧ E + Efus ∗ rrf + Sheat ρ + , Dt = d 0 Jq 0 FT (27) where Sheat is the distributed external heating source.
√
gz =
√ √
Finally, the entropy transport model is written as 71
gθ =
∂r ∂z
(31) ∂r ∂θ
2
+ r2
gφ = R0 + rcos(θ),
(32) (33)
IFAC LHMNC 2018 72 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
Fourier’s law is given by
X3 φ
√ gθ gφ J q = χ(t, z) √ dT, gz and the Fick’s law reduces to √ gθ gφ J q = γk (t, z) √ dµk . gz
δ
B
θ θ
z
z X2
R0
X1
a
0
z
Ip
Ω ∈ R3
∈ R2
Π ∈ R1
φ
where R0 denotes the large radius. The transformation coefficients can be recovered by solving the Grad–Shafranov equation and are considered for simplicity as external input parameters (Blum, 1989). As a result, the plasma is defined in one dimension with coordinate z such that z ∈ Π = [0, a] with a the small plasma radius.
β=
0
2π
√
0
gz gφ ββ dφ,
0
√
gφ αφ dφ
2π
√
gz gθ βφ dθ .
0
(36)
0
S heat (t, z) = Pheat (t)fheat (z), S k (t, z) = PMk (t)fMk (z),
(38) 2π
√
β gdθdφ.
(39)
0
and
2π 0
0
√ gθ dθ B θ = µH θ √ gz gφ
√ gz gθ dθ E φ = η(t, z)(j φ − j niφ ). √ gφ
(51) (52)
Using both reduction assumptions and following the approach proposed in (Vu et al., 2014), the three-dimensional model equations (22)(26)(25)(28) is reduced to a onedimensional model. From the quasi-static equilibrium assumption, it is assumed that the balance equation (26) is at steady state and is no longer used. The remaining balance equations (22)(25)(28) are reduced to a one-dimensional port-Hamiltonian model of the form ∂H ∂x = [J − R] + g(x)u(t), (53) ∂t ∂z
A1
2π
(50)
where fext , fheat (z), and, fMk (z) are continuous and differentiable deposition functions, and Pext (t), Pheat (t) and PMk (t) are the input total current, power and masses, respectively. We assume for simplicity, that the deposition functions are known, and that the manipulated variables are Pext (t), Pheat (t) and PMk (t).
The reduction methodology is power conservative, and the reduced variables result from their integration on the magnetic surfaces. The same analysis is applied to reduce constitutive relations. Hence the reduced Maxwell’s constitutive relations take the form
(46)
j ext = Pext (t)fext (z),
with the reduced variables
2π
Dφ (z)|z=a = Vloop ,
For the resistive diffusion equation, the reduced noninductive current denoted j ni = j bt + j ext is the sum of the reduced bootstrap current and external current. This last term is a distributed input while the plasma edge loop voltage Vloop is a boundary input. Distributed inputs j ext , S heat , and, S k result from linear combinations of all external current, heating and material sources, respectively. The distributions are defined by Gaussian deposition functions along the radial axis z. For control purposes, we assume that the control inputs are decoupled such that
Π
α=α β=
(45)
z=a
Alternatively, for α ∈ Λ0 and β ∈ Λ3 , the expression of the total energy reduces to √ α ∧ β gdzdθdφ = α ∧ β, (37) HII (α, β) = Ω
B θ (z)|z=a = Ip ,
Boundary conditions at z = 0 ensure the plasma smoothness and symmetry around the magnetic axis.
(35)
(44)
Boundary conditions of mixed type applied to the thermal and species diffusion equations, given as ∂ T ∂ρk = = 0, (47) ∂z z=0 ∂z z=0 ρk |z=a = 0 (48) = T0 . (49) T (z)
Π
2π
B θ (z)|z=0 = Dφ (z)|z=0 = 0,
where Vloop is the loop voltage applied to the plasma and Ip is the total plasma current.
Reduction methodology: Following (Vu et al., 2014), we observe that the total energy is defined as the integral of a sum of products between effort and flow variables defined in three dimensions. The reduction methodology uses the system geometry by successive integrations as depicted in Figure 1. As we are working in dimension 3, two cases of products between effort and flow variables are identified: either products of a 1-form with a 2-form (case I), or products of a 0-form with a 3-form (case II). For α ∈ Λ1 and β ∈ Λ2 , the expression of the total energy reduces to √ HI (α, β) = α ∧ β gdzdθdφ = α∧β (34) Ω
(43)
Considering the reduced one-dimensional model with toric coordinate index z ∈ [0, a], boundary conditions are given by Dirichlet conditions for the resistive diffusion equation
θ
Fig. 1. Reduction procedure to toric coordinates
with the reduced variables 2π √ α= gθ αθ dθ,
(42)
(40)
(41)
A2
72
IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018
Benjamin Vincent et al. / IFAC PapersOnLine 51-3 (2018) 68–73
73
model using axis-symmetry and steady state equilibrium of the momentum balance equation. This multi-physics portHamiltonian model consists in energy, species, entropy, and electro-magnetic balance equations. Couplings and nonlinearities resulting from the transport coefficients and the fusion reaction are the new features of the proposed model. Based on the obtained infinite dimensional model, control of Tokamak reactors with fusion reaction is to be considered using passivity-based control.
where x represents the reduced plasma state, where B, D, s, ω k denote the plasma magnetic, electric, temperature and density profiles expressed in the spatial domain Π, respectively. Note that this formulation differs from the representation (21) as dissipation is explicitly written using the phenomenological equations. For control purposes, it is more convenient to use the usual time differentiation instead of the material time derivative. A Lagrangian to Eulerian coordinates transformation is carried out when required. The total energy (10) is given by the integral H = Π H, where H is the reduced total energy density function. It consists of the sum of the electro-magnetic and internal energy densities 1 H= B ∧ H + D ∧ E + u(s, ω k ). (54) 2
REFERENCES Blum, J. (1989). Numerical Simulation and Optimal Control in Plasma Physics: With Applications to Tokamaks. Gauthier-Villars, New York. Boozer, A.H. (1992). Onsager symmetry of transport in toroidal plasmas. Physics of Fluids B, 4(9), 2845–2853. Braginskii, S. (1965). Reviews of Plasma Physics, volume 1, chapter Transport processes in a plasma, 205– 311. Consultant Bureau, New York. A. M. A. Leontovich. Duindam, V., Macchelli, A., Stramigioli, S., and Bruyninckx, H. (2009). Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach. Springer, Berlin Heidelberg. Felici, F., Sauter, O., Coda, S., Duval, B.P., Goodman, T.P., Moret, J.M., Paley, J.I., and Team, T. (2011). Real-time physics-model-based simulation of the current density profile in tokamak plasmas. Nuclear Fusion, 51(8), 083052. Marsden, J.E., Ratiu, T., and Abraham, R. (2001). Manifolds, Tensor Analysis, and Applications. SpringerVerlag, New York, 2nd edition. van der Schaft, A.J. and Maschke, B.M. (2002). Hamiltonian formulation of distributed-parameter systems with boundary energy flow. Journal of Geometry and Physics, 42(12), 166–194. Vincent, B., Hudon, N., Lef`evre, L., and Dochain, D. (2017). Burning magneto-hydrodynamics plasmas model: A port-based modelling approach. IFACPapersOnLine, 50, 13580 – 13585. Vu, N.M.T., Lef`evre, L., and Maschke, B. (2016). A structured control model for the thermo-magnetohydrodynamics of plasmas in tokamaks. Mathematical and Computer Modelling of Dynamical Systems, 22(3), 181–206. Vu, N.M.T., Lef`evre, L., Nouailletas, R., and Br´emond, S. (2014). Structure preserving reduction for thermomagneto plasma control model. In Proceedings of the 21st International Symposium on Mathematical Theory of Networks and Systems (MTNS), 716–723. Groningen, Netherlands. Vu, N.M.T., Nouailletas, R., Maljaars, E., Felici, F., and Sauter, O. (2017). Plasma internal profile control using IDA-PBC: Application to TCV. Fusion Engineering and Design. doi:10.1016/j.fusengdes.2017.02.074. In Press. Wesson, J. and Campbell, D.J. (2004). Tokamaks. Clarendon Press, New York.
Maxwell relations applied to the plasma in Tokamak reduced in one dimension gives two decoupled sub-models: a poloidal model (Dφ , B θ , H θ , E φ , and J φ ), and a toroidal model (Dθ , B φ , H φ , E θ , and J θ ). The well known resistive diffusion equation for the poloidal magnetic flux (Wesson and Campbell, 2004), when the electric equilibrium is achieved (Blum, 1989), results from the poloidal sub model, and leads to −1 ∂ Dφ ∂ Eφ 0 1 Rem (η) 0 = − 0 0 Hθ ∂t B θ ∂z − 1 0 (55) j ext − , 0
where Dφ , B θ , E φ and H θ denote the poloidal electric and toroidal magnetic fluxes, and the poloidal electric and toroidal magnetic field, respectively.
−1 (η) is derived from the plasma The dissipative element Rem resistivity η(t, z). Thus the resistivity operator is given by 1 Rem (η) = η(t, z). (56) A2
The external current source jext is the sum of the noninduced current that is an input source, and the bootstrap current a result of neoclassical transport arising in plasmas, leading to the structures √ ∂ω k ∂ 0 −1 µk 0 0 ρ g = − −1 ∂t 1 0 (γ ) 0 R J ∂z k k k 0 r (δ ω + ν k Mk ) + S k + fus m k , 0 (57) and √ ∂s ∂ 0 −1 T gρ ∂t = ∂z 1 0 J s FT (58) S heat + k µk S k . + σ s + T 0 5. CONCLUSION We presented a structured three-dimensional burning plasma model and derived a one-dimensional reduced 73