A Comparison of Velocity Skin Effect Modeling With 2-D Transient and 3-D Quasi-Transient Finite Element Methods

The analysis of the velocity skin effect (VSE) in electromagnetic launchers (EMLs) requires a 3-D transient finite element method, unlike magnetic skin and proximity effects. However, VSE is dominant at high speeds, and this creates convergence problems when moving or deformed mesh physics is used in a transient FEM in the 3-D analysis. Commercial finite element software cannot solve the electromagnetic aspects of such a high-speed application with a transient solver in 3-D. Although 2-D approximations can be used, such an approximation overestimates VSE resistance due to geometry simplifications. In this study, we proposed a novel quasi-transient 3-D FEM model where the air-armature region’s conductivity is varied to emulate the high-speed motion of the armature. Results showed that the 2-D approximation overestimates the VSE resistance by almost 40%. The proposed VSE model has been included in the EML model, and simulation results are compared for experimental results with different EMLs, EMFY-1, and EMFY-2 and showed good agreement.


I. INTRODUCTION
E LECTROMAGNETIC launchers (EMLs) utilize the Lorentz force to accelerate projectiles up to a few thousand m/s. Although the skin and proximity effects are well-known phenomena in electromechanical systems, the velocity skin effect (VSE) is an EML-specific phenomenon that occurs in the rails at the vicinity of the armature, as stated in [1]- [3]. To observe the VSE on the current distribution, transient finite element method (FEM) simulations are required to calculate the magnetic field distribution at the armature and rail interface, which is quite complicated to compute [4]. Moreover, due to large armature velocity, transient 3-D FEM has a convergence problem due to interpolation errors between stationary and moving nodes. Thus, special FEM codes are required since commercial software is only used for static mesh elements in simulations due to their formulations [4]. Manuscript  EMAP3D code uses the Lagrangian description and rail reverse motion techniques [5], which is also called the upwinding method to simulate an electrical sliding contact problem under transient conditions. EFFE software [6], [7] also has the same technique: rail moves backward rather than the armature movement. Quasi-static, time-harmonic FEM are also used to model velocity induced currents [8], [9]. Transient FEM simulations that have an armature movement have two computational domains: stationary and moving. Discontinuity in the field distributions at the intersection of the domains is handled using the continuity conditions. However, this is not the only strategy. Quasi-Eulerian methods, in which the material properties are swept in the direction of the motion instead of mesh movement, can also be used. Shatoff et al. [10] moved the rail conductivity backward to simulate the armature relative motion, and current distribution at the armature is targeted. Hundertmark and Roch [11] used time-and position-dependent conductivity functions in the air-armature region to simulate the transients due to the armature movement. In the previous study [12], we obtained a lumped VSE resistance through a similar conductivity function to enhance the simulation model's accuracy and verified with experiments.
In this work, we take a step further to studies about quasi-transient FEM in the literature. The quasi-transient model is used to calculate the accuracy of the 2-D models in VSE calculations. VSE resistances are obtained for armature speeds. The study is conducted within a 2-D quasi-transient FEM and verified with 2-D transient FEM, where the armature moves as a domain. Then, the VSE resistance is calculated in 3-D with quasi-transient FEM. The fitted VSE resistances are added to the 3-D static FEM simulations, such as in [12] and [13], individually as a figure of merit. The constructed model is compared with experimental results to evaluate the performances of VSE models. Both EMFY-1 and EMFY-2 launchers, which are shown in Figs. 1 and 2, respectively, are used in this study. The geometric parameters of those launchers are given in Table I. In Section II, the procedure to obtain VSE resistance and the simulation methods are explained. In Section III, a 3-D hybrid simulation model that uses fitted VSE resistances to examine their performances is explained. In Section IV, the simulation results and experimental data are compared, and the accuracy of the 2-D approximation in VSE analysis is discussed. 0093-3813 © 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.

II. METHODOLOGY
The block diagram of the conducted work is presented in Fig. 3. Three different simulation models are developed, which will be explained in Sections II-A and II-B in detail. The main aspects of the simulation models are given in Table II. In the analysis part, velocity-dependent VSE resistances are fitted with given velocity samples v, using two different models  simulation-test accuracy [13]. The described model will be explained in Section II-B. In this study, the hybrid model's objective is to increase simulation-experiment accuracy and create a figure of merit that is used to compare performances of VSE resistances obtained by 2-D and 3-D analyses.

