- Email: [email protected]
A PROCEDURE FOR THE COMPUTATION OF MAGNETIC FIELD, TEMPERATURE AND FORCE DISTRIBUTIONS IN THE DESIGN OF TOKAMAK TOROIDAL COILS
A PROCEDURE FOR THE COMPUTATION OF MAGNETIC FIELD, TEMPERATURE AND FORCE DISTRIBUTIONS IN THE DESIGN OF TOKAMAK TOROIDAL COILS
0
0
*
0
j
P. Molfino , G. Molinari , F. Rosatelli , A. Viviani
°
Università di Genova, 1st. di Elettrotecnica, viale Causa 13, 16145 GENOVA, GÈ
*
NIRA, via dei Pescatori 35, 16128 GENOVA, GÈ
Abstract The paper presents a formulation of the problem of coupled electric, magnetic and thermal fields in a TOKAMAK toroidal coil, oriented to give satisfactory results with reduced computer resources. This goal is achieved by the introduction of an equivalent two-dimensional coil, whose solution is related to the solution of the original three-di mensional problem. A package allowing the solution of the resulting set of partial diffe rential
equations
is also presented, its features are described and some results are
displayed and discussed.
1.
Introduction The design of toroidal field magnets for TOKAMAK machines requires the computation
of coupled electric, magnetic and thermal fields under non linear time-varying condi tions. The field unknowns interact each other essentially by modifying the material pro perties, course,
such the
as
the
presence
electric
and
of a thermal
thermal field
conductivity
is connected
and
the specific
to power
losses
due
heat. Of to eddy
currents in the conductors. The machine geometry is not symmetric, owing to the spaces between turns and to the presence of windows and insulators. However, a complete three-dimensional formula tion is not practical, since it is extremely heavy as to computation requirements. As a consequence, an equivalent two-dimensional formulation has been developed, by Which a solution oftsufficient precision can be obtained with a reasonable computation effort. This formulation computes magnetic field, current density and temperature confi gurations by introducing a fictitious current density to match together all the field quantities in an axisymmetric conductor geometry.
1107
In the following, the equivalent formulation is reported and the tools used to sol ve the field problem are described. Finally, some of the obtained results are displayed and discussed.
2.
Equivalent formulation The geometric parameters relevant to the toroidal coil of a Tokamak machine are
defined in fig. 1. A cylindrical coordinate system r, θ, ζ is used. An r-z-plane section of the coil
is shown in fig. la, whereas an r-θ-plane
section is shown in fig. lb.
The thickness of an elementary turn has been assumed as not constant to obtain maximum generality. · The equivalent
turn, taking place of an elementary turn in the two-dimensional
model, is represented in fig. lc.
b) r-6-plane section of the coil.
c) r-0-plane section of the equi]/ aient coil.
Let H and J represent the values of the magnetic field and the current density in the wire of fig. lb. It is assumed that H presents rotation symmetry. This can be accep ted as a reasonable approximation, even if the discrete nature of the winding deforms the magnetic lines of force, which are not circles centered on the magnet axis. This assumption makes this field compatible with the axisymmetric geometry of fig. le. As a
1108
consequence, the equivalent turn is the seat of a fictitious current density J*, related to H by the Maxwell equation: curl H
(1)
J can be expressed as a function of J* by computing the line integral of H on a circle centered on the coil axis (trace B in fig. la). Using eq. 1 and Stokes formula applied to a surface with rotation symmetry (trace s in fig. la), it is 21 f-H · —d l -fΓ2πΓΒ I ' Jo J A
ΓΒ.
J* · nds r άθ = 2π -N
J* · ri r ds
fB _ J · ή ds r άθ = 2ΝΘf J-e 1 A
J · n r ds
(2)
where N is the number of turns and
(3) e' = arctg <.* +
*)***
Since eq. (2) holds for any surface, it is possible to express J as a function of 7*, for small values of φ, in the form: j ,
π
ΝΘ1
j* * 21 - Ntg<|>
r
e+r
(4)
J* - k —£— J* e +r
An expression of Ohm law, including Hall e f f e c t ,
is
E+hy0JxH = pJ where p is
the
(5) electric
resistivity
and h
is the Hall constant. As a consequence,
expressing the electric field as function of H: curl E ■ - 1Ì-.-TT" = curl (pJ - h u n J x H) = curl u
0 ot
—T— (p J * - h u J * x H ) e+r o '
(6)
an equation for H can be obtained: curl
k — — (p curl H + h μυ0 H x curl H) e+r
3H
Λ
(7)
Eq.(7) reduce to a scalar equation since H has a component only in the Θ direction. The equation for the temperature T is: div(KgradT) -Cp-|^---p J 2
(8)
where κ is the thermal conductivity and c is the specific heat. P All parameters, like p, h, κ and c , can be expressed as function of temperature. P Eqs. (7) and
(8) form a set of nonlinear time-varying equations in two unknowns. The
boundary conditions are
1109
H = 0
on the external surface of the coil
(9)
— NI H = —
on the internal surface of the coil
(10)
on the surface of the coil
(11)
2ΤΓΓ
3T —— = 0 3n
where I is the current impressed in the coil. Condition (11) is an approximation valid for short time-intervals. For longer intervals, heat exchanges are to be accounted for. Owing to condition eq.(ll), the temperature equation eq.(8) presents an infinite number of solutions, differing from each other as to the reference level. The solution becomes unique if a thermal energy balance is made at each computation step, in order to evaluate the correct reference level. Since this balance makes the computation heavier, a common hypothesis is to assume that the heating process is locally adiabatic when the transient to be studied is very short, for instance under pulse conditions. This assump tion is also oriented to a more safe design of the coil, since the temperature level reached
is
higher.
For
this
reason,it
has
been
adopted
for the examples
reported
in this paper. 3.
Numerical solution The solution of the field problem defined above has been obtained by means of a
package,
developed
at
the
Electrical
Engineering
Department
of
the
University
of
Genoa/1/. The package allows the solution of boundary value problems defined by sets of partial differential equations and by general boundary and interface conditions under non linear time-varying conditions, provided that after time discretization the equa tions can be reduced to the elliptic type. The methods used are a generalized finite difference method allowing irregular grids /2/, or the finite element method. In the finite difference version the package allows to define the equations and the boundary and interface conditions as input data, together with material properties, user
defined
functions
of
the
unknown
solution
(such
as
forces, losses and other
quantities of interest for the user) and other functions for various purposes. This is achieved by an input processor, working like an interpreter, which is able to read and decode strings of characters written in a FORTRAN - Like Language, developed by some of the authors. In this language, the user can input data not only by numerical values, but even by general FORTRAN expressions. These expressions can contain numbers, parameters, space and time coordinates, unknown field quantities and their derivatives, user defined functions (material properties and other functions). The user can define parameters by special input cards, giving their names and their values. Once parameters have been defi
lilo
ned, they can be used in the expressions as FORTRAN variables. The user can also define material properties, giving their names and their expressions for each subregion into which the definition^ field of the problem is subdivided. Besides, the user can define functions of interest by giving names and expressions. Also user defined functions can be accepted in the FORTRAN expression after their definition. At
run
time,
the
input
is
interpreted
and
the
expressions
are analized.
If
all the data involved in an expression have a numerical value previously assigned, the expression is evaluated and the result is stored where necessary. On the contrary, if the expression cannot be evaluated, since it contains unknown values, i.e. coordinates or unknown field quantities, it is simplified as far as possible and then it is stored in a coded form as an expression, whose value will be computed by the program after the needed values have been obtained. This feature allows to define the problem to be solved by input data, without re compiling the package. The implementation of the axisymmetric model previously described has been done in a complete form, including Hall effect (which could give a small pertur bation
on
the
field
configuration
in
the
copper
conductors
owing
to the
magnetic
flux density of about 10 T and to the enormous current density), since this can be done without any difficulty. Other features of the package are segmented structure (35 segments, about 500 rou tines) and package virtual memory at FORTRAN level, to allow portability on machines without virtual memory; debugging facilities included in the package, with continuous tests about the consistency of the partial results; small· memory requirements. Run-time features are automatic grid generation, solution of the set of algebraic equations by both iterative or direct methods, time discretization by both implicit or explicit procedures, graphic output facilities. The package is under continuous develop ment: for instance, grid optimization procedures are under implementation /3/.
4.
Solution examples The equivalent formulation presented here has been applied to compute the elec
tric, magnetic and thermal field of a fusion machine. The machine has the dimensional values of FT, the Frascati T0KAMAK unit. The computation has been performed assuming that the applied voltage induces in the coil a current varying in time as given in fig. 2, under the condition of locally adiabatic heating of the conductors. Figs. 3 to 7 show equifunction curves relevant to magnetic field, H, absolute value of current density, J, temperature, T, radial and ver tical components of force density, f and f . r z
1111
All results refer to the time instant corresponding to the begin ning
of the
current
plateau. The
number of nodes used has been of 95 and the discretization val has
time inter
been of .1 s. The computa
tion time has been of 250 CPU se conds on a CDC
170/835
machine to
compute the whole transient for the time interval of fig. 2.
5.
Fig. 2. Time evolution of the applied current.
Conclusions The package described in this paper has been applied to solve the field problems of
TOKAMAK
machines
in some cases, using the formulation presented here. The experience
obtained till now has shown good stability of the results with respect to both the number of nodes
and
the
important
feature,
time since
discretization it
means
interval
that
used
results
in the computation. This is an
of sufficient
approximation
can be
obtained with reasonable computer resources, in spite of the hundreds of field solutions needed to describe the time evolution of nonlinear systems, especially if their geometry is more complex than that of the example shown in this paper. A very useful tool is also the possibility to change rapidly the formulation and the geometry of the problem by exploiting the features of the input processor, which allows to input data as functions of parameters and to obtain output data in the form nearest to the design requirements.
6.
References
/l/
P. Molfino, G. Molinari, A. Viviani : "A user oriented package for the solution of general
field
problems
under
time
varying
conditions",
IEEE
Transactions
on
Magnetics, 1982, vol. MAG-18, n. 2, pp. 638-643. /2/
G. Molinari, M.R. Podestà, G. Sciutto, A. Viviani: "Finite difference method with irregular grid
and transformed discretization metric", IEEE PES Winter Meeting,
1978, paper A 78 288/3. /3/
G. Molinari, A. Viviani: "Grid iteration method for finite element grid optimiza tion", Proc. Nat. Congress on Pressure Vessels and Piping, S. Francisco, 1979, vol. PVP-38, pp. 49-59.
1112
Fig. 3. Equif unction curves of the magnetic field.
Fig. 4. Equifunction curves of the current den sity.
Fig. 5. Equif unction curves of the temperature.
Fig. 6. Equifunction curves of the r component of the force density.
Fig. 7. Equifunction curves of the z component of the force density.
1113
0
0
*
0
j
P. Molfino , G. Molinari , F. Rosatelli , A. Viviani
°
Università di Genova, 1st. di Elettrotecnica, viale Causa 13, 16145 GENOVA, GÈ
*
NIRA, via dei Pescatori 35, 16128 GENOVA, GÈ
Abstract The paper presents a formulation of the problem of coupled electric, magnetic and thermal fields in a TOKAMAK toroidal coil, oriented to give satisfactory results with reduced computer resources. This goal is achieved by the introduction of an equivalent two-dimensional coil, whose solution is related to the solution of the original three-di mensional problem. A package allowing the solution of the resulting set of partial diffe rential
equations
is also presented, its features are described and some results are
displayed and discussed.
1.
Introduction The design of toroidal field magnets for TOKAMAK machines requires the computation
of coupled electric, magnetic and thermal fields under non linear time-varying condi tions. The field unknowns interact each other essentially by modifying the material pro perties, course,
such the
as
the
presence
electric
and
of a thermal
thermal field
conductivity
is connected
and
the specific
to power
losses
due
heat. Of to eddy
currents in the conductors. The machine geometry is not symmetric, owing to the spaces between turns and to the presence of windows and insulators. However, a complete three-dimensional formula tion is not practical, since it is extremely heavy as to computation requirements. As a consequence, an equivalent two-dimensional formulation has been developed, by Which a solution oftsufficient precision can be obtained with a reasonable computation effort. This formulation computes magnetic field, current density and temperature confi gurations by introducing a fictitious current density to match together all the field quantities in an axisymmetric conductor geometry.
1107
In the following, the equivalent formulation is reported and the tools used to sol ve the field problem are described. Finally, some of the obtained results are displayed and discussed.
2.
Equivalent formulation The geometric parameters relevant to the toroidal coil of a Tokamak machine are
defined in fig. 1. A cylindrical coordinate system r, θ, ζ is used. An r-z-plane section of the coil
is shown in fig. la, whereas an r-θ-plane
section is shown in fig. lb.
The thickness of an elementary turn has been assumed as not constant to obtain maximum generality. · The equivalent
turn, taking place of an elementary turn in the two-dimensional
model, is represented in fig. lc.
b) r-6-plane section of the coil.
c) r-0-plane section of the equi]/ aient coil.
Let H and J represent the values of the magnetic field and the current density in the wire of fig. lb. It is assumed that H presents rotation symmetry. This can be accep ted as a reasonable approximation, even if the discrete nature of the winding deforms the magnetic lines of force, which are not circles centered on the magnet axis. This assumption makes this field compatible with the axisymmetric geometry of fig. le. As a
1108
consequence, the equivalent turn is the seat of a fictitious current density J*, related to H by the Maxwell equation: curl H
(1)
J can be expressed as a function of J* by computing the line integral of H on a circle centered on the coil axis (trace B in fig. la). Using eq. 1 and Stokes formula applied to a surface with rotation symmetry (trace s in fig. la), it is 21 f-H · —d l -fΓ2πΓΒ I ' Jo J A
ΓΒ.
J* · nds r άθ = 2π -N
J* · ri r ds
fB _ J · ή ds r άθ = 2ΝΘf J-e 1 A
J · n r ds
(2)
where N is the number of turns and
(3) e' = arctg <.* +
*)***
Since eq. (2) holds for any surface, it is possible to express J as a function of 7*, for small values of φ, in the form: j ,
π
ΝΘ1
j* * 21 - Ntg<|>
r
e+r
(4)
J* - k —£— J* e +r
An expression of Ohm law, including Hall e f f e c t ,
is
E+hy0JxH = pJ where p is
the
(5) electric
resistivity
and h
is the Hall constant. As a consequence,
expressing the electric field as function of H: curl E ■ - 1Ì-.-TT" = curl (pJ - h u n J x H) = curl u
0 ot
—T— (p J * - h u J * x H ) e+r o '
(6)
an equation for H can be obtained: curl
k — — (p curl H + h μυ0 H x curl H) e+r
3H
Λ
(7)
Eq.(7) reduce to a scalar equation since H has a component only in the Θ direction. The equation for the temperature T is: div(KgradT) -Cp-|^---p J 2
(8)
where κ is the thermal conductivity and c is the specific heat. P All parameters, like p, h, κ and c , can be expressed as function of temperature. P Eqs. (7) and
(8) form a set of nonlinear time-varying equations in two unknowns. The
boundary conditions are
1109
H = 0
on the external surface of the coil
(9)
— NI H = —
on the internal surface of the coil
(10)
on the surface of the coil
(11)
2ΤΓΓ
3T —— = 0 3n
where I is the current impressed in the coil. Condition (11) is an approximation valid for short time-intervals. For longer intervals, heat exchanges are to be accounted for. Owing to condition eq.(ll), the temperature equation eq.(8) presents an infinite number of solutions, differing from each other as to the reference level. The solution becomes unique if a thermal energy balance is made at each computation step, in order to evaluate the correct reference level. Since this balance makes the computation heavier, a common hypothesis is to assume that the heating process is locally adiabatic when the transient to be studied is very short, for instance under pulse conditions. This assump tion is also oriented to a more safe design of the coil, since the temperature level reached
is
higher.
For
this
reason,it
has
been
adopted
for the examples
reported
in this paper. 3.
Numerical solution The solution of the field problem defined above has been obtained by means of a
package,
developed
at
the
Electrical
Engineering
Department
of
the
University
of
Genoa/1/. The package allows the solution of boundary value problems defined by sets of partial differential equations and by general boundary and interface conditions under non linear time-varying conditions, provided that after time discretization the equa tions can be reduced to the elliptic type. The methods used are a generalized finite difference method allowing irregular grids /2/, or the finite element method. In the finite difference version the package allows to define the equations and the boundary and interface conditions as input data, together with material properties, user
defined
functions
of
the
unknown
solution
(such
as
forces, losses and other
quantities of interest for the user) and other functions for various purposes. This is achieved by an input processor, working like an interpreter, which is able to read and decode strings of characters written in a FORTRAN - Like Language, developed by some of the authors. In this language, the user can input data not only by numerical values, but even by general FORTRAN expressions. These expressions can contain numbers, parameters, space and time coordinates, unknown field quantities and their derivatives, user defined functions (material properties and other functions). The user can define parameters by special input cards, giving their names and their values. Once parameters have been defi
lilo
ned, they can be used in the expressions as FORTRAN variables. The user can also define material properties, giving their names and their expressions for each subregion into which the definition^ field of the problem is subdivided. Besides, the user can define functions of interest by giving names and expressions. Also user defined functions can be accepted in the FORTRAN expression after their definition. At
run
time,
the
input
is
interpreted
and
the
expressions
are analized.
If
all the data involved in an expression have a numerical value previously assigned, the expression is evaluated and the result is stored where necessary. On the contrary, if the expression cannot be evaluated, since it contains unknown values, i.e. coordinates or unknown field quantities, it is simplified as far as possible and then it is stored in a coded form as an expression, whose value will be computed by the program after the needed values have been obtained. This feature allows to define the problem to be solved by input data, without re compiling the package. The implementation of the axisymmetric model previously described has been done in a complete form, including Hall effect (which could give a small pertur bation
on
the
field
configuration
in
the
copper
conductors
owing
to the
magnetic
flux density of about 10 T and to the enormous current density), since this can be done without any difficulty. Other features of the package are segmented structure (35 segments, about 500 rou tines) and package virtual memory at FORTRAN level, to allow portability on machines without virtual memory; debugging facilities included in the package, with continuous tests about the consistency of the partial results; small· memory requirements. Run-time features are automatic grid generation, solution of the set of algebraic equations by both iterative or direct methods, time discretization by both implicit or explicit procedures, graphic output facilities. The package is under continuous develop ment: for instance, grid optimization procedures are under implementation /3/.
4.
Solution examples The equivalent formulation presented here has been applied to compute the elec
tric, magnetic and thermal field of a fusion machine. The machine has the dimensional values of FT, the Frascati T0KAMAK unit. The computation has been performed assuming that the applied voltage induces in the coil a current varying in time as given in fig. 2, under the condition of locally adiabatic heating of the conductors. Figs. 3 to 7 show equifunction curves relevant to magnetic field, H, absolute value of current density, J, temperature, T, radial and ver tical components of force density, f and f . r z
1111
All results refer to the time instant corresponding to the begin ning
of the
current
plateau. The
number of nodes used has been of 95 and the discretization val has
time inter
been of .1 s. The computa
tion time has been of 250 CPU se conds on a CDC
170/835
machine to
compute the whole transient for the time interval of fig. 2.
5.
Fig. 2. Time evolution of the applied current.
Conclusions The package described in this paper has been applied to solve the field problems of
TOKAMAK
machines
in some cases, using the formulation presented here. The experience
obtained till now has shown good stability of the results with respect to both the number of nodes
and
the
important
feature,
time since
discretization it
means
interval
that
used
results
in the computation. This is an
of sufficient
approximation
can be
obtained with reasonable computer resources, in spite of the hundreds of field solutions needed to describe the time evolution of nonlinear systems, especially if their geometry is more complex than that of the example shown in this paper. A very useful tool is also the possibility to change rapidly the formulation and the geometry of the problem by exploiting the features of the input processor, which allows to input data as functions of parameters and to obtain output data in the form nearest to the design requirements.
6.
References
/l/
P. Molfino, G. Molinari, A. Viviani : "A user oriented package for the solution of general
field
problems
under
time
varying
conditions",
IEEE
Transactions
on
Magnetics, 1982, vol. MAG-18, n. 2, pp. 638-643. /2/
G. Molinari, M.R. Podestà, G. Sciutto, A. Viviani: "Finite difference method with irregular grid
and transformed discretization metric", IEEE PES Winter Meeting,
1978, paper A 78 288/3. /3/
G. Molinari, A. Viviani: "Grid iteration method for finite element grid optimiza tion", Proc. Nat. Congress on Pressure Vessels and Piping, S. Francisco, 1979, vol. PVP-38, pp. 49-59.
1112
Fig. 3. Equif unction curves of the magnetic field.
Fig. 4. Equifunction curves of the current den sity.
Fig. 5. Equif unction curves of the temperature.
Fig. 6. Equifunction curves of the r component of the force density.
Fig. 7. Equifunction curves of the z component of the force density.
1113