**2. Numerical Formulation**

The 3D numerical analysis of multi degrees of freedom electromechanical devices represents a challenging, still open problem, which poses several critical issues. Several commercial codes, typically based on the Finite Elements Method (FEM), provide packages for coupled electro-mechanical analysis, but they seem to be not very effective with multiple degrees of freedom and with sliding contacts. Moving conductors and sliding contacts usually require remeshing of the domains with consequent increase of computational times and potential numerical instabilities.

Integral Formulations (e.g., the Method of Moments) seem to work quite well with moving domains since they require the discretization of the active regions only (namely conductors and ferromagnetic materials), so avoiding the problem of coupling meshes with different speeds. Integral formulations implicitly enforce the far field boundary conditions, and are able to produce accurate results by using coarse discretization (when compared with those required by FEM).

The numerical investigation of the complete launch system has been performed by the research code EN4EM. It is based on an integral formulation, and it is under continuous development by the authors for investigating electromechanical systems [27–34].

EN4EM applied to the launch package is able to simulate the whole system considering the characteristics of the two devices and taking into account their "strong interaction". Figure 1 shows the steps of a conceptual flow chart of the numerical formulation: discretization in elementary volumes, writing and integration of Ohm's law, arrangement in an electric network and writing of the governing electro-mechanical equations. The usual notations for the electromagnetic quantities in Figure 1 are adopted according to the Table 1.

With reference to the inset in the bottom right of Figure 1, Ohm's law is written in the conductive elementary volumes where a uniform current distribution is assumed (row #1 of the inset). The additive properties of the integrals with respect to the integration domain allows expressing the fields and potentials in the *k*-th conductive elementary volume as a summation of contribution due to the currents in the other volumes (row #2 of the inset). Finally integration on the *k*-th elementary volume leads to an equation that can be seen as the voltage-current relationship of a branch that is a series connection of a resistor, an inductor coupled with other inductors, and a voltage generator controlled by the currents in other elementary volumes (row #3). The equivalent network of the device is built by connecting the terminals of the branches so obtained. In case of a multi-device system, the code is able to model separately the different devices by their equivalent networks. Then, all these networks are connected together according to the relative positions of the corresponding elementary volumes and to the presence of electrical contacts between them. The described procedure has been applied to model the launcher and the compulsator, obtaining a whole model composed of thousands of branches.

**Figure 1.** Conceptual flow chart of the numerical formulation. (Reproduced with permission from [30], IEEE, 2017).

**Table 1.** Notations.


Mesh analysis yields to the governing equations written in matrix form:

$$\mathbf{L}\left(\mathbf{C}(t)\right)\frac{d\dot{\mathbf{i}}}{dt} + \left(\mathbf{R} + \mathbf{K}\{\mathbf{C}(t), \dot{\mathbf{C}}(t)\}\right)\dot{\mathbf{i}} = \mathbf{e}(t) \tag{1}$$

The values of the elements of the matrices in (1) are function of the system configuration *C*(*t*) and its derivative . *C*(*t*) (i.e., the relative positions and velocities of the elementary volumes used to discretize the devices respectively). Coupling between electrical and mechanical equations is achieved by the terms *vk*(*t*) × *Bk*(*t*) that, once integrated on the elementary volumes, are assembled to form the matrix of the motional terms *K*, and by the terms *jk*(*t*) × *Bk*(*t*) which are integrated to provide the forces and the torques on the moving parts of the device. The mechanical equation for the armature of the launcher is:

$$m\ddot{\mathbf{x}}\_{\rm G} = F(\mathbf{i}, \mathcal{C}(\mathbf{t})) \tag{2}$$

while the equation for the compulsator are:

$$\mathbf{I}\_{\mathbf{G}}\dot{\boldsymbol{\omega}} + \hat{\boldsymbol{\omega}}\mathbf{I}\_{\mathbf{G}}\boldsymbol{\omega} = \mathbf{M}\_{\mathbf{G}}(\mathbf{i}, \mathbf{C}(t)) \tag{3}$$

In the above equations *F* represents the resultant force on the launcher armature and *m* is its mass, *M<sup>G</sup>* is the resultant torque on the rotor of the compulsator and *I<sup>G</sup>* is its inertia tensor. The coupled

differential equations for electrical and mechanical equilibrium are time varying, and the resulting system is nonlinear; integration is carried out as described in [27,28].

