Performance Evaluation of Dampers Under Dynamic Impact Based on Magnetic Field FE and CFD

It is necessary to install a damping energy absorbing structure for a certain type of naval gun test shell without projectile belt. Applicable high-performance dampers currently include viscoelastic colloidal damper and magnetorheological damper, and two kinds of special dampers for test shells have been designed and manufactured. In order to evaluate the mechanical properties of these two dampers under dynamic impact, a one-way coupling numerical method was proposed for simulation and analysis. Quasi-static experiments were used to verify the model based on magnetic field finite element (FE) and computational fluid dynamics (CFD). The model is also characterized by passive dynamic mesh and user-defined functions (UDF) based on C language. The results showed that the CFD models of both dampers were accurate. Though the colloidal damper is simpler in structure and lower in the cost, the MR damper has more significant damping and higher energy absorption ratio under impact. Compared with colloidal damper, MR damper is characterized by long working time, multiple vortices, chaotic streamlines and more uniform temperature distribution. When the power was continuously applied, the velocity of the impacted MR damper dropped to zero within 80ms, and the piston finally stopped at a stroke of 1.3mm. The colloidal damper reset within 8ms, and the reset velocity was 1.1m/s. In this study, the performance evaluation results of the two dampers were obtained, and the latest numerical methods for the mechanical performance analysis of colloidal dampers and MR dampers under dynamic impact were provided.


