Advanced Model of Electromagnetic Launcher

An advanced 2D model of electromagnetic launcher is presented respecting the influence of eddy currents induced in the accelerated ferromagnetic body. The time evolution of electromagnetic field in the system, corresponding forces acting on the projectile and time evolutions of its velocity and current in the field circuit are solved numerically using own application Agros2d. The results are then processed and evaluated in Wolfram Mathematica. The methodology is illustrated with an example whose results are discussed.


Introduction
Devices working on the principle of electromagnetic launching have been known more than 150 years. We can mention, for example, electromagnetic actuators of various constructions, parts of switching gears, electromagnetic guns, and many others.
But despite their relatively simple construction, their design and modeling is not an easy business because of the presence of nonlinear structural elements; some of them, moreover, move. Even when the physics of the system is known, the problem has not been solved completely so far because of strongly nonlinear interactions of involved magnetic field and two circuits (electric circuit and mechanical circuit).
For example, the first attempts to design and build a magnetic coil gun were only based on more or less sophisticated experiments, which is the way able to provide a functioning device, but mostly with poor op-eration parameters. This part of a rather old history is well summarized and documented in [1].
A more extensive development in modeling such devices began in eighties of the last century, which was possible due to achievements in hardware and software capabilities. But still, there appeared serious problems with the numerical solution due to strong nonlinearities involved in the task. Therefore, the first models were solved in the weakly coupled formulations, i.e., the particular equations describing the system were calculated separately. This approach could be considered only orientative and the results often exhibited unacceptable errors compared with experiments.
More reliable methodologies and also the first attempts to optimize selected structural parts started appearing in nineties of the last century and at the beginning of this century. Worth mentioning are several works solving the problem in the quasi-coupled formulation, which means that the mathematical model was solved iteratively as a whole, but specific nonlinearities were only respected in selected time steps (in the remaining time steps the system was considered linear). Information about these more advanced approaches can be found in [2], [3], [4], [5], [6], [7], [8] and others.
The presented paper represents an attempt to solve the problem respecting all nonlinearities in the system together with the influence of eddy currents due to both transformation and motion.

Formulation of Technical Problem
Consider the arrangement in Fig. 1. The system consists of a field coil carrying current of an appropriate time evolution, a ceramic or plastic leading barrel and a ferromagnetic body (projectile) that is to be accelerated. Geometry and physical properties of the above parts are supposed to be known, as well as the parameters of the feeding electric circuit.
The time variable magnetic field generated by the field coil produces magnetic force acting on the projectile that starts moving in the direction indicated by the arrow. The field coil must be switched off at the moment when the projectile is approximately in the middle of the coil. Then it continues moving just by its own inertia.
The aim of this work is to model the movement of the projectile (particularly its acceleration and velocity) with the help of computational software Agros2D [9], [10] and Mathematica 7 [11].