A. Analysis
Developing a realistic simulation model is critical for both the design and tactical operations of the EML. However, since the launcher operates in extreme conditions both electrically and mechanically, simulations are challenging. In this section, two different simulation models that are used in the analysis part are explained. The 2-D transient FEM and 3-D quasi-transient models are compared in terms of the calculated VSE resistance obtained from various armature speeds. In the 2-D transient FEM model, the armature is modeled with moving mesh elements. Therefore, there are discontinuous mesh nodes in the rail-armature interface different from static models. These mesh nodes are represented in Fig. 4(a). The continuity condition is applied in the interface to interpolate the field variables. Magnetic insulation is settled as the continuity condition. In Fig. 5, the current distribution in rails and the armature are given for steady and moving armatures. The current narrowing phenomenon is observed in the 2-D transient simulation model. The procedure used for modeling VSE resistance has two parts. First, 2-D transient FEM simulations are made at various constant speeds where dc excitation is used. Therefore, only dc and VSE resistances are calculated. Second, the evaluated VSE resistance is separated into two parts, transient and steady parts, which are used for fitting functions [13].
Proper use of the velocity feature in electromagnetic analysis needs deep physical insight. When the moving domain is of bounded extent in the direction of the motion, material properties vary in this direction [16]. It should be noted that the mesh movement at the sliding contact is not an option in 3-D. The underlying issue is that curl elements, which the interface uses for the magnetic vector potential being solved, give significant interpolation errors if the mesh is discontinuous. Thus, in the 3-D quasi-transient FEM, the armature mesh elements are not moved or deformed. However, the air-armature region's electrical conductivity is updated with time in the armature motion direction to emulate the movement. This approach can be called an Eulerian method since the material properties are modified to simulate the movement [10]. It should be   noted that, since the mesh elements are simply connected, field variables are continuous in all boundaries; thus, there is no need to introduce continuity conditions. In Fig. 6, the current distribution in the rail, which is calculated from 2-D transient and 2-D quasi-transient models, is given. It shows that the current distribution in the quasi-transient method is almost the same as the transient method in 2-D. The time and position-dependent conductivity function is given in Fig. 7. The armature is modeled as a rectangular prism for simplicity. Thus, the only spatial variation in the armature's electrical conductivity is in the direction of the armature movement.
Since all simulation models are time-dependent, proper time discretization is required. Variable-step solvers modify the step size during the simulation. They decrease the step size to increase accuracy when a model's states change quickly and increase the step size to avoid taking unnecessary steps when the model's states are changing gradually. Thus, fixed-step solvers have longer simulation times, but they are on the safe side. In this study, we used a variable-step solver with a minimum time-step penalty. Time-dependent solver convergence settled a value calculated from armature velocity to ensure the sudden conductivity jump. For example, for 1-mm maximum mesh size in the air-armature domain and 1000-m/s armature speed, time-dependent solver convergence is allowed to have 0.5-μs step size as maximum to prevent discontinuity due to movement. Such convergence limitation increases simulation time, but it enables us to see the transient due to the armature movement.
2-D models are built on the fact that field variables do not change in the symmetry direction. Thus, in the VSE investigation of EMLs, it means that the rails' current distribution does not change through rail height. However, the validity of the assumption is not determined in the literature. In Fig. 8, the current distribution in the middle surface of the rail when the armature has 1200-m/s speed is shown. Arrows denoted with letters used to show cross section positions are used in Fig. 9. The rail's middle surface is demonstrated with black lines in the rail cross sections in the same figure. An important argument that can be extracted from Fig. 9 is the current distribution is not homogeneous through rail height, which is also a sure sign that the error comes from the 2-D approximation of R vse . The current distribution, which is calculated from 2-D transient/quasi-transient FEM models, is similar to the current distribution at the rail's middle surface only. The outcomes obtained from the analysis part are listed as follows.
1) The verification of the quasi-transient model is made in 2-D. The current distributions and the resistance calculations are identical.
2) The current narrowing phenomenon due to the existence of VSE occurs mostly on the middle surface of the rails.