INDEX TERMS
Liquid density exponent n 2 Rheological index n 3 Gas polytropic exponent σ Electrical conductivity VOLUME  Magnetorheological (MR) dampers are usually filled with special magnetorheological liquid made of non-conductive liquid, tiny ferromagnetic powder and a small number of additives [1]. This special fluid shows the characteristics of common Newtonian fluid without the presence of magnetic field, however, when a magnetic field exists, the tiny metal powders form chains along the magnetic induction line, which will cause the fluid to switch to the characteristics of viscoplastic non-Newtonian fluid [2]. The fluid yield stress is determined by the magnetic flux density, which explains why the damping force of MR damper is controllable [3]. Compared with colloidal dampers and other dampers, MR dampers feature larger hysteresis loop range [4]. Currently, MR dampers are being developed rapidly and widely used in vehicle engineering [5], medical equipment [6], aircraft [7] and weapon industries [8].
Colloid dampers are usually filled with highly viscous and elastic colloidal fluid mainly composed of macromolecule silicon-oxygen organic polymer, such as silicone oil [9]. The material features excellent compressibility, chemical stability and damping viscosity due to its unique molecular helix structure and many molecular side groups. Currently, colloidal dampers are also used in a variety of industrial fields, including CRH trains [10], automobiles [11], light weapon [12], artillery [13], etc.
The geometric dimension and mass characteristics of test shell are almost the same as those of live shell, except that the belt diameter is intentionally reduced to be lower than that of the rifling. However, without the deformation and friction of the belt, excessive rebound velocity will cause the failure of closing breechblock, which may result in the failure of live shell simulation. Therefore, the damping structure inside test shell is necessary [14]. Fig.1 shows the positional relationship among test shell, barrel and damper.
The basis of studying MR dampers is the constitutive relationship of magnetorheological fluid, that is, the theoretical viscoplastic equation. In order to describe this rheological property, many viscoplastic equations have been proposed and used. Pipit et al. [15], Tahsin and Sevki [16], Yu et al. [17] and Chooi and Oyadiji [18] used the Bingham model and Herschel-Bulkley model to describe viscoplastic fluids in their research. Papanastasiou [19], Papanastasiou et al. [20] used an exponential function to describe the two-dimensional steady flow of viscoplastic fluid. However, the latest equations at present were composed of hyperbolic tangent functions and used in many simulations and calculations [21]- [24]. The studying colloid damper is the damping force model. The damping force proposed by Gokhan et al. [25] is a function of velocity, while that proposed by Gloria [26] is a polynomial function of velocity. Ahmet [27] established Maxwell model, Kelvin model, Oscillation model of thermoplastic polyurethane elastomeric buffer respectively, and used the experimental results for comparison. Ding and Ou [28], Liu [29] and others analyzed the flow of fluid through gaps and holes, and proposed that the damping force was a function related to velocity, geometric parameters, and material parameters. Jia [30] and Wang [31] analyzed the fluid in the gap based on the idea of differentiation, and also proposed a function related to velocity.
With the development of computers, finite element analysis and computational fluid dynamics simulation techniques are increasingly used for the analysis of MR fluids and MR dampers. Compared with the traditional one-dimensional analytical model, numerical simulation is more accurate and available for better visualized parameters. For example, parameter contours, streamlines, non-Newtonian viscosity region, vortex distribution, etc., which are all of great significance to the research of MR damper [32]. However, there are still many difficulties ahead in CFD simulation at present: (a) The problem of passive motion with implicit boundary conditions. (b) The definition of non-Newtonian action region following the moving mesh. (c) Coupling magnetic field FEA and CFD. (d) Treatment of viscoplastic constitutive equation.
The multiphysics model was established by David et al. [33]. The magnetic field and fluid were coupled in COMSOL, besides, experiments and simulations were compared under sinusoidal periodic motion. Zekeriya and Tahsin [34] coupled the magnetic field and fluid in ANSYS/ CFX, and compared the experiments and simulations of uniform quasi-static motion and sinusoidal periodic motion respectively. The non-Newtonian viscosity was written in a special language (CCL). Wael et al. [32] coupled the magnetic field and fluid in COMSOL and ANSYS/Fluent, and also performed sinusoidal periodic motion analysis. Tong et al. [35] established a fluid mesh model for colloidal dampers with holes and without gaps. Sun et al. [36] established a laminar CFD model with double ended and gap type damper, and compared it with the experimental data of the drop hammer impact.
Di et al. [37] established a three-dimensional fluid-solid coupling model of the colloid recoil device and verified it by firing experiments.
From the above literatures on the two kinds of dampers, it could be concluded that the early researches were mostly based on simple one-dimensional analytical models which were gradually replaced by more accurate numerical simulation methods. However, so far, in the researches of magnetic field FEA and CFD, only sinusoidal periodic motion or uniform quasi-static motion has been applied. These are all active motion problems with ''explicit'' boundary conditions. While the problem to be studied in this paper is the simulation of MR damper under a certain velocity impact, which belongs to the ''implicit'' passive motion problem with unknown boundary conditions in advance. No relevant reports on such problems have been found in the literature, therefore, this study is expected to fill the gaps in the literature.
Among them, as a complex problem of multi-physics coupling and interdisciplinary, the MR damper is characterized by electromagnetic field generation, magnetorheological fluid properties change, air spring action, fluid-solid coupling, etc. A one-way coupling numerical simulation method was proposed in this paper to analyze and compare the mechanical properties of the two dampers. Firstly, the finite element model was established in COMSOL/ Multiphysics and the magnetic circuit was calculated. Then the magnetic field analysis results were input into ANSYS / Fluent, and the damper mechanical characteristics under dynamic impact was obtained by CFD calculation.
In this study, innovative numerical method combining magnetic field FE and CFD was proposed to solve the passive motion problem of MR damper. For the first time, the mechanical properties of colloidal damper and MR damper with the same size were evaluated through experiment and simulation.

II. THE NUMERICAL APPROACH A. STRUCTURE AND GEOMETRY
The main external dimensions of the colloidal damper and the MR damper have to be close to make sure that they could be installed in the same test shell. Fig.2 shows the structural principle for colloidal damper and MR damper. The colloidal damper is mainly composed of piston rod, piston, cylinder, limit block and viscoelastic fluid material. Due to the inconsistent dimensions at both ends of the piston rod, the internal volume of cylinder becomes smaller when the piston and piston rod move towards the compression direction. The colloidal fluid is compressed, and in this way the elastic potential energy is stored. At the same time, the colloidal fluid flows through the annular gap to generate a damping force. When the piston and piston rod move towards the recovery direction, the elastic potential energy is released; besides, the colloidal fluid flows through the annular gap in the opposite direction, continuing to dissipate energy. On the other hand, MR damper is mainly composed of piston rod, magnetic fluid piston, air piston, cylinder, limit block, coil, compressed air and  magnetorheological fluid. The diameter of both ends of piston rod is equal, and the elastic restoring force mainly comes from compressed air. Besides, there are two sets of coils in the piston, which are sealed by the casing. The wire passes through the inside of piston rod and is connected to an external DC power supply. When the coils are energized, the magnetic field distribution is approximately the same as that shown in the figure. As the result, the magnetic field causes the ferromagnetic metal particles in magnetorheological liquid to form chains along the magnetic induction line, which drives part of the fluid in the annular gap to show a great yield stress. The fluid passing through the gap will generate a large damping force. The dimensions of the two dampers are shown in Table 1.
The installation diagram of colloidal damper and MR damper is shown in Fig. 3. The damper cylinder is fixed on projectile by limiting structure and stop screw, and the damper piston rod is fixed on cartridge by stop pin. When the relative velocity and displacement between projectile and cartridge are generated, the damper starts to work. The initial fluid pressure of the colloidal damper and the initial pressure of the MR damper air spring can assist them to reset. From the structural point of view, the obvious advantage of colloidal damper lies in its simple structure. The air spring, DC power supply and control system of MR damper will result in the increase of manufacturing cost and use cost.

