Three-dimensional numerical simulation of multi-physical coupled environment during electromagnetic propulsion

. Electromagnetic propulsion technology is a new propulsion technology which uses electromagnetic force to push objects into high or ultra-high speed. However, during propulsion, the armature and launch load are affected by harsh multi-physics environment. The coupled effect was often ignored or over idealized in previous numerical simulations, which leads to large errors between simulation and experiment. In this paper, a three-dimensional multi-physics coupling environment simulation model is established. The distribution of electromagnetic field, structure field and temperature field in armature from static to high speed are effectively calculated. Moreover, the coupled effect between physical fields and the influence of key parameters on the simulation results are revealed. In conclusion, this model can reproduce the physical field distribution of armature and rail during electromagnetic propulsion. The time-varying resistance and inductance and the velocity skin effect (VSE) are the key parameters affecting armature movement. This method reveals a feasible path for the simplification of multi-physical model and the mitigation of extreme environment.


Introduction
As a new type of launching platform, electromagnetic railgun has the characteristics of high initial velocity, low pollution and large muzzle kinetic energy compared with conventional platform [1][2][3]. It can thrust the carrier into 2000 m/s or higher speed in a few milliseconds, which can easily break through the speed limit of traditional chemical energy platform [4]. However, during the launch time, the in-bore harsh multi-physical field environment is one of the factors restricting its development [5]. The electromagnetic force is produced by the interaction of in-bore strong magnetic field and current flowing through armature, which promotes the high speed movement of armature and launch load [6]. The transient electromagnetic field and force field produce Joule heat and friction heat respectively. It will affect the conductivity of conductors, which affect the contact status between rails and armature [7]. Therefore, it is of great importance to study the multi-physics coupling effects and analysis method for better designing the launch load system.
The sliding electrical contact between armature and rails during launch has always been the main problem restricting the calculation of three-dimensional transient electromagnetic force of electromagnetic launcher [8]. Hsieh et al. [9] developed a special simulation software EMAP3D for electromagnetic propulsion process research. The magnetic diffusion process of moving conductor was described by Lagrange equation, and the coupling analysis of electromagnetic field, force field and temperature field was realized. Based on the hybrid finite element method and boundary element method, Lin Q. H. et al. [10] established the simulation model of electromagnetic field and temperature field, and realized the multi-physical field coupling analysis

Theoretical model of multi-physical field environment during electromagnetic propulsion
The multi-physical field environment during propulsion is closely related to the movement of the armature. In this section, the equivalent circuit model is established to simulate the pulse current at each time. During launch, the pulse current passes through the rails and armature and generates the magnetic field. The interaction between magnetic field and current produces electromagnetic force. The force field includes the transient stress field and the transient dynamic field. It is assumed that the launch load is stationary relative to the armature, and only the transient dynamic field is considered. For the temperature field, it is mainly caused by the joule heat and the friction heat. This section will analyze the quasistatic physical model, and write programs to achieve Iterative multi-physical field coupling model.

The input current model
Nowadays, the energy supply system of electromagnetic railgun is mainly consisting of multiple capacitor banks. The equivalent circuit model of electromagnetic launch system is shown in Fig. 1. The capacitor bank is composed of numbers of capacitor energy storage modules. One or more capacitor banks are connected to the confluence plate at the tail of the electromagnetic launcher. The rails can be equivalent to inductance and resistance varying with time, the armature can be simplified to resistance . The inductance and resistance in the whole discharging circuit are composed of capacitor group and electromagnetic launcher. The capacitor group is composed of multiple encapsulated sub-circuit modules, and each module is composed of different number of capacitor modules. The appropriate input current is realized by controlling the trigger time of different sub-circuit switches. The inductance and resistance in the electromagnetic launcher will change with the movement of the armature, thus affecting the current at the subsequent time [20].
During the RLC circuit discharge, = − ⋅ / represents the input current, = = − / represents the voltage at both sides of the resistance, = / = − / represents the voltage at both sides of the inductance. According to Kirchhoff 's law, there are: By solving capacitor voltage and loop current, the time of maximum current is = = ⁄ . So the discharge current of single capacitor bank is as follows: where is the initial charging voltage of capacitor, and is the total resistance in the circuit, mainly including the internal resistance of the line and the real-time resistance of rails and armature during movement. The resistance is related to the moving distance of the armature, and the resistance gradient of the rails is usually taken as 0.1 mΩ/m [21]. is the total inductance in the circuit, including load inductance and real-time inductance. The inductance is also related with the movement of the armature. Generally, the rails' inductance gradient are 0.51 μH/m [22]. Therefore, the equation of the discharging multiple capacitor banks is: (3)