It is worth to mention that the currents in the branches of the equivalent network are not fictitious currents. The current in a branch of the equivalent network is in correspondence with one component of the current density of an elementary volumes as shown by the insets in the upper right and in the lower left corners in Figure 1. By the above described equivalent network the distribution of all the relevant electromechanical quantities can be evaluated everywhere in the whole system. In particular, the expression of the force by the term *jk*(*t*) × *Bk*(*t*) allows to evaluate its distribution in each elementary volumes and to determine the contribution to the total dynamical action in materials that are usually heavily stressed when used in EML technology.

Simultaneous analysis of interacting devices simply requires building a bigger equivalent network constituted by the properly connected equivalent networks of the components. This implements a "strong interaction" between devices for which a "strong coupled" analysis is performed.

If the devices can be assumed magnetically uncoupled (i.e., leakage magnetic fluxes from a device that links the others are negligible) independent networks are built, linked together at the common terminals and simultaneously solved. It is worth to remark that the topology of the resultant network is the same than that of the more general case of magnetically coupled devices. The only difference lies in the filling of the inductance matrix *L* and of the motional terms matrix *K*, which are more densely populated in the latter case.

Coupling of the equivalent network so far described with external circuits is straightforward and this will be performed in the next section to investigate the "weak interaction" between the rail launcher and the compulsator.

The presence of sliding contacts has been taken into account by introducing an auxiliary network [20]. Considering that the motion in the rail launcher is characterized by a straight trajectory, all the possible contacts between elementary volumes on the rails and on the armature are known a priori. A new branch is set for each couple of volumes (one on the inner parts of the rails and the other on the faced outer surfaces of the armature) that can have a contact during the motion.

Figure 2 shows an example of a set of auxiliary branches that take into account the sliding contacts. The circuit elements in the auxiliary branches are function of the shared portion of the faces (if any) between the two volumes. When the shared portion is zero, the auxiliary branch opens. The proposed model of the sliding contacts is coherent with the adopted formulation based on an equivalent network of the entire system.

**Figure 2.** Example of the auxiliary network used for taking into account the presence of sliding contacts. (**a**) All the branches are shown. (**b**) Only the branches that are not open circuits are evidenced.

Finally the introduction of the equivalent network allows the use of advanced analysis techniques for the evaluation of the sensitivity of the response with respect to parameter variation [35]. This represents an important tool in gradient based optimization processes.

EN4EM was validated by comparison with results produced by other numerical codes and with experimental data [5,6,27,28,31,32].

#### **3. Weak Interaction Analysis**

In this section the analysis based on the "weak interaction" between rail launcher and compulsator will be critically reviewed and the difficulties of the equivalent circuit extraction will be discussed.

Figure 3 shows the active parts of a single-phase, two-pole compulsator with selective passive compensation provided by a discontinuous conductive shield. The inner part of the device consists of two stationary field coils which produce a magnetic flux density distribution whose axis is in the vertical direction (Figure 3b). The two series connected armature windings are located on the rotor which is the outer part of the machine. Figure 3c shows the stationary compensating shields that are in the central part of the machine. Figure 3a shows a 3D view of the device. A more detailed description is reported in [29].

**Figure 3.** An example of the modeling of a compulsator by EN4EM. Snapshot of (**a**) the 3D view; (**b**) cross section of the machine; (**c**) discretization of the compensating shield. (Reproduced with permission from [29], IEEE, 2017).

When using a compulsator to feed a rail launcher, its rotor has to be preliminarily driven to a proper angular speed at no load conditions; part of the stored kinetic energy of the rotor will be delivered to the railgun armature. At the firing instant, the armature windings of the compulsator are connected to the rails of the launcher. During the launch, the rail armature accelerates and at the same time, the rotating part of the compulsator decreases its angular speed. The rates of change of the speeds of both the moving parts (and consequently their positions) are not known a-priori and substantially depend on the delivered current.

Since the simultaneous 3D analysis of the two electromechanical devices can be very time consuming, a common practice is to substitute one or both the device with equivalent circuits. This procedure may lead to significant error in the current flowing in the devices because of the difficulties in the determination of the topology and of the parameter extraction of the equivalent circuit.