B. PROPERTIES OF FLUID MATERIALS
Defined as a viscoelastic fluid, the colloidal fluid material is a polyorganosiloxane with a chain structure of different degrees of polymerization. It shows both the viscosity of liquid and the elasticity of solid, and the viscosity is much higher than that of general fluids.
As shown in (1), liquid compressibility is described by the simplified Tait compressible liquid state equation [39], and the bulk modulus of elasticity is used to measure the compressibility of a liquid. The smaller the bulk modulus of elasticity is, the easier the liquid is to be compressed.
where, n 1 is the liquid density exponent, E 0 is the reference bulk modulus, ρ is the liquid density, ρ 0 is the liquid reference density, and p is the liquid pressure. The constitutive equation of viscosity for generalized non-Newtonian fluids can be described by the Ostwald de Waele empirical formula, which is a power law model [38]. The model considers viscosity to be a power function of shear rate, as shown in (2).
where, µ is the dynamic viscosity,γ is the shear rate, K is the consistency coefficient, and n 2 is the rheological index. When n 2 < 1, it indicates that the fluid shows shear thinning behavior, and belongs to pseudoplastic fluid in non-Newtonian fluid. Most industrial polymer fluids are pseudoplastic, and the colloidal fluid used in this damper is a pseudoplastic fluid as well. The MR fluid material is a viscoplastic fluid with a yield strength under the action of magnetic field. The relationship between yield strength τ y and magnetic flux density B is generally obtained through experiments. The MR damper uses LORD's product MRF-140CG magnetorheological fluid, and the relationship between magnetic flux density and yield stress was fitted by the Hill sigmoidal function, as shown in (3). Where, α is 117024.8, β is 0.96766, and η is 1.69363.
According to the theoretical viscoplastic equation shown in [32], [33], the dynamic viscosity of Bingham fluid is defined by the hyperbolic tangent function, as shown in (4).
where, α and ζ are numerical factors that control the growth of the non-Newtonian viscosity at very low shear rates, and µ p is the plastic viscosity. The hyperbolic tangent function is available to make the flow curve more stable at very low shear rates. Figure 4 shows the comparison between the two materials defined above and a typical rheological curve. Compared with the Standard Bingham model, the Tanh Bingham model improves the continuity of rheological properties at low shear rates.
It is of no significance that the compressibility of MR fluids also needs to be considered, as shown in (1). This is because that the pressure change at the moment of impact is extremely large, and if the fluid compressibility is not considered, the pressure calculation will be abnormal. The main parameters of the two fluid materials are shown in Table 2.

C. AIR SPRING
The air in the MR damper is at certain initial pressure, and its main function is to provide a restoring force to reset the compressed damper autonomously. The mechanical properties of air spring are generally described according to the gas polytropic process, as shown in (5). Table 3 shows some  parameters of the gas spring.
where, F s is the spring force rested on piston rod, P 0 is the gas initial pressure, A air is the cross-sectional area of air chamber, W 0 is the initial volume of air chamber, x is the displacement of piston rod, and n 3 is the gas polytropic exponent.