Three-dimensional magnetic diffusion model of electromagnetic propulsion process
As the wavelength of input current is larger than the size of the railgun. The displacement current can be ignored. It is studied as a quasi-static system. According to the differential form of Maxwell equation and the Ohm's law: In three-dimensional calculation, it is assumed that there is no ferromagnetic conductor in the rail and armature, and the conductivity and permeability are independent of time , then is used instead of , and the magnetic diffusion equation of the conductor in the three-dimensional model is: where ⋅ / is the velocity diffusion term. Assuming that the armature is stationary and the rail has -axis reverse velocity relative to the armature, the magnetic diffusion equation in the armature region can be expressed as: The magnetic diffusion equation of armature and rail is solved by setting boundary conditions and initial value , and the current density values in -, -and -axis are obtained as follows: Ignoring the magnetic induction intensity of the rails and armature in the -, -axis, the current density in -, -and -axis are: Since the input current is a low-frequency current, the generated magnetic field is similar to a stable magnetic field. Assuming that the package of the projectile is made of non-metallic materials, the magnetic induction at the inspection point of the launch load is calculated by using Biot-Savart's law: where is the current distribution area, is the vector diameter of the source point (current element ), is the vector diameter of the field point. Therefore, the magnetic field distribution in front of the armature can be calculated according to the vector diameter of the source point and where ℎ is the height of the rail, is the skin depth, is the initial position of rail, is the axial length of armature. The magnetic flux density generated from the source point of each surface of the armature to the central axis is: where ℎ is the height of the armature, is the skin depth, is the initial position of rail, is the axial length of armature.

Three dimensional dynamic model of armature
The force of armature can be expressed according to Newton's second law: where is the displacement of armature, is the axial electromagnetic force on the armature, is the friction force on the armature. According to the electromagnetic propulsion equation, the electromagnetic force on the armature can be expressed as: The axial electromagnetic force can be expressed as: The friction is directly proportional to the pressure of the armature on the rails: In the structural field, the electromagnetic propulsion received by the armature is essentially provided by the Lorentz force, which is mainly divided into the axial propulsion of the armature and the component force to the rail :

Three dimensional thermal diffusion model
According to Fourier's law and thermal equilibrium law, the three dimensional thermal diffusion equation is: The first item of the right type is the heat generated by Joule heat, and the second item is the heat generated by friction heat.
The joule heat of rails and armature is: The friction heat is:

Multi-physical field coupling simulation model
The electromagnetic propulsion process is divided into a large number of microseconds in sequence. In a single time period, the circuit simulation module is used to assign capacitor bank parameters and trigger threshold settings. Then the input current is introduced into the quasi-static electromagnetic model to calculate the electromagnetic field. The variation of magnetic flux density and current density is simulated, and the results are introduced into the mechanical model to update the velocity and displacement. Finally, the results of electromagnetic model and mechanical model are introduced into the thermal model to update the temperature distribution. The last time solution of the previous time period is taken as the initial value of the next time period for iteration to realize the simulation of the whole propulsion process. According to the multi-physical field simulation model, the current distribution, magnetic field distribution, temperature distribution and velocity can be obtained.

Multi-physical coupling simulation of electromagnetic propulsion process
In this section, the simulation model is established based on the magnetic field module, PDE module, structural mechanics model, solid heat transfer module and circuit module in the commercial software COMSOL Multiphysics. On this basis, the secondary development is carried out, and the multi-module coupling operation is realized by programming.

Simulation of capacitor discharge circuit
According to the pulse power structure in Fig. 1, the circuit simulation module is used to build the simulation circuit. The capacitance value is 2 mF, the initial charging voltage is 9.5 kV, the wave-modulation inductance is 48 μH, the line internal resistance is 10 mΩ, the load equivalent inductance is 2 μH, and the load equivalent resistance is 20 mΩ. The package sub-circuit modules P01 to P04 represent four capacitor banks. The number of pulse power modules and trigger time of each capacitor bank are shown in Table 1. The charging voltage of each PFU module is 9.5 kV. The equivalent resistance of each transmission cable is 7 mΩ, and the equivalent inductance is 0.6 μH. The resistance and inductance values in armature and rail are determined by the calculation of structural mechanics and PDE module in the following part. The total displacement of armature position is obtained by adding the calculated displacement and the initial displacement. The resistance gradient of the rail is 0.1 mΩ/m, and the inductance gradient is 0.51 μH/m for calculation. The voltage drop caused by the increment of inductance and resistance is calculated and introduced into the main circuit to output the current value at the next stage. With the iteration, the current distribution of electromagnetic propulsion is obtained.