Referring to the compulsator, when looking for a lumped equivalent circuit, we have to distinguish between the two compensation techniques usually adopted. When compensation is achieved by the use of shorted discrete coils, we can consider the self and mutual inductances of all windings (field, armature and compensating ones). Considering that the values of the mutual inductances depend on relative angular positions only, the lumped equivalent circuit of the compulsator is able to give accurate results; its building is a matter of evaluation (experimental or by computations) of self and mutual coefficients and EN4EM can be used as an extraction tool.

Things go in a different way when non-uniform (discontinuous) compensating shields are used. In these devices, the extraction of the values of the equivalent internal inductance *L* is not as straightforward as with shorted discrete coils, since it is function of the position of the rotor, as well as of the speed of the rotor, which is not known and varies during the system operation.

In fact, the speed of the rotor imposes the frequency of the electrical quantities in the shield, which in turn determines the path and the amplitude of the induced currents and their shielding effects which contribute to the internal equivalent inductance. Figure 4 shows an example of eddy currents distribution on one of the discontinuous compensating shields of the compulsator above described whose operating conditions are reported in [29]. In particular the compulsator was loaded with a simple lumped R-L equivalent circuit representing the rail launcher. The rectangle in Figure 4 represents the flattened cylindrical surface shown Figure 3c.

**Figure 4.** Eddy current distribution on the conductive compensating shield. (Reproduced with permission from [29], IEEE, 2017.)

The equivalent internal inductance of the compulsator, as well as the electromotive force at the terminals of the machine, are functions of the distribution of the currents on the shields. Expressing these functional dependencies represents a challenging problem.

Similar difficulties arise when looking for an equivalent circuit of the rail launcher. Figure 5 shows a rail launcher fed by an ideal voltage source whose waveform is given by *e*(*t*) = 50 1 − *e* −*t* τ *V*. Details about the geometry are in Table 2. In this figure, 20 elementary volumes located near the launcher breech are shown. The device was analyzed by EN4EM and the current density waveforms in the evidenced volumes are reported in Figure 6. The solid curves are associated with the labelled sections in reported in the inset of Figure 5; the greatest current density occurs in the section #5 in the inner part of the rail. The peak values became smaller as the outer boundary of the rail is approached. The reported waveforms put into evidence the uneven current density distribution in the rail section for most of the launch time, the ratio between the maximum and the minimum value is greater than two. This behavior is a consequence of the skin and proximity effects; also, the velocity skin effect has an heavy impact on the current distribution, which is a function of the rate of change of the electrical quantities imposed by the generator and the rate of change of the mechanical quantities i.e., the speed of the armature [36,37]. Both these quantities change during the launch.

**Figure 5.** A sketch a rail launcher. The inset put into evidence the cross section of half of the upper rail. The labelled rectangles represent the cross section of the inner layer of the elementary volumes, i.e., those involved in the sliding contacts with the armature.


**Table 2.** Description of the devices.

The basic equivalent circuit of the railgun, here shown in Figure 7, is usually composed of three elements. A resistor and an inductor, both varying with the distance travelled by the armature, and another resistor, that takes into account the motional induced electromagnetic force, is related to the inductance gradient and its resistance depends on the speed of the armature.

In the light of the above discussion, these circuital components should be function of the frequency which influences the effective current density distributions in both the rails and the armature, as well as of the time, because of the speed and of the distance travelled by the armature. Some lumped equivalent circuits of the rail launchers adopt simplified expressions of the form: *R*(*z*, *f*) = *Ra* + *R*<sup>0</sup> + *R z* and *L*(*z*, *f*) = *L*<sup>0</sup> + *L z*, where *R* = *dR dz* is the rail resistance gradient, *<sup>L</sup>* <sup>=</sup> *dL dz* is the rail inductance gradient, *Ra* is the armature resistance and *R*<sup>0</sup> and *L*<sup>0</sup> are the resistance and the inductance due to the connection wires [1]. All these expressions discard the dependence of these parameters on the frequency and are

not able to take into account the velocity skin effect. On the other hand, the formal definition of the parameters as shown in Figure 7 and their estimate pose complex problems.

**Figure 6.** Time waveforms of the *z*-axis component of current density in the rail at the feeding end. Uneven distribution occurs for about 80% of the launch time.

**Figure 7.** Basic equivalent lumped parameter circuit of a rail launcher. The presence of f (frequency) among the independent variables means that the circuit parameters are function of the rate of change of the electromagnetic phenomena.