D. STATIC MAGNETIC FIELD FE MODELING
The finite element model of the magnetic field in MR damper was established in the multiphysics simulation software COMSOL/ Multiphysics, and the calculation of electromagnetic field was based on Maxwell's equations [6], as shown in (II-D), including the equation of continuity for constant electric charge density, Ampere's law, and Faraday's law.
where, H is the magnetic field intensity, J is the electric current density, A is the magnetic vector potential, B is the magnetic flux density, the relationship between magnetic field intensity and magnetic flux density satisfies the (7). µ r and µ o are the material permeability and the vacuum permeability, respectively. The multi-turn coil domain model in COMSOL is described by equations (8).
where, σ is the electrical conductivity, E is the electric field intensity, J e representing the external contribution of the electromagnetic coil to the ambient current density, N coil is the number of coil turns, I coil is the current in coil wire, and is the magnetic flux. The piston head and cylinder are made of 1045 carbon structural steel, and the piston rod is made of low permeability metal. The coil wire adopts the copper wire following AWG 24 standard. According to the symmetry, a two-dimensional axisymmetric mesh was generated, as shown in Fig. 5. There are 70375 quadrilateral cells in total in the computational domain.

E. FLUID CFD MODELING
In the CFD calculation of MR damper, it is necessary to consider turbulence and energy transfer, which involves fluid control equations, as shown in (9). The MR fluid in the damper follows three basic physical laws, i.e., the law of mass conservation, the law of momentum conservation and the law of energy conservation corresponding to three fluid governing equations respectively, namely continuity equation, Navier-Stokes equation and energy equation [40].
where, ρ is the fluid density, V is the fluid velocity vector, f b is the body force, p is the differential pressure, µ is the dynamic viscosity, h is the specific enthalpy, λ is the heat conductivity coefficient, φ is the viscous dissipation term, and S h is the source terms. Two-dimensional axisymmetric quadrilateral meshes were still used in fluid computing domain, and dynamic meshing layering technique was used for the movement of piston boundary. When the cells are compressed to a certain extent, the layer of cells close to the wall will collapse and merge. On the contrary, when the cells are stretched to a certain extent, the layer of cells close to the wall will split into a new layer. The time when the mesh changes depends on the splitting factor α s and the collapse factor α c . In order to realize the meshing layering technique, the Mesh Interface needs to be considered when generating the mesh. As shown in Fig. 6, the walls on both sides of piston were defined as rigid body moving walls, and only the fluid domain mesh was generated. The cells at the gap were densified, containing 95403 quadrilateral elements in total. VOLUME 8, 2020 When simulating the sinusoidal periodic motion and uniform quasi-static motion of damper, the boundary conditions of moving walls are explicit, that is, the velocity and displacement of moving walls at each time step are determined before calculation. However, when simulating the impact motion of damper, the boundary conditions of the moving walls are implicit, that is, the force of fluid on the piston determines the piston motion state, and the piston motion state also affects the flow of fluid. This is a passive motion problem of fluid-solid two-way coupling. Following the main principle to obtain the external forces and moments of the rigid body through numerical integration of pressure and shear stress over the moving surface, six-DOF module embedded in Fluent was used to solve the problem, and then to deal with the dynamic characteristics of the object [41]. In numerical calculations, the forces and moments on boundary are used for solution at each time step, and explicit Euler formula is used to describe the calculation of velocity and angular velocity, which is the basic theory of six-DOF module, as shown in (10).
where, t is the current time, t is the time step, v and ω are velocity and angular velocity of the center of rigid body respectively, F and M are the force and moment on the surface of rigid body respectively, and m and I are the mass and moment of inertia of the rigid body respectively.