Mathematical Model
The basic mathematical model of electromagnetic launcher consists of a nonlinear partial differential equation describing the magnetic field in the system (taking into account the eddy currents due to the velocity of motion) and two highly nonlinear ordinary differential equations for the time evolution of the field current and velocity of motion of the projectile. The corresponding equations are presented and in the following text.
The electric circuit (see Fig. 2) obeys the equation with the initial conditions Here, R denotes the overall resistance of the circuit, L is the inductance of the field coil that is a nonlinear function of the position s of the projectile, its instantaneous velocity v, and instantaneous value i of the field current, C stands for the capacitance of the battery, and U C0 represents its initial voltage. Finally, L 0 = L ( s start , v = 0, i = 0) is the initial inductance of the system, s start denoting the starting position of the projectile.
The time evolutions of velocity and trajectory of the projectile are governed by the equations where m is the mass of the projectile, v is its velocity in, F em is the electromagnetic force acting on it and F dr denotes the sum of the drag forces (that are given by friction in the leading barrel and aerodynamic force). The quantity F em , analogously to the inductance L, is also a nonlinear function of the instantaneous position s of the projectile, its velocity v and instantaneous value i of the field current. Both the above quantities, F em and L have to be calculated from the current distribution of magnetic field in the system.
The distribution of magnetic field may be described using several formulations based on various potentials. In our case, the magnetic vector potential A is used as the principal quantity. The full equation for this quantity may be written in the form [12] curl where µ denotes the magnetic permeability, γ is the electrical conductivity, and J ext is the field current density in the field coil. The boundary condition along a sufficiently distant boundary is of the Dirichlet type.
The term −γ ∂ A ∂t − v × curl A represents the total eddy current densities induced in the system. The first term −γ ∂ A ∂t is eddy current densities due to transformation, the second term γ v × curl A is eddy current densities due to the movement. Their effects act against the magnetic field and lead to a slight deceleration of the projectile.
The magnetic force F em acting on the projectile is determined using the magnetic Maxwell stress tensor T em . The vector of the force is in our case given by the integral [13] where the integration is performed over the whole boundary of the projectile. The magnetic Maxwell stress tensor itself is given by the relation where H and B are the vectors of magnetic field strength and magnetic flux density, respectively ( B = curl A, B = µ H), I stands for the unit matrix and symbol ⊗ represents the dyadic product.
The inductance L occurring in Eq. (1) can be calculated from the expression where W m is the energy of magnetic field of the system that may be calculated using the formula where V is the volume of the definition area.
As the arrangement may practically be considered axi-symmetric, both quantities A and J ext have only one nonzero component in the tangential direction, while the trajectory s, velocity v and both forces F em and F dr have only one nonzero component in the xdirection (see Fig. 1).

Numerical Solution
The task was solved numerically by means of our own application Agros2D and commercial software Mathematica 7. The code Agros2D (that is completely free) is based on a fully adaptive higher-order finite element method that is intended for numerical solution of 2D nonlinear and nonstationary multiphysics problems described by a set of partial differential equations. It is characterized by a number of quite unique features such as finite elements up to the 10th order, efficient multimesh technology, hanging nodes of any level, combination of various elements including curvilinear elements (for approximation of curvilinear boundaries and interfaces) etc. The evaluation of the results and some auxiliary computations were carried out in Mathematica 7 using a lot of own procedures and scripts. The whole algorithm consists of the following steps.

Inductances and Magnetic Forces
Computation of the dependencies of the inductance L and x-component of the force F em,x on the field current i and position x for particular velocities v. In this way we will obtain two sets of nomograms (for particular velocity) from which we can easily find the values of both above quantities for values i, x and v x (and, where necessary, also their derivatives). The final values of L and F em,x are for any current, position of the projectile and its velocity calculated using linear interpolation in the above nomograms.

Circuit Equations
Equation (1) and Eq. (2) are solved simultaneously but successively, in time steps.
• Solution of Eq. (1). First it is necessary to modify the term d(Li)/dt.
For further processing, we determine the second derivative of term Li with respect to time: After another differentiating Eq. (1) with respect to time t and substituting for d 2 (Li)/dt 2 we obtain Introducing the differences instead of the derivatives and substitution g = di/dt leads to a system of difference equations in the form: and with the initial conditions Here, a denotes the acceleration, k is the number of step and ∆t k stands for the k-th time step.
• Solution to Eq. (2). No modification of this equation is necessary.
• The values L, dL dx and F em,x can easily be determined from the above nomograms using an appropriate interpolation or extrapolation technique. On the other hand, the nomograms must be appropriately dense.

Numerical Problems Connected with Solution
The principal complication in the numerical solution is computation of the nomograms for higher values of the field current. High field currents bring about either partial or complete oversaturation of the ferromagnetic projectile. But this oversaturation may be evaluated only approximately, because the magnetization curve of material in the relevant domain can only be estimated. This leads to non-estimable errors of results, and, moreover, the related iterative processes are often accompanied by various undesirable phenomena such as oscillations, and generally, a very poor convergence rate, which means a long time of computations.

Illustrative Example
The mathematical model was tested on a portable three-stage electromagnetic accelerator depicted in Fig. 3 that was completely developed and built by the second author. Both the model and measurements were realized just on the first stage of the device.
The principal aim of the task was to compare the measured and calculated velocity and trajectory of the projectile.

