Nuclear Engineering and Design 59 (1980) 401 4 1 0 © North-Holland Publishing Company
FINITE ELEMENT ANALYSIS OF EXPERIMENT ON DYNAMIC BEHAVIOR OF CYLINDER DUE TO ELECTROMAGNETIC FORCE KenzoMIYA, Kosei HARA and Yoneho TABATA Nuclear Engineering Research Laboratory, Faculty of Engineering, University of Tokyo, Tokai.mura, Ibaraki-Prefecture, Japan Received 2 January 1980
An experiment on the dynamic behavior of a cylinder resulting from a transient electromagnetic force was conducted. This experiment can be considered associated with the dynamic deformation of a first wall or a vecuum vessel of a controlled thermonuclear tokamak type reactor. The dynamic strain of the cylinder was measured under a transient strong magnetic field using a non-inductive strain gage and a shielding room. A finite element method was applied to the analysis of the experiment by solving Maxwell's equations and the equation of motion of the cylinder simultaneously. The numerical results are in good agreement with the experimental results.
1. Introduction
tions of the components due to these forces. On the other hand a great deal of work has already been done with respect to the magnetic forming of metals [3].M0st of those works are associated with the plastic deformation of the metals, which is usually difficult to analyse without referring to numerical methods such as the finite element method or the finite difference method. Some simple problems were investigated experimentally and analytically, and the results show fairly good agreement. For example, Forrestal and Overmier [4] examined the response of a circular elastic ring to a cosine distributed impulse over half of the ring circumference experimentally and analytically. Comparison of the measured membrane strain with theory was very good. Walling and Forrestal [5] also examined the response of a circular 6061T6 aluminum ring deformed in the plastic range and got very good agreement between experiment and prediction. Lal and Hillier [6] gave a theoretical analysis of electromagnetic forces exerted on a tube and showed that time-varying tube inductance, circuit resistance and irreversible plastic work are significant parameters of the electromagnetic forming machine. The model studied by them is similar to one adoped in the present paper. The method o f analysis adopted here is the finite element method. Application of the method releases
Electromagnetic forces have been applied to metal forming for many years. Recently a new problem associated with electromagnetic forces has appeared in tile engineering field of magnetic confinement fusion reactors [1 ]. A realization of the magnetic confinement fusion power plant will require a very strong magnetic field of about ten teslas at the locations where plasmas of high temperature a n d densityishall be confined. Sources of the transient high magnetic field are the plasma current itself, carrying a few mega amperes flowing in a toroidal direction, and a current in poloidal coils for Ohmic heating of the plasma in an initial stage. Eddy currents induced in metal components by these transient electrical sources interact with the strong toroidal field as well as the transient field produced by these sources; resulting in very large electromagnetic forces in the metal components such as the first wall, vacuum vessel, blanket structure and the toroidal field coils [2]. From the viewpoint of making an optimal design of a fusion power plant a great effort should be made to minimize the electromagnetic forces on structural components. Studies done up to now in the engineering field of fusion reactors almost totally exclude such a sort of investigation of transient electromagnetic forces and the deforma401
402
K. Miya et aL / Dynamic behavior of cylinder due to electromagnet force
us from some assumptions needed when trying to solve even a considerably simpler problem analytically. In other words, the finite element method can calculate a more detailed and accurate representation of the behavior of a specimen subject to transient electromagnetic forces and assures the possibility of analyse a large deformation or the plastic work of metals. It is usually rather difficult to measure a dynamic deformation under a time varying strong magnetic field when using a common strain gage. The reason is that it is very hard to eliminate completely noises arising from various sources such as induced current in measuring devices, difference of ground levels and electromagnetic waves. In this study some contrivances were made. Most parts of the measuring system were contained in a shield room with an insulated transformer. And non-inductive strain gages were used to minimize the induced current.
Jo = current density supplied by external source, u = displacement of conductor. Eq. (4) corresponds to Faraday's law of induction, namely the electric field intensity E is produced by the changing magnetic field. Eqs. (1), (2) and (4) can be simplified if the vector potential A is introduced. Eq. (4) suggests mathematically the following relation between the vector potential and the magnetic induction B = rot A ,
(6)
because eq. (6) automatically satisfies eq. (2). A substitution of eq. (6) into eq. (4) gives, E = - bA/Ot.
(7)
If eqs. (1), (3) and (6) are combined, it is seen that the vector potential satisfies the following equation rot rot A = # d .
(8)
2.2. Finite element formulation o f axi-symmetric
2. Finite element analysis of dynamic electromagnetic force
field
As is well-known, an electromagnetic field is governed by Maxwell equations supplemented by Ohm's law. In the case of the electromagnetic force discussed in this paper, the displacement current can be neglecte neglected. Thus Maxwell equations are,
The cylindrical coordinates (r, O, z)are adopted here as the geometry of the experimental device is axi-symmetric, as will be shown later. In such a case the current density J is in the tangential direction only, and is a function of r and z. This requires that the vector potential should also have only a tangential component, denoted by A if eq. (8), or the law of Biot-Savart, is considered. And thus,
rot H = J ,
(1)
J
div B = 0 ,
(2)
Jo = (0, io, 0),
s = ~e,
(3)
A = (0, A, 0).
rot E = - a B / ~ t ,
(4)
From eqs. (5) and (8) we get
2.1. Maxwell equations
and Ohm's law is; J = o [ e + 3u/3t
× B] +Jo,
where H = magnetic field intensity, B = magnetic induction, E = electric field intensity, = permeability, o = electrical conductivity, J = current density,
(5)
=(0,/,0),
-r p ~rLr~r
[L
] at
--(rA) p~z 2
"(~-~XB)o]+]°}r:O'
(9)
The term (au/at) X B)o in eq. (9) is the tangential component of the current density induced by the motion of the conductor in the magnetic field of B. Tiffs equation determines the distributed induced current density of the cylinder shown schematically in fig. 1.
K. Miya et al. / Dynamic behavior of cylinder due to electromagnet force Copper Foil Shield ~lilil
403
Cylindricol Shell radius = 24. 5rnm I thickness • 1.0 mm height = 1 5 0 . 0 mm
Strain Gage
l=r°--I Condenser Bank
I
\Coil
Shield
Room
F i g . 1.
An application of the Galerkin method to eq. (9) would be possible if the unknown vector potential is defined element by element in the usual manner as
A(r, z, t) = IN(r, z)] (A(t)) e
(10)
where (A(t)) e is the vector of nodal vector potential at time 't' related to an element named 'e', and [N(r, z)] is the so-called shape function. Here, taking [N(r, z)] as a weighting function, the application of the Galerkin method to eq. (9) produces the next equation,
f__f[N]T{~;(l~r(rA)) r .)1}
~lis equation can be represented simply as follows; [P} {aAlat} + [a] (A} = {J},
(13)
where pp
[e]
= J Jar~N] T [N] dr dz,
[Q]
=ff 1{~ (r[U]T) 1r orfl~.(r[U]) 0 [N]T 3 } + 0z ~ (r/N]) drdz,
(13(a))
(13(b))
+1-~-~z2 82 CrA ) + r/o
Ve
+ro _ a A + L at
X
drdz = 0 , (11)
o
where the matrix [N] T is a transpose of [N] and lie is a domain of the element 'e'. Eq. (11) can be transformed as follows by means of partial integration considering that IN] should be required to vanish on a boundary.
Eq. (13) can be applied to a field of a vacuum where the electrical conductivity 'o' is zero. On the other hand, a finite element equation of motion of a body may be introduced by applying the principle of virtual work together with D'Alembert's principle as follows [7];
tM]~ 3t
{u)+ [K] {u} = (F),
(14)
where
bz j d r d z
+
bA ffo[Nlrr-~-idrdz
[MI = f f [K]
2rrr[Nl r P [NI dr dz,
=ff 27rr[ClT[D] [C] dr dz,
(F} =f f 2nr{[N]T (g X B)} drdz,
(14(a)) (14(b)) (14(c))
and [M] is the mass matrix, [K] is the stiffness ma-
K. Miya et al. I Dynamic behavior of cylinder due to electromagnetforce
404
trix and (F} is the vector of nodal Lorenz force, p is the density of the conductor, [C] is the displacement-strain matrix and [D] is the elasticity matrix. A damping effect is omitted in eq. (14). Eqs. (12) and (14) must be solved simultaneously in order to determine the dynamic electromagnetoelastic field. The vector potential, current density, displacement and nodal force at time t = n • At, corresponding to the nth step, are denoted by {A }n, {J}n, {u }n and {F}n, respectively. A well-known approximate relation between the nth and (n + 1)th vector potentials is:
i{0 /
{A }n +1 = {A )n + ")'
n
+ (1 - 7)
/0:/1 n+l
between displacement, velocity and acceleration:
(o 11
t t 0/.2 Jn + t at2 Jn +IJ
07)
(18)
Substitution of these two equations in eq. (14) gives
([3/] +
aP[K] } [at
jn+l
• At,
= {F}n+ 1 - [K] {u}n + At -~ n (15)
where 7 is a constant and At is the time mesh. With ~' = 0.5, eq. (15) reduces to the usual trapezoidal formula. From eq. (13) we get,
la ul ]
+ kAt 2 t-
7rj. I.
(19)
Test results showed that the term ((au/at) x B)o is very small in comparison with Jo and thus the term was not included in the following calculation.
3. Experimental procedure = [pl-l((J)n+l - [Q](A }n+I).
aA}
~-n+l
3.1. Experimental device
Substitution of these two equations into eq. (15) gives the following equation:
In order to compare the numerical results with experimental ones, the dynamic tangential strain of the cylinder caused by a repulsive force between the cylinder and the solenoid coil inside of it was measured through a strain amplifier made by us for this experiment. The experimental device consists of three parts; an electric power supply, a loading equipment and a measuring system as shown in fig. 1. In fig. 2 are shown the loading equipment (matching impedences and the cylinder) on the front and a condenser bank on the back. The condenser bank (capacitance c = 8 X 5/aF), when charged up to 30 kV, can produce a peak current of 30 kA with a rise time of 20/as. When the gap switch is closed, the current through the solenoid coil produces an impulsed magnetic field, and thus the eddy current in the cylinder which causes the vibration of the cylinder as a result of Lorenz force J X B. The discharged current was measured by a schunt in every experiment. The solenoid coil was made by winding copper wire of 2.5 mm
/1
= Att [P] - 7 [ Q ]
)
{A)n + 7 ( J ) n + ( 1 - ' ) ' ) ( J ) . + l . (16)
The transient electromagnetic field can be obtained by solving eq. (16) successively. The eddy current density induced in the cylinder can be obtained by substituting the instantaneous vector potential in eq. (7), and then the resultant equation in eq. (5). This eddy current density is substituted in eq. (14(c)) to calculated the electromagnetic body force. The Newmark's/~ = ¼ method [8] was applied to successive solutions of eq. (14). It gives the following relations
K. Miya et al. / Dynamic behavior of cylinder due to electromagnet force
405
Fig. 2.
radius with a vinyl insulation around a supporter of bakelite. Its height and diameter are 0.15 m and 0.043 m respectively, and its inductance is 6.6/~H. The cylinder is made of brass (70% Cu, 30% Zn). It has 0.049 m inner diameter, 0.15 m height and 0.001 m thickness. The elastic modulus and Poisson's ratio are 1.37 X 10 ~° Pa and 0.374 respectively, and the electrical conductivity 28%.
enclosed by a filament of it is minimized. A dummy gage was put very near the active gage to cancel induced currents in the two gages in the differential strain type amplifier. Furthermore, the strain gage is shielded by copper foil in anticipation of the effect of the magnetic shielding. The lead wires of the gages were twisted to make the receiver loop as small as possible.
3.2. Strain measurement
4. Numerical analysis and experimental results In measuring the strain in a high, transient electromagnetic field, special care was taken to proper!y ground and shield the equipment. The grounding had to be done so as to reduce the noise made by the high electric potential near the measuring point. This was accomplished by containing the measuring system, which has an isolated power source, in a shielded room and groundin~ it near the measuring point so as. to minimize the noise voltage from ground currents flowing through a common impedence. As to the shielding, it is needed to reduce the induced current in the strain gage and lead wires due to the changing magnetic field. The best way to do this is by decreasing the area of the receiving loop. To this end, a non-inductive strain gage developed by Kyowa Dengyo Co. Ltd. was used. A schematic figure of the gage is shown in fig. 1. It can be seen that a loop
4.1. Numerical results
The cylindrical specimen with dimensions is shown in fig. 1. As shown there the 24 turn solenoid coil was around a grooved bakelite bar and set inside the cylinder the eddy current flowing in the tangential direction only and the transient magnetic field must be analysed with a computer code based on eqs. (16) to (19). A triangular mesh division of a half area of the cylinder, coil, and air is shown in fig. 3. Only a half area is required because the eddy current and the magnetic field are symmetric with the lateral center line (r-axis). This means that the boundary condition 8A/Sz = 0 should be applied on the lateral center line. This boundary condition corresponds to an adiabatic one. Vector potentials on the z-axis are always zero
406
K. Miya et al. / Dynamic behavior o f cylinder due to electromagnet force (mm) 4O0.0
¢.(z)
150.0
tO0.C 84.(
c)
75 .( in
66.C 55.C 44J
::.'.'S&~ 22 .C t f.O ~(r)
O.C J
O. 0
20.0 24.5 25.5
50.0
1500
I
4000 (ram)
Fig. 3.
for an axi-symmetric distribution. The two remaining boundaries of the domain in fig. 3 require the vector potential on them be zero as they are located far from the external source (the solenoid). The shape of the current supplied to the solenoid from the condenser bank is also shown in the figure, and was used as the input for the numerical calculation of the eddy current. In this case the permeabilities of the cylinder, the coil and the air are all equal, but the electrical conductivities are zero except for the cylinder. The conductivity of the coil should be zero since the external current in the coil includes the effect of induction. The model of the coil in fig. 3 does not correspond to the actual one, but a magnitude of the external current was determined in such a way that the total ampereturns are equal for both cases. In fig. 4 is shown an example of the numerical results obtained for the eddy current at the center of the cylinder. The transient external current is shown upper most. The numerical result obtained from an another computer code based on an inductance method using elements divided in the axial direction only is also shown there. A complete agreement between
both results from the finite element method and the inductance method (which is not stated here) cannot be expected. This is because the thickness of the cylinder was divided into four layers, and the solenoid was approximated with a shell for the case of the finite element method; while the thickness was not divided and the solenoid was as it was in the case of the inductance method. Relative errors between both methods are within 10%; showing a validity of the code based on the finite element method. It goes without saying that a time when the induced current density shows a maximum corresponds exactly to one at which the external current is at a maximum. The distribution of the induced current in the axial direction is shown in fig. 6 at the time when it reaches the maximum. Numerical results by the inductance method are also shown in the figure with a parameter for the number of mesh divisions in the axial direction. The abscissa runs from the end of the cyliner to the center, corresponding to a half length of the cylinder. The induced current density obtained from the inductance method converges to an accurate value as the number of mesh divisions increase. Upper points cor-
407
K. Miya et al. /Dynamic behavior o f cylinder due to electromagnet force
-
Imax :14300A ;~ To
r,,. ,..,. o
_1
5
:
q
I0
20 30 40 TIME (/~sec)
50
60
0.5
4J
,¢
%
o
>t~-O.5 la.i
O I-
1o"
zo
k
30
I
I
I
40 I
TIME (/z see)
50..'/~tz~/'-
,~
~\ by inductancemethod / 7 ",~ "-..,I '% sr S
-
-1.0 r~ t~ D O
a -1.5
t/J
7 ~ --2.0
~%N
ssSSIS
~ ~ e m e n t
method
~g. 4.
. abcissa: 1 ms/div. ordinate: 0.5 V/div.
respond to locations nearer to the coils while lower points to locations farther from the coils. The ratio of the mean amplitude of the oscillation to the mean in the range near the center is about 1.5%. The induced current density near the end of the cylinder is a little lower than that near the center. This is caused by the leakage of the magnetic flux from the end. In fig. 5 are given two examples of the experimental results with the same external currents (14 300 A)but different time scales. 4.2. Comparison o f numerical and experimental results
In fig. 7 are shown numerical and experimental results for comparison. Dynamic strains from analysis and experiment do not agree completely. However, periods from analysis and experiment are both nearly equal to 40/as. This period can be calculated approxi-
Fig. 5.
mately for a thin cylinder following eq. [10], rp = 27rR(p/E) '~ ,
(20)
giving Tp = 40.6/as with R = 0.025 m, p = 8560 kg/m 3 and E = 1.28 X 1011 Pa. The period obtained from the experiment is 40/as and that from the numerical analysis is 39/as. Both of these periods agree quite well with each other. On the other hand the distribution of the induced current in the cylinder wall must usually be confined to a radial depth less than the wail thickness. The extent of penetration of the magnetic flux, and so the current, should be provided by so-called skin
408
K. Miya et al. /Dynamic behavior o f cylinder due to electromagnet force
depth 6 given by, 6 = (2t~//.tco)1/2 ,
(21)
where K is the resistivity and co the angular frequency of the first cycle of the transient external current. Introducing the physical constants of the material used here in eq. (21), namely K = 6.14 X 10 -a,/~ = 4rr X 10 -7 and co ~-8.3 X 103, a result of 6 =3.4 X 10 -3 (m) is obtained. Thus the skin depth given by eq. (21) is larger than the wall thickness of 1 × 10 -a (m). Therefore, the cylinder does not tend to act as a magnetic shield confining the magnetic flux to the area of the gap between the coils and the cylinder. However, when the skin depth is smaller than the thickness, the current 11 in the single-turn tube is calculated as follows [6],
t~ <
% x
g D
E
R2 I , = I2N2 ( R - h/2) 2 '
20
40
Distonce,
60
where I2, R2, N2 are the current flowing in the solenoid coil, its radius and the number of turns on the coil, respectively; and h is a thickness of the cylinder. Baines et al. [11 ] calculated the uniform radical pressure at a sufficient distance from both ends of the tube as follows considering the tube and the coil as infinite parallel plane conducting sheets under the condition of no-leakage of the magnetic flux through the tube wall (6 is small compared to its thickness),
75
~ (x103m)
Fig. 6.
1000
experimental
result
/ N ~
P = laoN2Il12/2L z , C
l
i
l
I
I
(23)
where L is tile length of the tube. Dynamic displacement due to this magnetic pressure can be approximately given by a solution of the following equation of motion of the thin cylinder [10],
i 0
'~O - 5 0 C
.E
(22)
a2u/at 2 + co2u
40C
2
=f(t),
(24)
where 20C
/&
g
-40(
•
i
nume r e s u l t (by F . E . M . ) ~ . /
-60C
Middle p o r t ~'~
Edge port
Time
( bL s e c )
Fig. 7.
]
J
1 co = - ~ ( E / p ) l a ,
(25)
and f(t)is given by eq. (23) taking into account the shapes of the time varying currents I i ( t ) and I 2 ( t ) The maximum strain predicted from eq. (24) is about 710 #. This gives an over estimate since no-leakage of the magnetic flux is assumed. Peak strains in the oscillation of the tangential strain shown in fig. 7 are not constant because the axial distribution of the electromagnetic force is not uniform, as suggested by the
409
K. Miya et al. /Dynamic behavior o f cylinder due to electromagnet force
results shown in fig. 6. The wave shown in fig. 7 involves the effect of a wave propagating in the axial direction, thus resulting in the slightly complicated wave shape shown in fig. 7. Peak strains given by the finite element analysis are between 300/a and 650 tx as shown in fig. 7. The experimental results show that those strains are between 600/a and 1000/x where the peak of the current shown there is 14 300 A. As stated before, agreement as to the period is quite good between the finite element analysis and the experiment, but as to the strain it is not so good. As technical problems such as the noise reduction and rate dependence of material constants are not completely resolved, a further study will be needed. In fig. 8 is shown the relation between the peak strain and the external current of the solenoid coil. As the magnetic pressure is proportional to the product of the currents 11 and 12, with I1 to 12 given as shown in eqs. (23) and (22) respectively, the induced strain should be proportional to a square of the external current 12. On the abscissa is plotted various external currents obtained from the relation between the charged voltage and the current. Experimental data are closely approximated by parabollic curves. Three parabolas correspond to the peak occurred at different times as shown in fig. 8.
5. Concluding remark The finite element method rhay be one of the most powerful methods to analyse the rather complicated static and dynamic deformation of metallic components due to an electromagnetic force. An example of the finite element scheme on the axi-symmetric problem, which couples the electromagnetic field with the equation of motion of the axi-symmetric body, is shown in this paper. An application of the scheme to plastic work of metal components has not yet be done up to now, but seems very important, and possible, from the viewpoint of metal forming. A measurement system of the dynamic strain particularlyunder a transient strong magnetic field, is also very important in the fields of metal forming and structural design of a fusion power plant. The experiment reported here is considered useful for the establishment of such a measurement system.
Acknowledgement This work has been done with a financial support from the Controlled Thermonuclear Reactor Committee (Chairman, Prof. Akira Sekiguchi), Faculty of
200C
:=L
C
iL ,llVllVVVi
150C
--
IOO3
C
500
-.-
.//./
•
, , =
I
I
tO Coil
#5
current Fig. 8.
,
Imox
20 (kA)
25
i
410
K. Miya et al. / Dynamic behavior of cylinder due to electromagnet force
Engineering, University of Tokyo. The authors would like to thank Prof. Y. Ando, University of Tokyo for his encouragement to them. The authors also would appreciate the discussion with Dr. F.C. Moon, Cornell University and Dr. J.E. Akin, University of Tennessee on this work, and thank Mr. K. Someya for his sincere help with this experiment.
References [1] For example, Proc. 7th Symp. Engrg. Problem of Fusion Research, Knoxville, TN, USA, Oct. 1977, pp. 13331397. [2] B. Misra and V.A. Maroni, Thermal hydraulic analysis of two fusion reactor first wall/blanket concepts, ibid, pp. 1459-1464.
[3] F.C. Moon, Mechanics Today (Pergamon Press, Oxford, 1978) V01. 4, pp. 307-390. [4] M.J. Forrestal and Overmier, AIAA Journal 12 (5) (May 1974) 722-724. [5] H.C. Wallingand M.J. Forrestal, AIAA Journal, 11 (8) (August 1973) 1196-1197. [6] G.K. Lal and M.J. Hillier, Internat. J. Mech. Sci. 10 (1968) 491-500. [7] W.K.H. Panovsky, Classical Electricity and Magnetism (Addison-Wesley,Reading, 1961). [.8] O.C. Zienkiewicz, Finite Element Method in Engineering Science (McGraw-Hill,New York, 1971). [9] N.M. Newmark, A Method of Computation for Structural Dynamics (ASCE, 1959). [10] K. Miya and J. Silverman, Thermomechanical dynamic behavior of the solid wall of a pellet fusion reactor, to appear in J. Nucl. Technol. (1980). [11] K. Baines, J.L. Duncan and W. Johnson, Proc. Instn. Mech. Engrg. 190 (1966) 348-362.