B. Validation
The validation section's primary goal is to find a figure of merit for VSE resistance calculation methods. For that purpose, the fitted VSE resistances from the previous section are embedded in the 3D-static FEM coupled with PPS circuitry, which is called the hybrid FEM model. Then, several experiments are conducted to explore the precision of the simulations. Before going into the comparison study results, a brief explanation of the hybrid FEM model is given.
EMLs have a high aspect ratio due to long rails, resulting in a large number of mesh elements. In order to reduce  computational load, partial rail modeling can be used. The rail portions that are not modeled in the FEM model can be modeled as lumped circuit parameters. Ceylan et al. [13] used a hybrid model for simulation of EMFY-1 EML. One of the geometries that we used in the 3-D hybrid FEM model is shown in Fig. 10. This geometry belongs to the EMFY-1 launcher, which has a 25 × 25 mm caliber. The armature is modeled stationary at preload distance, 0.7 m. The structure of the hybrid model, which we used, is demonstrated in Fig. 11.
It should be noted that the hybrid FEM model has a stationary armature. Therefore, the current distribution of the rails did not narrow due to VSE. The current distributions in the rail when 3-D static and 3-D quasi-transient methods are used are illustrated in Fig. 12. The physical effects, such as inductance and resistance variations due to armature movement, are added to the simulation environment using lumped parameters. One of the mentioned physical effects is VSE resistance, which is denoted as R vse , which is calculated from 2-D transient and 3-D quasi-transient analyses.
Other lumped parameters that we used to model the rail portion which related with armature regions are defined as exterior inductance L ext , interior inductance L int , ac resistance R ac , and back EMF resistance R emf . L ext is used to model inductance contribution in the air region between rails. However, L int is used to model inductance contribution in the rails and air region between rails, excluding the exterior region. Total inductance contribution due to p amount armature movement can be modeled, to sum up, that inductance (1) The back EMF resistance, R back EMF , formulation can be extracted from (2) to (5). EML can be considered as single turn coil. Therefore, linkage flux λ(t) can be calculated from the inductance and rail current I (t) Fig. 10. Geometry used in 3-D FEM analysis for EMFY-1. The integration area that was used in the inductance calculation is shown as purple. The area between rails is called exterior boundary and the inductance, which is related in this area called L ext . The area inside the rail is called the interior boundary. The interior boundary is used in rail inductance calculation, which is called L int . Fig. 11. Hybrid FEM model used to simulate EMLs. Resistances and inductances between the PPS network and 3-D FEM body use to model the armature movement. With that structure, not only the computational complexity is reduced but also the total number of mesh numbers is reduced as all the rail length is not required to be meshed.
Back EMF due to armature movement is modeled as a resistance. Although it behaves as a voltage source, this assumption, from a mathematical perspective, does not pose a problem. R ac can be calculated using the joule loss in the rail portion volume V. In (6), the resistivity of the rails is denoted as ρ rail . J rail and . represent the current density inside the rails and l 2 -norm, respectively The only lumped parameter that cannot be computed directly from the hybrid FEM model is R vse . This parameter is imported from the studies that were done with 2-D transient and 3-D quasi-transient methods. Since all the other parameters are kept constant, the error that comes from R vse calculations is directly transferred to the outputs of the hybrid FEM simulations. We avoided full discussion of the 3-D Hybrid FEM model since we do not want to duplicate  with previous works. Interested readers can access detailed discussion about other parameters in [12] and [13].

III. RESULTS
In Fig. 13 Fig. 14. It should be noted that R vse values are computed between 0.2-and 0.5-m ranges. The reason is that rail portioning is used in both models for reducing computational load. Since they showed a linear relationship with rail position, we extrapolate them for other ranges. After exploiting lumped parameters, including R vse , which obtained 3-D quasi-transient and 2-D transient FEM results, the hybrid simulation model is tested with experimental results to investigate the differences. The rail currents are compared to evaluate accuracy of R vse fitting since it is a significant contributor [13] to the total EML resistance.
The open area launch tests of EMFY-1 were conducted in 2018 [13]. EMFY-2 was tested indoors in 2019. Experiment parameters of EMFY-1 and EMFY-2 are given in Table III  show that the 3-D quasi-transient method is the most consistent with experimental results.
The overestimation of R vse with the 2-D method can be explained using the experimental results. At the later stages of the launching process, the EML circuit can be modeled as a passive RL circuit. The reason is the fact that, after the discharging of the last PPS unit, all free-wheeling diodes, D 1 to D N , get into conduction mode. Thus, capacitors are disconnected from the EML, as shown in Fig. 19. The rail current decay rate depends only on the EML inductance, L EML , and resistance, R EML . The 2-D transient simulations have faster decay than experimental results due to the overestimation of R EML . It should be noted that discrepancy at the late stage of   the experiments links up with R vse dominance at the late stage of the launch [13].

IV. DISCUSSION
Although the quasi-transient method can model velocity-dependent resistance calculation in EMLs, there  are some limitations of the proposed method. One of the disadvantages of the quasi-transient method is that the armature cannot be analyzed in detail due to geometry limitation due to the conductivity function definition. In this study, the armature is modeled as a rectangular prism, as shown in Fig. 20, but they have more complex geometries in actual cases. Another disadvantage is that the time-dependent solver convergence is limited when the quasi-transient method is used. To prevent sudden conductivity jumps between mesh elements, time-steps are locked at a specific value, which increases the simulation time. 2) The hybrid FEM model simplifies the computational problems of EMLs, such as the number of mesh elements and the requirement of the armature motion in the simulation. Although the field distribution at the sliding contact does not contain the current narrowing phenomenon, the VSE's resistive effects are added to the simulation environment, which improves losses and efficiency calculations.
3) The 2-D model overestimates VSE resistance. The reason for that is the current narrowing phenomena occur in the middle of the rail more dominantly. 2-D approximation assumes that the current distribution is identical in all rail height directions. However, this is not the case in the actual rail current distribution. The exploited VSE resistance is a critical lumped parameter for the overall system simulations. The simulation-experiment coherence is enhanced with the proposed method. The same procedure is applicable with different caliber EMLs or actuators, which has VSE on their conductors.