Multi-physics coupling simulation
The parameters of rail and armature are constructed according to the model in [23]. The height of rail is 124 mm, the width is 85 mm, and the length is 12000 mm. The commonly used U-shaped armature is selected for armature, and the total width is 100 mm. The 1/2 model in Fig. 3 is 50 mm, and the contact area between armature and rail is 124 mm×100 mm. The rail and armature are divided by tetrahedral mesh. In order to reduce the calculation amount, the adaptive mesh is adopted. The mesh of armature area, the back side of rail and the contact area is densely set. Since the influence of front-end position is small, the mesh is sparse. The density, conductivity, relative permeability, thermal conductivity and specific heat capacity of rails and armature are set in Table 2. (2) In order to intuitively represent the division of boundary conditions, the boundary conditions are shown in the plane. Since the rail and armature are symmetrical about the center of -axis, the 1/2 model is used for analysis. When the armature slides at a high speed relative to the rail, the VSE will be caused, so that the magnetic field and current are concentrated on the inner edge of the guide rail and the rear edge of the armature. The boundary 1 are set as = 0. The boundary 2 is set as = , where the current line density is = /ℎ, is the input current, ℎ is the rail height. The left rail's boundary 3 is set as / = 0, the armature's boundary 3 is set as / = 0.

Analysis of numerical simulation results
The pulse current obtained by calculation is shown in Fig. 5(a). It can be found that the total current synthesized by the sequential discharge of four pulse capacitors, and the trigger setting time of each pulse valley is close to that of the capacitor bank. In the simulation, the peak current reaches 1.15 MA, and the duration is 12 ms. Pulse current mainly formed three steps, the first step is the rising period, the time is 0-0.1 ms, the current is mainly provided by the first capacitor bank, the rising trend is close to sine function. The second step is the platform period, the time is 0.1-0.9 ms, and the current is provided by all capacitor banks. The current curve in this step forms multiple pulse peaks due to the different trigger time. The last step is the drop period, which is synthesized by the discharge phases of multiple capacitor banks. . It can be seen that the current density is mainly distributed on the surface of the rail and the armature, the contact area between the rail and the armature. Due to the skin effect, the current density at the armature edge is large. With the increase of velocity, the peak current density is mainly distributed in the inner contour of the armature and the contact area between the armature and the rail, and the peak current density reaches about 1.2×10 10 A/m 2 , as shown in Fig. 5(c1) and Fig. 5(d1). Moreover, the peak current density is multiple relative to that in other regions, which is consistent with the simulation results in [24]. The reason for this phenomenon is the combined action of skin effect and velocity skin effect.
Figs. 5(b2) -(d2) represent the distribution of magnetic induction at each time. It can be found that the magnetic induction is mainly distributed in the inner surface of the rail and armature, and the peak value reaches 20 T. From the -plane section, there is a shallow diffusion depth near the armature. The VSE resulting in current concentration in the tail of the armature. As shown in Fig. 5(d), the armature velocity reaches 1680 m/s at 11 ms, and the magnetic induction in the inner surface of the armature gradually spreads forward with time. Fig. 5(b4)-(d4) show the temperature response of the armature and rail at three moments, and the calculation accuracy is sufficient to simulate the millimeter-level or even more precise movement in the initial propulsion process. At the initial time, the temperature is mainly distributed on the inner surface of the armature and the rail. The temperature mainly due to joule heat generated by the current density, and a hot spot is formed at the contact point between the armature tail and the rail, with the peak temperature reaching 440 K. Then the speed increases and the friction heat ratio increases gradually. The hot line is formed on the contact surface of the armature and the rail. The temperature reaches 450 K and the peak temperature of the hot spot reaches 700 K. In Fig. 5(d3), the temperature on the hotline is 600 K and the temperature on the hot spot is 1100 K. It can be seen that the initial supply of temperature during electromagnetic launch is mainly provided by Joule heat generated by large current. With the increase of velocity, the VSE appears. The heat energy is mainly concentrated on the contact point between armature tail and rail, and there is also a high temperature distribution on the contact surface. The axial and radial acceleration distribution obtained by Eq. (19) in Fig.6(a). The axial peak acceleration is 19000 g, and the radial peak acceleration is 8000 g. The radial pressure is expressed in the form of acceleration for comparison with the axial acceleration. The radial pressure is converted into friction force during propulsion, which affects the movement of the armature. The axial acceleration lasted about 8 ms in the peak stage, and the maximum difference was 6000 g. The velocity and displacement obtained after integration are shown in Fig. 6(b), For the whole propulsion process, the armature is accelerated to nearly 1700 m/s in just 12 ms, with the propulsion distance of only 11 m. Based on Eq. (9)(10)(11)(12)(13)(14), the magnetic field distribution of central axis position is solved, which provides a reference for the design and protection of precision circuit components. As shown in Fig. 7, the simulation results show that the magnetic induction decreases exponentially with the increase of the distance from the armature. The peak magnetic flux density at 10 mm from the armature reaches 1.4 T, while the magnetic flux density at 210 mm is only 0.1 T. Therefore, the position of the magnetic sensitive device is as far away from the armature as possible, and the direction is parallel to the magnetic field direction as far as possible, so that the influence of the pulse magnetic field on the device is minimized.