F. DEFINITION OF USER DEFINED FUNCTION
The moving of the piston of MR damper will cause the influence region of magnetic field to change, therefore it is necessary to track and capture non-Newtonian region. Unlike explicit boundary model, the ''special viscosity'' region of implicit boundary model is not available for direct definition and dependent on real-time coordinates of moving walls. Therefore, user-defined functions (UDF) were written based on C language and compiled in the solver. The ''DEFINE_PROPERTY'' macro was used to define the dynamic viscosity of calculation domain. First, the ''Get_Domain'' macro and the ''Lookup_Thread'' macro were used to obtain the pointer of moving wall, on which the ''F_CENTROID'' macro was used to access the coordinate information of moving wall and store it in an array. Then,    In order to update the range in real time as piston moves, the circle center coordinate is a function of piston displacement. The red color of CFD is the non-Newtonian viscosity region.
information was accessed from the array, and the shear rate was obtained by the ''C_STRAIN_RATE_MAG'' macro. Finally, the ''C_CENTROID'' macro was used to access all mesh cells in the computational domain and assign the non-Newtonian viscosity function to the cells affected by magnetic field. The piston only moves in the specified direction, and there are air spring force and limit block, therefore it is necessary to define a special UDF file for six-DOF module. According to the header file six_dof.h, the ''DEFINE_SDOF_ PROPERTIES'' macro was used to define the properties of rigid body, and the code can continue to be written  into the above UDF file. The ''SDOF_MASS'' macro was used to define the mass, and its value takes the mass of projectile. The ''SDOFO_1DOF_T_P'' macro was used to define one translational degree of freedom. ''SDOFO_LOC, SDOFO_ MIN, SDOFO_MAX'' macros were used to define displacement limits. Finally, the ''SDOF_LOAD_F_X'' macro was used to write (5).

III. RESULTS AND DISCUSSIONS A. RESULTS OF MAGNETIC FIELD FE ANALYSIS
When the current I = 4.0 A was input into the coil, the magnetic flux density contour in the calculation domain is shown in Fig. 7. The figure shows that the maximum magnetic flux density is at the piston inner corner, and most of the magnetic induction lines pass through fluid into cylinder from both ends of piston. This is because that the magnetic permeability of the fluid and the piston rod is low, and the magnetic flux is more distributed in the magnetic material. The magnetic flux density distribution on the center line of the annular gap was output, as shown in Fig. 8. The curve shows that there are two peaks in the magnetic flux density at both ends of the annular gap, and the maximum value reaches 1.96 Tesla. In the middle and outside of the two peaks, the magnetic flux density decreases rapidly to nearly zero. The curve was output to provide UDF file for next CAD analysis.
It is of no significance that the main effect region of magnetic flux intensity in fluid domain is shown in Fig. 9. Being different from [32], the main region affected by magnetic field is not entirely in the annular gap, and there is overflow at both ends of gap. The calculation result of MR damper may be affected, therefore, it is necessary to write the ''arc'' region at both ends of annular gap in UDF code, which, combined with the magnetic flux density distribution curve shown in Fig. 8, can more accurately describe the non-Newtonian viscosity.

B. QUASI-STATIC CFD RESULTS AND EXPERIMENTAL VERIFICATION
Firstly, the quasi-static motion of the damper was simulated, which is a typical experimental process [1], [3], [14] and carried out on the INSTRON5982 universal testing machine, VOLUME 8, 2020 FIGURE 13. Velocity contours and streamtraces of the colloidal damper. The total time is 8ms, the maximum displacement is reached at 2ms. as shown in Fig. 10. Either one of the piston rod and cylinder was fixed, and the other was loaded and unloaded at a uniform velocity. Pressure sensors and displacement sensors were used to collect resistance and displacement data. Fig. 11 shows the mechanical characteristics of a uniform quasi-static motion under the same conditions. It is shown in the figure that the overall resistance of colloidal fluid damper is greater, and the resistance curve is steeper at the same velocity. The resistance curve of MR damper is relatively smoother, and the area of the hysteresis loop is larger. In addition, when comparing the hysteresis loops obtained from experiments and simulations, it shows that the model is accurate. Using quasi-static experiments to verify the model is the basis for next impact simulation.