Input Data
The principal dimensions of the coil of the first stage, leading barrel and projectile are shown in Fig. 4. Their values are: D 1 = 36 mm, d 1 = 8 mm, L 1 = 50 mm, D 2 = 6.75 mm, L 2 = 52 mm and initial position of the end of the projectile x start = −38 mm.
The voltage U C0 of the capacitor battery is 350 V. The parameters of the electric circuit were measured for frequency 156 Hz, because the shape of the corresponding current corresponds to the shape of real current in the series RLC electric circuit. For the above frequency, the capacitance C = 7.11 mF, the total resistance of the circuit R = 0.145 Ω and inductance of the coil L = 0.220 mH. The coil contains 203 turns wound in 7 layers, each having 29 turns. The conductor is of circular cross section and its diameter is 1.6 mm. The leading barrel is made of brass with a longitudinal gap, so that it has no effect on electromagnetic quantities.
The projectile is made of material Vacoflux 48 (produced by German company Vacuumschmelze GmbH & Co). After manufacturing, the surface of the projectile was hardened, but the conditions of hardening are not known. That is why we used the curve µ r = µ r (B) of its primary magnetization that is depicted in Fig. 5. Electric conductivity of this material is approximately 2.5 MS·m −1 and its mass density is 8120 kg·m −3 .

Results
For the accuracy of results the decisive role is played by the correctness of the time evolution of current i in the field coil. This current was measured and also modelled using Eq. (1). The results are presented in Fig. 6.
The blue line (measured) shows the measured current that was switched off at the moment when the projectile reached the middle of the coil. The black line (calculated 1) shows the evolution for circuit parameters given in subsection 5.1. (here, however, the current was not switched off). And just for comparison, the red line (calculated 2) shows another calculated evolution for circuit parameters determined by their measurement at the industrial frequency, i.e., not quite correctly. Obviously, the accordance of the blue and black lines is extremely good.
Important are the nomograms of the inductance L, its derivative with respect to coordinate x, and mag-  In fact, the figure contains two nomograms. The transparent one shows the case when the relative permeability of the projectile corresponds to Fig. 5, but ends at the value µ r = 4.6 for B = 1.5 T (for lower values of µ r there appeared severe problems with the convergence). The real nomogram, however, is the lower one. For its construction, a special restricting function had to be proposed, being able to simulate even lower values of the relative permeability. In case of higher values of the field current i, the projectile becomes strongly oversaturated and its relative permeability drops to 1.
The analogous nomograms for the derivative of inductance L and magnetic force F em,x are depicted in Fig. 8 and Fig. 9.
The nomogram containing the derivative dL/dx is not as smooth as it should be. The visible oscillations are caused by numerical errors due to problems with the convergence of the iterative processes in the domain of a high saturation of the projectile. A better level of smoothness could be achieved by using interpo-    lation polynomials of higher orders, but at a cost of a substantially longer time of computation. On the other hand, many tests proved that these oscillations do not affect the resultant evolution of the field current too significantly.
The time evolution of motion of the projectile was calculated using Eq. (2), where the drag force is represented by the aerodynamic resistance. Its value is given by the formula where ρ is the density of air, S is the area of the cross section of the projectile and C is a constant that depends on the shape of the projectile. In our case   Fig. 11 shows the time evolution of its velocity. The time step ∆t = 10 −6 s. The red lines are the dependences without respecting eddy currents in the projectile, the blue ones represent them. The differences are, however, very  small. The eddy currents decelerate the velocity of the motion, but only very slightly.
The average muzzle velocity measured on the device was 28.9 m·s −1 . The difference from the calculated values is about 1 m·s −1 (over 3 %), which can be considered quite excellent.

Conclusion
The paper is aimed at numerical modeling of operation characteristics and parameters of an electromagnetic accelerator.
The proposed mathematical model consisting of one partial and two ordinary strongly nonlinear differential equations was solved numerically in the quasi-coupled formulation, with pre-calculated nomograms of the inductance of the system, its derivative and magnetic force acting on the projectile. The computations had to be carried out extremely carefully because of poor convergence due to oversaturation of projectile at high values of the field current.
The model proved to provide realistic results whose agreement with the measured data is excellent. Next work will be focused on modeling all three stages of the electromagnetic gun.