Comparison and discussion
The main influencing factors in the complex environment are the velocity skin effect and the changes of resistance and inductance during movement. The elaboration and demonstration of their influence are discussed in the following parts.

Influence analysis of input current
The change of resistance and inductance caused by the moving armature during launch. The actual current waveform is formed by controlling the sequential discharge of multiple capacitor banks. The equivalent current method is used to simulate the current during launch [25], and the influence of time-varying circuit parameters of the launch device is ignored. The equation of the current can be seen in Eq. (23): In this section, the equivalent resistance and inductance is changed with the velocity and displacement of the armature, and the coupling of the structural field and the circuit simulation model is realized. The simulation results are compared with three-step current and those without considering the transient changes of resistance and inductance, as shown in Fig. 8.
The method in this paper (time sequence current) is compared with two common equivalent input current methods. One is the equivalent current (without considering the moving armature), and the other is simplified three-step current. The comparison between the current curve obtained in this paper and the other two methods is shown in the Fig. 8.
At the rising step of the current, the variation characteristics of the three input currents are basically the same. In the platform step, the difference between the three-stage current and the method in this paper is the largest. It can be found that the error between the peak value and valley value reaches about 0.2 MA. The capacitor bank discharge sequence generated by the ripple phenomenon was ignored with the method of three-step current. The change of resistance and inductance was ignored with the method of the equivalent current, and the difference is also existed. The amplitude at the initial step is slightly larger, and is slightly smaller at the later step. The dynamic characteristics calculated by using the three methods are shown in the Fig. 9. It can be found that the acceleration value calculated by three-step current cannot show the ripple phenomenon in sequential current calculation. The maximum difference in acceleration is 5000 g in 1-9 ms, while the result of three-step current is a stable waveform. The equivalent current has little difference from the time sequence current. The value in 1-4 ms is larger, and the value in 4-12 ms is smaller, with an error of 1000 g. The comparison of velocity is shown in Fig. 9(b). The velocity obtained by the three-stage current is the largest different from that in this paper, the value is up to 260 m/s, and the value between equivalent current and that in this paper is up to 50 m/s.