C. DYNAMIC IMPACT CFD RESULTS
After verifying the accuracy of the two damper models, a dynamic impact simulation was performed. The dynamic characteristics of the damper loaded with projectile mass under a velocity impact of 7.4 m/s was solved. For this passive dynamic mesh problem, only an initial velocity needs to be given to the rigid body boundary, and the subsequent movement was automatically solved by the Six-DOF module with fluid-solid bidirectional coupling. Fig. 12 and Fig. 13 show the results of fluid velocity and streamline distribution during the impact. The total working time of the MR damper is 80 ms, and the maximum displacement is reached at 4 ms. However, the total working time of the colloidal damper is only 8 ms, and the maximum displacement is reached at 2 ms. During the movement of MR damper, many vortices are generated in the fluid, and the streamlines are extremely irregular. The larger vortices appear behind the direction of piston motion. Cavitation may occur at the vortex core. Though its influence on the mechanical properties of MR damper needs to be studied, it is not within the scope of this article. In contrast, the streamlines of the colloidal damper are relatively stable, no vortex is generated, and the region of maximum velocity remains in the annular gap.
The kinetic energy absorbed by the damper is mainly converted into the thermal energy of the fluid due to the viscosity of liquid, and the internal friction will generate heat. Fig. 14 and Fig. 15 show the fluid temperature distribution when the initial temperature is 300 K. The high temperature fluid is mainly produced in the annular gap where the shear rate is the largest. The temperature distribution of the MR damper is more uniform, and finally the average temperature on the right side is higher. The temperature of the colloidal damper is distributed near annular gap, and eventually most of the high temperature liquid flows to the left. This phenomenon is caused by the higher overall viscosity, weak fluidity, and short working time.
The kinematic data of the two dampers were compared, as shown in Fig. 16. At the same initial velocity, the maximum stretching distance of MR damper reaches 8.2mm, which is longer than that of colloidal damper of 6.5 mm. The compression process of the colloidal damper lasts 2.2 ms and the recovery process lasts 5.6 ms. The compression process of the MR damper lasts for 4 ms, while the recovery process lasts for about 80 ms. The working time of MR damper is much longer than that of colloidal damper, which greatly reduces  the rebound velocity. It is of no significance that, as shown in Fig. 16(a), the MR damper cannot be reset completely, and it almost stops moving at a displacement of 1.3 mm, this is caused by the properties of viscoplastic fluid. When rebounding, the fluid shear rate gradually decreases, and it can be seen from (4) that when the shear rate decreases to a certain extent, the shear stress of fluid will increase significantly. When the air spring force is not greater than VOLUME 8, 2020 the fluid resistance, the piston will slow down until it stops. However, this phenomenon will not affect the application of MR damper in the test shell. Disconnect the power supply to restore the liquid to Newtonian fluid, and the damper will be completely reset. The final velocity of the colloidal damper is 1.1 m/s, which shows that the energy absorption rate of the colloidal damper is 97.8%, and that of the MR damper is 100%. Fig. 17 shows the hysteresis loops obtained from impact analysis, and it can be seen that the impact curves of the two dampers are quite different, besides, the resistance curve of colloidal damper is more stable. The resistance peak value of MR damper is higher, reaching 110 kN.

IV. CONCLUSION
A one-way coupling numerical simulation method combining magnetic field FE and CFD was proposed to evaluate the mechanical properties of the two special dampers installed in the naval gun test shell. Aiming at the passive problem of dynamic impact, which is an implicit boundary, a CFD model involving multi-physics was established. After calculation and analysis, some conclusions were obtained.
(1) The dynamic characteristics of MR damper and colloidal damper were analyzed based on the combination of FEA and CFD numerical method. The results of quasi-static experiment and simulation have a high matching degree, which proves that the model is accurate. The energy absorption rate of colloidal dampers under impact is 97.8%, while the MR damper has a higher energy absorption rate, which is almost 100%. However, the structure of MR damper is more complicated, and with higher manufacturing and use cost.
(2) The analysis method of MR damper under dynamic impact was proposed and applied for the first time. In CFD analysis, the key points are to track and capture the non-Newtonian viscosity region, and to deal with the passive dynamic mesh. Based on the ''double peak'' magnetic flux density distribution obtained by FE analysis, the non-Newtonian viscosity of viscoplastic fluid was defined by hyperbolic tangent function. The ''circle equation'' was added to capture the non-Newtonian viscosity region accurately. At the beginning of the impact, the liquid compressibility was considered, and the force of air spring was described by gas polytropic equation. The mathematical models considered above were all written into user-defined functions, which are based on C language.
(3) During the dynamic impact process, multiple vortices with different sizes appeared in MR damper, while the streamline of colloidal damper was stable, and there was almost no vortex generation. The heat of damper was mainly generated in the annular gap. Due to the improved fluidity, the final temperature distribution of MR damper was more uniform, and the temperature of colloidal damper was concentrated near the annular gap. Compared with colloidal fluid damper, MR damper has higher resistance peaks, longer recovery time, and lower recovery velocity. In the case of continuous power supply, the MR damper will not be completely reset, but the circuit control system is available to solve this problem.