Influence of VSE on simulation results
The movement of the armature relative to the rails cause the VSE, which takes the magnetic line out of the contact area, thus counteracting the slow magnetic diffusion into the conductive material of the rail and armature. In this paper, the partial differential equations of rail and armature are established to simulate the relative movement, and the results of each step are iterated to the next step to realize coupling. The velocity skin effect is the main factor affecting the current distribution during launch. In the previous discussion, it can be found that the magnetic field, current and temperature distribution of high-speed motion are significantly different from those of low-speed and static processes.  Fig. 11. The distribution of current density and magnetic induction in moving state As shown in Fig. 10, the magnetic flux density distribution at the -interface of the contact area between the rail and the armature, and the current density distribution between the armature and the rail. The distribution of magnetic flux density and current density under the influence of velocity, as shown in Fig. 11. It can be found that in the static process, the diffusion depth of magnetic flux density is 14 mm, and gradually decay. During the movement, the flux density concentrates on the surface and the diffusion depth is only 2.5 mm. The static current density is mainly distributed in the rail and the armature edge, and reaches the highest at the boundary position of the armature throat. The distribution of current density in the motion process is significantly different from that in the static process. When the VSE is not considered, the current density is distributed on the rail. The peak current density is affected by the skin effect, and the peak current density is concentrated in the armature and the rail corner, reaching 7×10 9 A/m 2 . When considering the VSE, the current density on the rail surface increases to 7×10 9 A/m 2 , and the peak current density at the contact position between the rail and the armature reaches 9×10 9 A/m 2 . It can be found that the velocity skin effect will lead to the concentrated distribution of current density and the surface of the rail and the contact position with the armature, and form a concentrated distribution. It mainly affects the current density distribution on the rail and the current density distribution at the connection position between the rail and the armature.
The dynamic characteristics of armature considering and not considering VSE are compared in Fig. 12. It can be found from Fig. 12(a) that the -axis acceleration distribution obtained by the two simulation models is basically coincident within 0.002 s. With the increased of the velocity, the distribution of current density and magnetic induction are affected by VSE. Thus, the electromagnetic force is changed. The error between two simulation models is gradually increased in 0.004-0.012 s, and the peak error reaches about 4000 g. In contrast, the error of -axis acceleration obtained by the two simulation methods is relatively low, indicating that VSE has little effect on axial pressure. The comparison of velocity between the two simulation methods is shown in Fig. 12(b). The maximum velocity error reaches about 150 m/s. Therefore, the VSE should be considered in the simulation. a) Acceleration b) Velocity Fig. 12. Dynamic characteristics of armature considering and not considering VSE

Experimental verification
In order to verify the reliability of the simulation method proposed in this paper, we compare the simulation and experimental results in [26].
The influence of VSE is simulated by changing the rails' structure model and the contact area between rails and armature in the reference. Moreover, the electromagnetic field and structural field environment during launch are measured by the projectile-borne storage testing system. The method in the reference did not realize the magneto-mechanical coupled effect during launch. The electromagnetic field and structural field were not calculated iteratively, and the influence of fieldcircuit coupling was not considered. a) Velocity b) Displacement Fig. 13. Comparison of velocity and displacement between method in reference, this paper and experiment Based on the theoretical model and numerical simulation method in this paper, the input current, structure and material parameters of armature and rails in the reference are established. The differences between the two simulation and experiment results are compared and analyzed in Fig. 13 and Fig. 14. It can be found that the simulation results in this paper cannot completely reproduce the experiment results, as some simplified factors in the simulation and the non-ideal factors in the experiment. However, compared with the simulation results in the reference, the error between the simulation and the experiment results becomes smaller. The error between the simulation of this paper and the experiment is reduced to 100 m/s in velocity and 0. 5  displacement. This also reflects the reliability of the method in this paper. The comparison of the magnetic field distribution during launch is shown in Fig. 14. The distance between the sensors and the armature is 60 mm. According to the parameters in this paper, the magnetic induction is calculated by the simulation method in this paper, as shown in the black curve in Fig. 14. Considering the influence of ferromagnetic shell material, there is no ripple phenomenon caused by pulse timing current, and the simulation results in this paper are closer to the experiment results. The error of peak magnetic induction is reduced to 0.6 T.

Conclusions
In this paper, the mathematical model of three-dimensional electromagnetic field, structural field and temperature field in electromagnetic launch process is established. The secondary development of finite element software is carried out to restore the current distribution, armature motion characteristics and magnetic field and temperature distribution in the electromagnetic launch process. The influence of VSE on magnetic field diffusion and current density distribution is revealed, and the magnetic field distribution in front of armature is deduced based on Biot-Savart law. The contrastive analysis of the simulation results is given and the conclusions are summarized as follows: 1) The iterative calculation of circuit and multi-physical fields is realized by using the method of time sequence circuit model. The measured value or equivalent value is no longer required for the input current. Compared with the method of equivalent current and three-step current, the higher accuracy is obtained in this paper.
2) The VSE changes the distribution of current density and magnetic induction between the armature and the rail. This phenomenon leads to the concentrated distribution of current density and magnetic induction on the armature tail. Moreover, the dynamic characteristics are also influenced, and become more obvious with the increase of velocity.
3) The thermal field is mainly contributed by joule heat and friction heat. The VSE forms a hot spot in the contact area of armature tail, and the temperature is as high as 1200 K. Meanwhile, the effect of friction heat cannot be ignored. With the increase of the velocity, the hot line between the rail and the armature becomes more obvious.