Adaptive magnetorheological fluid energy absorption systems: a review

In the last two decades, magnetorheological (MR) fluids have attracted extensive attention since they can rapidly and continuously control their rheological characteristics by adjusting an external magnetic field. Because of this feature, MR fluids have been applied to various engineering systems. This paper specifically investigates the application of MR fluids in shock mitigation control systems from the aspects of three key technical components: the basic structural design of MR fluid-based energy absorbers (MREAs), the analytical and dynamical model of MREAs, and the control method of adaptive MR shock mitigation control systems. The current status of MR technology in shock mitigation control is presented and analyzed. Firstly, the fundamental mechanical analysis of MREAs is carried out, followed by the introduction of typical MREA configurations. Based on mechanical analysis of MREAs, the structural optimization of MREAs used in shock mitigation control is discussed. The optimization methods are given from perspectives of the design of piston structures, the layout of electromagnetic coil, and the MR fluid gap. Secondly, the methods of damper modeling for MREAs are presented with and without consideration of the inertia effect. Then both the modeling methods and their characteristics are introduced for representative parametric dynamic models, semi-empirical dynamic models, and non-parametric dynamic models. Finally, the control objectives and requirements of the shock mitigation control systems are analyzed, and the current competitive methods for the ideal ‘soft-landing’ control objectives are reviewed. The typical control methods of MR shock mitigation control systems are discussed, and based on this the evaluation indicators of the control performance are summarized.

There are four operational modes of actuators based on MR fluids (Wang and Liao 2011, Gołdasz and Sapiński 2017), i.e. flow mode, shear mode, squeeze mode, and pinch mode.As shown in figure 4(a), the operational mode in which the MR fluid flows between two surfaces of fixed polar plates is called the flow mode or the valve mode.The shear mode is shown in figure 4(b) where the two polar plates are not fixed but constantly moving during the operation, and the moving direction of the polar plates is perpendicular to the direction of the applied magnetic field.In the squeeze mode, which is shown in figure 4(c), one polar plate in contact with the MR fluid is fixed, while the other moves along the direction of the magnetic field.Figure 4(d) displays the pinch mode, in which the two polar plates are fixed, and the flow of MR fluids results from the pressure difference in the flow direction (Elsaady et al 2020b).In the pinch mode, the ratio between the pressure difference and the flow velocity of MR fluids increases when the magnetic field is strengthened, which is different from in the valve mode where the ratio between the pressure difference and the flow velocity of MR fluids remains constant (Imaduddin et al 2013).Actuators developed based on MR fluids can work in one, two, or three operational modes (Phu and Choi 2019), and most MR fluidbased energy absorbers (MREAs) work in the flow mode.MREAs in the squeeze mode can be used for the applications requiring a small range of piston strokes, but the MREA force in the squeeze mode is larger than that in the shear mode or the valve mode, because of the larger compressive yield stress (Gong et al 2014).Usually, the force of MREAs can be increased by adopting mixed operational modes (Yazid et al 2014).

MREAs
MREAs have advantages of continuously adjustable force, wide dynamic range, fast response, simple structure, and low power consumption, thus they have been broadly used (Ahamed et al 2018) in the fields of automobile suspensions (Batterbee and Sims 2007, Bai et al 2013, Sun et al 2019, El Majdoub et al 2020), train suspensions (Liao and Wang 2003), vibration control of bridges (Gordaninejad et al 2002, Zhao et al 2019, Xu et al 2020, 2022) and seismic protection of buildings (Xu and Li 2008, Bitaraf et al 2010, Li et al 2013, Xu et al 2014, Caterino et al 2022, Huang and Li 2022).These applications have been extensively studied, and most of MREAs for vibration control work at low speeds (<1.0 m s −1 ).Furthermore, the theoretical system and technical methods have been developed for MREAs in the structural design (Zhu et     According to the form of excitations, Shou (2020) defined shock and vibration conditions.The shock signal can be regarded as a step signal approximately, while the vibration signal can be regarded as sinusoidal signals or random signals.
For both low-speed vibration control and high-speed shock mitigation, shock mitigation control systems based on MREAs needs to take states of the controlled object, such as the  velocity, the displacement, and the acceleration of payloads as the input of the controller.The system controller then calculates the desired MREA force for the current state of the controlled object.Some MR shock mitigation control systems are also equipped with force sensors to instantaneously detect the force output by MREAs, and the recorded data is used as the input of the system controller.Figure 5 shows the dynamic adjustment process of the MREA force in MR shock mitigation control systems.Using the inverse dynamic model based on the dynamic model of MREAs, the MREA controller solves the command voltage that needs to be provided to the current driver.Finally, the current driver provides the required current for the MREA to realize the optimal control of the MREA force.
The force provided by MREAs is limited to a certain range.By increasing the applied current from the minimum value allowed by MREAs to the maximum value, the MREA force increases from the lower limit to the upper limit.The size of this limited range affects the performance of MR shock mitigation control systems, and the force ranges provided by MREAs are mainly determined by the structural design.To enhance the control performance on vibrations/shocks, the MREA force should be accurately adjusted according to appropriate rules, depending on different input excitations.Hence, for structural design and parameter optimization of MREAs, it is important to develop a dynamic model that can accurately describe the characteristics of MREAs and design an effective control method for the MR vibration/shock control system.

MREA design and characteristics description
As the first step for the structural design of MREAs, it is necessary to establish an analytical model based on structural parameters and the piston velocity by analyzing the relationship among the piston velocity, the applied current, and the MREA force.The most common optimization goal is to maximize the dynamic range of the MREAs.Constraints mostly include geometric sizes, the MREA force at zero or no current (i.e.passive MREA force or uncontrollable MREA force), and the requirement of magnetic saturation.Then the objective function can be defined according to the performance demands.When the initial structural parameters are given, the unknown parameters can be optimized by the corresponding method, and the optimized parameters can be further adjusted with fluid dynamic analysis and magnetic circuit simulation.Finally, all the parameters in the analytical model can be determined, and the structural design of MREAs can be completed.
Based on the in-depth research into the flow characteristics (Xu and Sun 2021, Martínez-Cano et al 2022, Sharmili et al 2023) and the manufacture of MR fluids (de Falco Manuel et al 2023, Kubík et al 2023), the optimal design of MREAs operating for low-speed vibrations has been extensively studied (Nguyen et al 2007, Bagherkhani and Mohebbi 2022, Liu et al 2022), but literature on the optimal design of MREAs under high-speed impact conditions is not much.
Based on the Bingham-plastic (BP) nonlinear flow model with minor losses proposed by Mao et al (Mao et al 2008(Mao et al , 2009)), Singh et al (2014) designed a double-ended MREA to describe the characteristics of the MREA more accurately in case of impacts.With the constraint of the maximum allowable passive MREA force, the optimization problem is solved by using the weighting method and the dual objective function constraint method.One of the objective functions is designed to maximize the controllable MREA force, and this objective function is used as the cost function.In the second objective function, the maximum allowable passive MREA force is constrained by a parameter.In a design of the MREA for aircraft landing gears, Batterbee et al (2007) chose the conventional oil-gas EA structure to preliminarily estimate a set of design parameters, including the length and the outer diameter of the MREA, and then the optimal design was proposed considering the quasi-static model and the magnetic circuit design theory.For the bi-fold MREA used in helicopter landing gears, Saleh et al (2019) collected constraints of the geometric size, the minimum allowable passive MREA force, and the requirement of magnetic saturation.The genetic algorithm was adopted to find the interval where the optimal solution of each design variable is located, and then the nonlinear sequential quadratic programming was used to find the optimal solution.
The model of MREAs can be divided into static MREA models and dynamic MREA models.Static MREA models can be used for the structural design of the MREA, and dynamic models can be employed to design the MREA controller.Based on dynamic models, the mathematical relationship between the applied current and the MREA force can be established for the MREA in the system controller.The system controller calculates the desired MREA force according to the working condition, and the MREA controller calculates the current input, which is required to provide the MREA force equal to the desired MREA force calculated by the system controller.Therefore, the accuracy of the dynamic model to describe the characteristics of the MREA will affect the performance of the control system due to the deviation between the MREA force and the desired force.
As a common way to build up a dynamic MREA model for describing the dynamic characteristics of MREAs, elastic elements, damping elements, friction elements, and other mathematical operators with special meanings are utilized in a specific way to establish the mapping between the input and the output.Models built in this way are also called parametric models, and classical parametric models include the Bingham model (Stanway et al 1987), the nonlinear bi-viscous model (Stanway et al 1996), and the Bouc-Wen model (Ismail et al 2009), etc.
In order to further improve the accuracy of the dynamic MREA models, especially the accuracy when describing the dynamic characteristics of MREAs at high speeds, researchers developed classic parametric models.Spencer et al (1997) proposed a modified Bouc-Wen model, in which the three constant items in the Bouc-Wen model were redefined with cubic polynomials of the applied current, and the accuracy of the modified Bouc-Wen model was enhanced.Bai et al (2015) proposed a restructured model using the Bouc-Wen hysteresis operator, which is advantaged when the structure of modified Bouc-Wen models is too complex.In addition to the developments of classic parametric models, aiming at improving the prediction accuracy of dynamic characteristics of MREAs under high-speed impact conditions, dynamic models based on resistor-capacitor (RC) operators (Chen et al 2018, Bai et al 2019a) were also proposed in recent years.
Based on both testing data analysis and device working principles, the mapping function between the output and the input is established with physically meaningless parameters to describe the dynamic characteristics of MREAs.The models based on mapping functions with parameters of no physical meaning are also called non-parametric models (Shou 2020).Choi et al (2001) employed polynomials to establish a mapping function to describe the input and output characteristics of MREAs.With the developments of machine learning technology, neural network models were also used as nonparametric models to predict the strong nonlinear characteristics of MREAs under high-speed impact conditions (Boada et al 2011, Shou 2020, Bahiuddin et al 2021).The neural network models include feed-forward neural networks trained by extreme value learning machines (Bahiuddin et al 2021), neural networks models based on recursive lazy learning methods (Boada et al 2011), and neural networks models using the adaptive differential evolution algorithm to optimize parameters for the support vector regression (Shou 2020).In addition, Bai and Tang (2021) combined a dynamic model based on the RC operator with a neural network model and proposed a dynamic model based on the dynamic RC operator.Parametric models provide clear physical meanings of model parameters, and these non-parametric models show superiority in prediction accuracy.

Control methods
As relatively mature methods, conventional control methods including the Skyhook control (Yao et al 2002), the Sky-Ground hook control (Ahmadian and Pare 2000) In contrast to low-speed vibration problems, the high-speed impact problem has characteristics such as large impact loads, high impact speeds, short event time, and uncertainties (Shou 2020).A suitable control method is necessary for good performance of shock mitigation controls.At present, a variety of control methods have been applied to different MR shock mitigation control systems.
Current control methods for aircraft landing gear systems are mainly aimed at improving energy absorption efficiency, emphasizing the absorption of more impact energy in the same piston stroke.Han et al (2019) applied Skyhook controllers to the aircraft landing gear.Due to the insufficient adaptability of Skyhook controllers with fixed gains in different landing conditions, a hybrid controller was proposed to enhance the impact mitigation performance of the landing gear.Considering the energy absorption efficiency of 100% as the control objective, Yoon et al (2020) proposed a control method to achieve higher efficiency of absorbing impact energy than the Skyhook controller.The adaptive sliding mode hybrid control was applied to the aircraft landing gear by Luong et al (2020), and the energy absorption efficiency was also enhanced by comparing it with the Skyhook controller.
Hu et al (2012) applied the PID control, the on-off control, and the single-input fuzzy control in an MR gun recoil system, respectively.Li et al (2019) proposed an optimal control fuzzy-compensated for MREA in a gun recoil system.
Compared with the optimal control method or the twodimensional fuzzy controller, this hybrid control method has higher energy absorption efficiency and better shock mitigation stability.
Maeda et al (2019) used an MREA landing gear to replace the conventional plastic deformation landing gear on a lunar lander with a mass of 180 kg and proposed a control method based on the Lyapunov function to reduce the roll rate of the lunar lander on the uneven ground.Both the simulation and the prototype experiment showed that the lander equipped with MREA landing gear can effectively prevent the lander from overturning while landing and improve the adaptability to the landing terrain and the landing capability for complex landforms.Wang et al (2019) proposed a new structure using the bypass MREA as the primary strut of landers and analyzed the feasibility of applying this landing gear structure to a heavyduty lander with a mass of 1350 kg.A fuzzy controller was developed for MREAs with the main control objective of preventing the lander from overturning.The acceleration of the lander cabin, the derivative of acceleration to time, the pitch angle, and the roll angle of the cabin were selected as the input of the controller, and the voltage values of the four MREAs were taken as the output of the fuzzy controller.The simulation proved the effectiveness of the MREA applied to heavy-duty landers, and the capability of landing on complex and unexplored terrains was significantly improved compared with the lander based on honeycomb aluminum structures.
In shock mitigation control systems, the acceleration transmitted to the protective structure during the shock mitigation is required to be less than the allowable acceleration because the protective structure and the occupants are greatly affected by the transmitted impact.At the same time, the peak value of acceleration should be as small as possible during the shock mitigation, and the acceleration should be decreased monotonically after reaching the peak value without any fluctuations.In order to meet the above requirements, the available stroke of the MREA should be utilized to the greatest allowable extent to reduce the transmitted impact loads as much as possible.
In the adaptive MR shock mitigation control system, 'softlanding' is often taken as the control objective to achieve optimal control performance.Soft-landing means that when different impact excitations are exerted on the shock mitigation control system, the piston velocity drops to zero exactly when the piston stroke reaches the maximum (Wereley et al 2011).At this time, the impact energy of the payload is completely absorbed.Compared with the case where the velocity of the payload is reduced to zero before reaching the maximum stroke, the acceleration curve of the payload at soft-landing is smoother and the peak value of acceleration is smaller when the control objective of soft-landing is achieved.
The MR impact mitigation control system and the related technologies are shown in figure 6 to realize the soft-landing.Considering both the maximum allowable acceleration of the payload and the boundary condition of the soft-landing, according to which the piston velocity is zero when the piston stroke is the maximum, the system controller calculates the desired controllable MREA force required based on the states of the protective structure.The inverse dynamic model can be employed by the MREA controller to calculate the applied current that needs to be provided to the electromagnetic coil of MREAs, and the current driver finally provides the required current for the MREA to achieve the control objective of soft-landing.
In recent years, there has been research on soft-landing control in single-degree-of-freedom (SDOF) systems.Wereley et al (2011) proposed the optimal Bingham (Bi) number control method to realize the soft-landing control of the drop hammer (i.e., the payload under drop-induced impact).The control method of optimal Bi number was subsequently applied to a SDOF system considering the MREA force control in the rebound stroke (Singh and Wereley 2013), and the time lag phenomenon of MREA response was also taken into account in the follow-up research (Choi and Wereley 2015).In order to meet the limitation of the deceleration threshold of the payload on specific occasions, the minimum duration deceleration exposure (MDDE) control method was proposed by combining the maximum constant deceleration (MCD) control and the optimal Bi number control (Wang et al 2021a), and the quadratic damping phenomenon of the passive MREA force at high impact speeds was also considered (Wang et al 2021b).In addition, Bai and Yang (2021) used the constant force control method to realize the soft-landing control test of the drop hammer.Li et al (2023a) studied the influence of dynamic models on the control performance of the SDOF MR shock mitigation control system.The system controller was established by using the dynamic model based on the RC operator and the Bingham model, respectively.These two system controllers combined with the constant force control method or the optimal Bi number control in pairs, and the four combined SDOF MR shock mitigation control systems were tested.The velocity-acceleration conversion ratio (V-ACR) and the average velocity change rate (AVCR) were also proposed as two evaluation indicators.Li et al (2023b) also proposed a two-degree-of-freedom (2DOF) optimal Bi number control method, which was applied to a shock mitigation control system of the 1/4 automobile suspension.As the research on the soft-landing control of the SDOF system is being continuously improved, the soft-landing control of systems, such as MREA-based gun recoil systems, aircraft landing gears, and planetary landers, will become one of the main directions in the research field of the MR impact mitigation control in the future.
According to figure 6, in order to realize the reasonable performance of MR shock mitigation control systems, it is necessary to optimize the structure of MREAs according to the requirements of the control system (Zhu et al 2012).The structural design of MREAs should be guaranteed to provide the MREA force that the controller need, and the structural design of MREAs can be guided by precise analytical models.When the structure of MREAs is designed and confirmed, a dynamic MREA model is required to correctly predict the strong nonlinear characteristics of MREAs during the impact event and to accurately calculate the required current input based on the desired MREA force of the controller.Finally, an effective control method is needed to achieve the control objective of soft-landing under different impact conditions.
As shown in figure 6, this paper will introduce the classic structural configuration of MREAs and the structural design criteria of MREAs for high-speed shock mitigation control in the second section.Subsequently, the static and dynamic MREA models of MREAs, as well as the dynamic models of MREAs suitable for high-speed impact control will be introduced in the third section.Finally, the current control methods for the MR shock mitigation control system will be introduced in the fourth section, which will focus on the control method with soft-landing as the control objective.The related evaluation indicators of the shock mitigation control performance will be also reviewed.Figure 7 is a typical structural schematic diagram of MREAs.For a conventional EA, the force is generated when the piston runs and compresses the lower chamber.Although the fluid can flow from the lower chamber to the upper chamber through the gap, the volume of fluid flowing into the upper chamber is much smaller than the volume swept by the piston.The fluid in the lower chamber is compressed and the pressure rises, while the pressure of the fluid in the upper chamber reduces.The viscous pressure drop caused by the movement of the piston is generated.In addition, the gas compensation chamber is separated at the bottom of the MREA by a floating piston in the lower chamber of the MREA.The chamber is utilized to compensate for the volume change in the cylinder caused by the movement of the piston.When no current input is applied, the MREA can be regarded as a passive EA, and its force, F passive , can be expressed as:

Designs of MREAs
where F gas is the spring force generated by the accumulator; ∆P η is the viscous pressure drop due to the movement of the piston; A P is the effective piston area, which can be written as: where r p is the radius of the piston and r pr is the radius of the piston rod.
The viscous pressure drop, ∆P η , in equation ( 1) can be written as: where ρ is the density of MR fluids; f is the Darcy friction factor; L = (L a + L c ) is the total length of the MR fluid gap; L a is the total active length of the MR fluid gap where the MR fluid will be activated to produce the MR effect; L c is the length of all electromagnetic coils in MREAs, which is also the total length of the non-activated area in the MR fluid gap; h is the annular fluid gap, and v f is the fluid velocity in the annular duct or valve, which can be expressed as: where ẋp is the piston velocity and A h is the cross-sectional area of the MR fluid gap.
According to equation (1), the cross-sectional area, A h , of the fluid gap can be adjusted for the passive EA during the working process, which will increase the structural complexity.When a magnetic field is applied perpendicular to the direction of flow, the MR fluid will transform from a Newtonian fluid to a non-Newtonian fluid (Gedik et al 2012), and the yield stress increases as the magnetic field increases.The pressure drop due to MR yield stress, ∆P MR , can be controlled by the applied current.
Combining equations (1)-(3), the total force of MREAs, F total , can be defined as: where ∆P MR is the pressure drop due to the MR yield stress, which can be written as: where τ y is the yield stress of the MR fluid, which can be taken as a function of the magnetic field strength (and the applied current).

Configurations of MREAs
The single-ended MREA shown in figure 7 has the advantages of lower zero-field frictional force, less required piston rod space, and ease of installation and arrangement.But the damping capabilities of the single-ended MREA are, to some extent, reduced by the spring effect of the accumulator.A doubleended MREA has been developed to achieve better mechanical performance as shown in figure 8.When the MREA is in the compression stroke, the fluid amount in the lower chamber pumped by the piston of single-ended dampers is larger than the empty volume of the upper chamber because of the volume occupied by the piston rod in the upper chamber.In order to compensate for this volume difference, the accumulator is configured at the lower chamber.For double-ended MREAs, since the piston rod runs through the entire cylinder of MREAs, there is no volume difference between the lower chamber and the upper chamber.Therefore, there is no need to set up an accumulator.As seen from figure 8, the removal of the floating piston and the gas compensation chamber allows the MREA to increase the effective stroke of the piston at the same axial size while the structure of MREAs was also simplified.In addition, double-ended MREAs do not have a spring effect caused by the accumulator (Yang et al 2002) and have better mechanical performances.Because of the accumulator, single-ended MREAs have both MREA force and spring force, and the spring force can be harmful to some shock mitigation applications because of rebounding.However, since the piston rod runs through the entire cylinder of the double-ended MREA, it is necessary to reserve a space for the movement of the piston rod.Therefore, compared with single-ended MREAs, double-ended MREAs require a larger axial dimension.
Figures 9(a)-(c) show the monotube MREA, the twin-tube MREA, and the bypass MREA, respectively.In figure 9(a), it is a monotube MREA with two electromagnetic coils, and this MREA consists of an MR fluid reservoir, an accumulator, a piston and a piston rod, etc.The monotube MREA has a simple structure (Ebrahimi 2009, Zhu et al 2012).Compared with the monotube MREA, the twin-tube MREA has an inner cylinder and an outer cylinder, as shown in figure 9(b), and the structure of the outer cylinder can protect the internal parts of this twintube MREA.The inner cylinder is filled with an MR fluid, the piston/piston rod assembly moves in the inner cylinder, and the outer cylinder is also partially filled with an MR fluid.The MR fluid can flow between the inner and outer cylinder through a valve assembly at the bottom of the inner cylinder, which can be used to balance the volume change in the inner cylinder caused by the movement of the piston.The increase/decrease of the volume occupied by the piston rod in the inner cylinder is equal to the volume of the MR fluid flowing into/out of the outer cylinder through the bottom valve assembly.The flow of the MR fluid through the valve assembly also generates an additional force, which increases the passive MREA force.
Since the energized electromagnetic coil generates additional heat energy, the temperature of the MR fluid will rise, and an excessively high temperature will deteriorate the property of the MR fluid (Dong et al 2017, Kariganaur et al 2022).In order to control the temperature of the MR fluid, the performance of heat dissipation should be considered when designing the MREA.Due to the outer cylinder, the heat dissipation of the twin-tube MREA is not as good as that of the monotube MREA (Ebrahimi 2009).When working vertically (if MREA is installed vertically), on the other hand, the long-term downtime of MREAs may result in the sediment of MR fluids, which will easily cause the bottom channel to be blocked, then the MREA is unable to work normally.
In order to increase the controllable MREA force, an MREA with double electromagnetic coils was developed based on the structure of MREAs with only one electromagnetic coil, as shown in figures 9(a) and (b).According to equation ( 6), the length of the magnetic pole in the MR fluid gap can be increased by adding a second set of electromagnetic coils, and the total active length of the MR fluid gap where the MR effect occurs can be subsequently increased.At the same time, the strength of the external magnetic field acting on the MR fluid is also increased, thereby improving the yield stress, τ y .Increasing the number of electromagnetic coils and the MR yield stress, τ y , can result in an increase of the pressure drop due to MR yield stress, ∆P MR .As the controllable MREA force increases with the increase of ∆P MR , the dynamic range of MREA also increases.The dynamic range is defined as the ratio of the total MREA force at field-on status to the total MREA force at field-off (i.e.no field) status.So, the larger the dynamic range, the greater the controllable range of the MREA force (Bai et al 2019b).Some MREAs for shock mitigation control systems use three or more sets of electromagnetic coils.
In order to further improve the performance of heat dissipation, a bypass MREA with the electromagnetic coil placed outside the cylinder was developed based on the monotube MREA, as shown in figure 9(c).The advantage of this configuration is that the heat generated by the electromagnetic coil can be directly dissipated into the atmosphere, avoiding the decrease in the apparent viscosity of MR fluids (Hitchcock et al 2007).While the disadvantage is that with the same volume and electromagnetic coils the bypass MREA generates a smaller MREA force than the MREA with internal coils (Grunwald andOlabi 2008, Zhu et al 2012).
According to equations (3) and (4), by increasing the crosssectional area of the MR fluid gap, A h , the fluid velocity in the annular duct or valve, v f , can also be reduced, thereby achieving the purpose of reducing the viscous pressure drop, ∆P η .When the radial size is allowable according to the requirement for MREAs in MR shock mitigation control systems, the structure shown in figure 9(c) can be adopted, and the MR fluid gap and the electromagnetic coil can be arranged outside of the cylinder.The cross-sectional area, A h , can be increased with this structure to reduce the passive MREA force and to improve the dynamic range without reducing the controllable MREA force (Carlson 2009).
Most of MREAs used in typical MR vibration or shock mitigation control systems adopt the monotube structure with a valve-in-piston-head configuration (Zhu et al 2012), and it is common to incorporate two or more sets of electromagnetic coils.The structure of multiple electromagnetic coils can increase the volume of MR fluid activated by the magnetic field, thereby increasing the controllable MREA force within the limited size constraints (Goncalves 2005).Adopting monotube structures with a valve-in-piston-head configuration makes the design of MREAs more compact, and also makes the passive force of MREAs smaller at zero current, thus the controllable MREA force is relatively large when the applied current is applied.
In order to improve the landing performance of the aircraft, Luong et al (2020) developed a monotube MREA for the landing gear, as shown in figure 10.The MREA works in the flow mode, and two sets of electromagnetic coils are integrated with the piston rod.The diameter of the piston rod of this MREA is large, so a gas compensation chamber (accumulator) and a fluid chamber filled with an MR fluid are set in the hollow part of the piston rod.The sliding bearing installed on the piston moves with the piston, dividing the space of the cylinder into upper and lower chambers, and the lower chamber communicates with the chamber in the piston rod, then both of them serve together as the lower chamber of the MREA.It makes the axial dimension of this MREA more compact than that of the MREA shown in figure 7 with the same MREA force and piston stroke.As the piston moves, the MR fluid flows from one chamber to the other through two annular gaps arranged in the piston.There is a relief valve in the piston to make the As a result, the force of this MREA is asymmetric, and the MREA force of the extension stroke is smaller, which meets the requirements on characteristics of the MREA force for the landing gear system of aircrafts.
MREAs are also used in the artillery/gun recoil system to solve the problem that the force of the passive EA gradually decreases with the piston velocity.MREAs are capable of producing increasing controllable MREA force to compensate for the decrease of the passive MREA force, thereby minimizing the recoil stroke of the artillery/gun.Li et al (2019) developed an MREA for gun recoil systems as shown in figure 11, which works in the shear-flow mixed mode.The monotube configuration was adopted in this MREA, where three sets of electromagnetic coils were integrated with the valve-in-the-pistonhead configuration.Bai et al (2013) proposed the internal bypass structure MREA to improve the control performance of shock mitigation for the automobile suspension system, as shown in figure 12. Through theoretical modeling and numerical simulations, the principle and the mechanical behaviors of the internal bypass MREA were further studied (Bai et al 2019b), followed by prototype experiments under high-speed impact conditions (Bai et al 2018).The MREA adopted a twin-tube structure configuration with multiple sets of electromagnetic coils, and the bypass gap between the inner cylinder and the outer cylinder was used as the MR fluid gap.When the piston moves, the MR fluid flows between the upper and lower chambers through the bypass.Five sets of electromagnetic coils were evenly wound on the outer wall of the inner cylinder, so that the MR effect can be generated in the entire MR fluid gap, thereby increasing the MREA force.
The single-ended MREA has the advantages of a simple structure, compactness, and easy installation.The drawback is the necessity of a high-pressure accumulator.The doubleended MREA does not require an accumulator, resulting in better damping characteristics.But, for a longer piston stroke, more space is required.The monotube MREA has the advantage of a simple structure but is more susceptible to being damaged by external impacts.The twin-tube MREA has a more robust structure, and its accumulator requires lower gas pressure, resulting in better mechanical performance.However, its structure is more complex.The bypass MREA has the advantages of simplicity in structure design and assembly, as well as ease of maintenance.Designing a bypass MREA based on a passive EA is less challenging (Zhu et al 2012).In addition, the bypass structure can achieve a larger dynamic range, providing greater damping force to meet the needs of shock mitigation control.However, it is less compact due to the bypass duct, which is also more susceptible to being damaged.

Design strategy of MREAs for shock mitigation
To further analyze the structural design of MREAs suitable for shock mitigation control, figure 13 shows the MREA force and the piston velocity of three different MREAs.The structure of the unicoil MREA is the same as that in figure 8.The multicoil MREA has a twin-tube structure with multiple sets of electromagnetic coils, and the electromagnetic coils are increased based on the structure shown in figure 9(b).The internal bypass MREA is illustrated in figure 12.According to equation ( 5), when the spring force generated by the accumulator and the minor losses are negligible, the total MREA force is composed of the controllable MREA force, A P ∆P MR , and the passive MREA force, A P ∆P η .By adjusting the applied current, the yield stress of the MR fluid, τ y , can be controlled in the range of zero to the maximum value, so that the controllable MREA force can be tuned from zero to the maximum value.According to figure 13, the controllability of the internal bypass MREA is the best among the three MREAs, because it provides a larger range of MREA force at most piston velocities.Bai et al (2019b) introduced concepts, such as the dynamic range and the controllable velocity range, to evaluate the controllable characteristics of different MREAs.It can be deduced from figure 13 that when the passive MREA force is almost the same, the greater the controllable MREA force, the greater the dynamic range.However, at high piston velocity the controllable MREA force will be restricted by the allowed maximum MREA force due to the structure of MREAs, and the maximum controllable MREA force cannot actually be achieved.As a result, the dynamic range of MREAs decreases as the piston velocity increases.In order to further improve the dynamic range, it is necessary to reduce the passive MREA force, F η .As shown in figure 13, the passive MREA force increases with the piston velocity.When the passive MREA force reaches the maximum allowable MREA force at the piston velocity of ẋpmax in figure 13, the controllable MREA force, A P ∆P τ , has to be zero to prevent the MREA from being damaged.If the passive MREA force, F η , can be reduced, the MREA force range will be further increased at the same piston velocity, and the allowed maximum piston velocity, ẋpmax , will also increase, which is beneficial for the shock mitigation control.
The controllable velocity range is defined as the range of the piston velocity between ẋpmin and ẋpmax .Here, ẋpmin is the minimum piston velocity at which the total MREA force reaches the allowed maximum MREA force, F thr , when the applied current, I, is the maximum value.While ẋpmax is the maximum allowable piston velocity when the applied current, I, is zero (Bai et al 2018).The piston velocities, ẋpmin1 , ẋpmin2 , and ẋpmin3 are the lower piston velocity bounds of the three MREAs for a specific MREA force, respectively.In order to improve the performance of MREAs, the larger range of controllable velocity can be expected to be acquired with the smaller velocity, ẋpmin .The multicoil MREA has a smaller range of controllable velocity, which can only provide the desired damping force equal to F thr when the piston velocity is greater than or equal to ẋpmin2 .
On the basis of the analysis above, the diameter of the piston rod in the MREA shown in figure 9 is relatively large, which reduces the effective area, A P , in equation ( 2) and the fluid velocity, v f , in equation ( 4), thereby achieving the purpose of reducing the viscous pressure drop, ∆P η , in equation (3).Although this design increases the ratio of the controllable MREA force to the total MREA force, the value of the controllable MREA force decreases with the decrease of the effective area, A P .To enhance the controllable MREA force, two electromagnetic coils are used, then the active length of the MR fluid increases compared to that of a single electromagnetic coil.At the same time, the yield stress, τ y , also increases with the strengthening magnetic field-induced pressure drop, ∆P MR , in equation ( 6) increases.This has a greater effect on the controllable MREA force than the reduction of the effective area, A P .Therefore, the controllable MREA force can be effectively increased by this structure.The total length of the MR fluid gap increases with the increase of the total active length of the MR fluid gap, making the numerator in the expression of viscous pressure drop in equation (3) increase, resulting in an increase in the viscous pressure drop and an increase in the passive MREA force.However, the decrease of the average flow velocity in the MR fluid gap has a greater effect on the viscous pressure drop than the increase of the activation length of the MR fluid does.This is because the viscous pressure drop, ∆P η , has a quadratic relationship with the flow velocity, v f , but a linear relationship with the activation length of the MR fluid.Thus, the passive MREA force shows a downward trend.
The method of structural design adopted by the MREA in figure 11 is similar to the MREA in figure 10, and the number of electromagnetic coils is increased to three.Both the structures shown in figures 9 and 10 adopt the valve-in-piston-head configuration, which will increase the radius of the piston.The radius of the piston rod, r pr , needs to be increased to make the effective area of the piston, A P , small enough, and then the mass of the MREA increases.In addition, the length of the piston will be increased in the valve-in-piston-head configuration with multiple electromagnetic coils, resulting in a reduction in the maximum piston stroke.Larger axial size and mass are required to guarantee the constant maximum of the piston stroke (Bai et al 2019b).
As shown in figure 12, the internal bypass MREA can solve the problem of the structures as shown in figures 9 and 10 to some extent.The internal bypass MREA provides a larger MREA force and a larger piston stroke with smaller mass and axial size, while also ensuring a larger controllable velocity range and larger dynamic range for the MREA.
When the axial dimension of the piston is the same as the structure of the MREA represented in figure 8, the internal bypass MREA arranges the electromagnetic coil on the outer wall of the inner cylinder to facilitate the arrangement of multiple electromagnetic coils, and then the activation length of the MR fluid is increased.Compared with the MREA structure modified by increasing electromagnetic coils on the basis of figure 9(b), the axial dimension of the piston is reduced significantly, so that the piston stroke can be greatly increased with the same activation length of the MR fluid.In addition, the effective piston area, A P , of this MREA is smaller than that of the structure in figure 9(b), and the annular valve gap, h, of this MREA is also comparatively small.According to equations (3) and (4), the above combined effect makes the passive force of the MREA in figure 12 equal to that of the MREA derived from the structure in figure 9(b).Since the decrease of the valve gap, h, in the structure shown in figure 12 has a great effect on the controllable MREA force than the decrease of the effective area, A P , does.Thus, the maximum controllable MREA force of the MREA in figure 12 is larger than that of the MREA structure in figure 9(b).As shown in figure 13, before the piston velocity reaches the velocity, ẋpmin2 , the dynamic range is larger than that of the other two MREAs, meanwhile, the controllable velocity range is also larger.Therefore, the internal bypass MREA shown in figure 12 is more suitable for shock mitigation and vibration control systems, and the mass of this MREA will be lighter.
To improve the control performance of MR shock mitigation control systems, the response time of the system should be reduced in the structural design of MREAs.Koo et al (2006) proposed the definition and the experimental measurement method of response time and analyzed the main factors influencing the response time of MREAs.According to the conventional definition of system response time, only the response time of the actuator is considered.For the MR shock mitigation control system, Bai et al (2018) pointed out that the response time of the feedback sensor, the system controller, and the MREA controller in MR shock mitigation control systems should not be ignored when mitigating a short duration of the impact event.The response time is composed of the activation time and the rise time.The activation time is the starting time of the output signal.The rise time is defined as the time when the output signal reaches a certain proportion of the steady state, starting from the end of the activation time.The response times of the system, t total , and each unit, t rj , can be respectively expressed as: where the response time of the system, t total , is defined as the sum of the response time of the jth unit in the control system.t rj is the response time of each component such as the feedback sensor, t r0 , the data acquisition and control algorithm calculation unit, t r1 , the current driver, t r2 , and the MREA, t r3 .Also, t aj and t risej are the activation time and the rise time of the jth unit, respectively.The activation time of the jth unit depends on its own characteristics which reflects the real-time controllability of this unit, while the rise time of jth unit is affected by the rise time of the j-1th unit.The activation time of the jth unit starts during the rise time of the j-1th unit.Thus, the overlapped time between the activation time/rise time of the jth and the activation time/rise time of the j-1th unit should be subtracted when calculating the response time of the jth unit, which also makes the definition of the response time of the jth unit only reflect the real-time controllability of the jth unit.In the MR shock mitigation control systems, the feedback sensor is the beginning of the signal flow chain, and the activation time of the feedback sensor, t a0 , is zero.Gu et al (2016) conducted a systematic study on how to reduce the response time of various MREAs from the aspect of electromagnetic coil layout, and the research pointed out that the inductance and driving electronics of MREAs are the main factors producing response delay time.The response time of MREAs can be effectively reduced by some methods of structural design like dividing a large electromagnetic coil into multiple parallel small electromagnetic coils and introducing an electromagnetic coil for magnetic field quenching to offset the magnetic hysteresis effect in the yoke and MR materials when the magnetic field decreases.Strecker et al (2015) also studied structural design methods to reduce the response time of MREAs.Ding et al (2023) considered the response time for the MREA-based automobile suspension system and proposed the time-delay-dependent H∞/H2 optimal control method.
Based on the above analysis, the following methods of structural design can be adopted to expand the force range provided by MREAs: (i) increasing the number of electromagnetic coils to increase the active length of the MR fluid; (ii) using inner bypass structure to reduce the radius of the piston by decoupling the electromagnetic coil from the piston; (iii) increasing the radius of the piston rod; (iv) reducing the annular valve gap (i.e. the MR fluid gap) based on the comprehensive application of the first three methods; (v) considering the response time in MREA design, modeling as well as control.

Brief introduction
According to the modeling method, MREA models can be divided into analytical models and dynamic models.The analytical models are based on the constitutive model of MR fluids, adopting the Navier-Strokes equation and the continuity equation, according to corresponding boundary conditions and initial conditions, to calculate the output MREA force.For example, an analytical model of MREAs in flow mode can be expressed by equations ( 1)-( 6).Since the analytical models include the rheological characteristics of MR fluid and the geometric parameters of the MR valve, they are mainly used for the preliminary structural design and the parameter optimization of MREAs.Analytical MREA models are further divided into quasi-static models that ignore the inertial effect of MR fluid and unsteady models that consider the inertial effect of MR fluid.When the MREA is subjected to high-speed impacts, the influence of the inertia effect cannot be ignored, so the unsteady model is often used to provide guidelines on the structural design of MREAs for shock mitigation applications.This review work will focus on the unsteady models of MREAs.
The dynamic MREA models refer to the models that are established based on experimental data and accurately predict the dynamic characteristics of MREAs.It is necessary to acquire the optimal parameters of mechanical MREA models through an algorithm to minimize the error between the model simulation results and the experimental data.This process is called parameter identification.Since the dynamic MREA models are based on a series of experimental data of MREAs, they are also called phenomenological models.Dynamic MREA models are divided into parametric models that equate the MREA with a combination of linear/nonlinear springs, dampers, and hysteresis operators, and non-parametric models with parameters that have no physical meaning.
Parametric models are composed of elastic elements, damping elements, friction elements, and other hysteresis operators with special meanings.Commonly used parametric models include Bingham-model-based dynamic models, Bouc-Wen hysteresis operator-based models, Dahl hysteresis operator-based models, bi-viscous models, hyperbolic tangent function-based models, and RC operator-based dynamic models.These models are mainly used to describe and predict the MREA force-velocity hysteresis characteristics of the MREA under sinusoidal excitations (i.e.low-speed vibration conditions), but there are few studies on the parameter models applied to MREAs under high-speed impact conditions.It is a great challenge to establish an accurate analytical model or a parametric dynamic model to predict the highly nonlinear behavior of the MREA under high-speed impact conditions.In this case, the turbulence and the fluid-solid coupling effect should be considered in dynamic MREA models, which will increase the difficulty of establishing the model.
In the mapping relationship established by non-parametric models, the parameters have no actual physical meaning.The non-parametric models, such as polynomial models, neural network models, and fuzzy logic models, have the advantage of strong robustness, and they are suitable for linear, nonlinear, and hysteresis systems.In order to describe the characteristics of MREAs under high-speed impact conditions, Shou (2020) proposed a hybrid dynamic model of MREAs by using support vector regression and employed an adaptive differential evolution algorithm to optimize the regression parameters.
The dynamic models of MREAs are the foundation for the MR controller, in which the applied current is determined and provided to the MREA by taking the desired MREA force as input.The dynamic MREA models can also be used for simulation analysis of the MR control systems.

Analytical MREA models
The analytical models of MREAs are divided into quasi-static models and unsteady models.The quasi-static models ignore the fluid inertia effect under the assumption that the flow of MR fluid is a quasi-static flow.Contrary to the quasi-static models, the unsteady models take the fluid inertial effect into account, and the MR fluid is considered to be an unsteady flow.The quasi-static model can describe the characteristics of MREAs in low-speed vibrations, but fails to accurately describe the characteristics of MREAs under high-speed impacts.The unsteady model divides the total force of MREAs into parts of MREA force resulting from the viscous flow, the MR effect, the minor losses, the inertia effect, and the spring force for MREAs equipped with an accumulator.The quasistatic model does not consider the inertial force part, and some quasi-static models also ignore the minor losses part.Based on whether the inertia force calculation process assumes the acceleration of the MR fluid at each point within the flow area is the same, the unsteady models can be classified into two categories: models with area-averaged acceleration and models with non-averaged acceleration.Because the flow velocity of the MR fluid increases rapidly when the MREA is subjected to high-speed impact, the inertia effect should not be ignored at this time, however, the quasi-static models ignoring the inertia effect are still the most used analytical models for MREAs.In contrast, there are relatively few studies on the unsteady models of MREAs (Zhang et al 2010, Shou et al 2019).Besides the influence of inertial effects and other factors on analytical models, some analytical MREA models also consider the influence of temperature on the characteristics of MR fluid to achieve more precise control of MREAs (Li et al 2021).

Quasi-static model.
As the impact speed increases, the dynamic range of MREA decreases sharply, and the controllability is severely reduced.This phenomenon is mainly related to the minor losses in the MR fluid gap, and minor losses refer to the additional pressure loss due to the change in the shape of the MR fluid gap.The change in the shape of the MR fluid gap can result in the changes in the flow direction and the flow area of the MR fluid gap, and flow separation, vortices or wake vortices will happen in the fluid.Then the additional pressure loss and the energy loss will be generated, and this is the effect of minor losses.It has a linear relationship between the minor losses and the square of the velocity.In low-speed vibration conditions, the effect of minor losses is negligible.However, in high-speed impact conditions, the passive force of MREA increases greatly owing to the effect of minor losses, while the controllable MREA force remains unchanged, and then the dynamic range of the MREA is significantly reduced.There are many different types of reasons for minor losses, such as entrance or exit effects, gradual/sudden expansions or contractions, bends, valves, and fittings.
Based on the BP constitutive model, Mao et al (2008) established a nonlinear model considering the minor losses for MREAs under high-speed impact conditions, which is known as the BP constitutive model with minor losses (BPM) model.Subsequently, Mao et al (2013) considered the influence of wall roughness on the performance of MREAs, and a new model was developed based on the quasi-static BP model using the Darcy friction coefficient.The model can be expressed as: where F f is the friction force; sgn (•) is the symbolic function; ∆P ml is the pressure drop due to the minor losses; ∆P η , ∆P MR and ∆P ml can be respectively written as: where K m_i stands for the ith minor loss coefficient in the flow system and v i is the average fluid velocity corresponding to this minor loss coefficient.
Mao et al (2014) applied the proposed model to design MREAs.The value of the Darcy friction factor in the active length of the MR fluid gap was considered to be different from that in the non-activated length, and the fluid velocity was also considered to be different.Therefore, the viscous pressure drops of the MR fluid in the active area and the nonactivated area should be calculated separately.Thus, ∆P η can be rewritten as: where f p and f c are the Darcy-friction factors for the fluid flow in the MR gap and the coil gap, respectively; h c is the thickness of the coil gap; v fc is the average fluid velocity in the coil gap.v f and v fc can be calculated as: where A g and A c are the cross-sectional areas of the active MR gap and the effective coil gap, respectively.The Darcy-friction factors, f, is dependent on Reynolds number, Re: where ε is the average roughness of the pipe wall; D h is the hydraulic diameter, whose value is twice the annular valve gap for a parallel-plate channel; α can be calculated by the Reynolds number: The Reynolds number is defined by: where η is the post-yield viscosity.
(17) The roughness Reynolds number is defined by: Then the viscous pressure drop can be calculated.Prototype experiments demonstrate that the model is smoother in the transition region, and the dynamic characteristics of MREAs can be described better in the low Reynolds number region, the high Reynolds number region, and the transition region.It is worth noting that when Mao et al (2014) applied the BPM model to design the MREA, the pressure drop due to the minor losses in equation ( 11) was rewritten as: where K SE and K SC are the coefficients of sudden expansion and sudden contraction, respectively.They were determined by using the empirical formula: In the fabrication of MREAs, the electromagnetic coils are generally encapsulated with epoxy to prevent them from being polished by the MR fluid and worn.The surface roughness in electromagnetic coils will be significantly increased by the epoxy.Considering the characteristics of MR fluids at high flow rates, Singh et al (2014) established a completely rough wall model, in which the Darcy friction factor is a constant independent of the Reynolds number.The characteristics of MREAs can be well described with this proposed model at high Reynolds numbers, and at the same time, the complexity of the model is reduced.The Darcy friction factor was expressed as: 3.2.2.Unsteady model using area-averaged acceleration.
Considering the influence of minor losses and inertia effect, Mao et al (2013) used the area-averaged acceleration to calculate the MREA force due to the inertial effect and established an unsteady BMP model based on the BP constitutive model and the unsteady Bernoulli equation, which can be expressed as: where ∆P inertia is the pressure drop due to inertia effects; A pr is the cross-sectional area of the piston rod; P gas is the pressure of the accumulator.
The pressure drop due to the inertia effect proposed by Mao et al (2013) can be used for the calculation of typical MREAs with multiple flow channels, which can be expressed as: where L ni is the length of the ith passive fluid gap; a f is the instantaneous average acceleration of the fluid flow in the MR valve; a ni is the instantaneous average acceleration of the fluid flow in the ith passive fluid passage in the MREA.a f and a ni are regarded as a variable related to time.
Mao et al (2013) compared the unsteady-BPM model with the previously proposed quasi-static BPM model (Mao 2011), and a series of impact experiments were carried out on an MREA.It is found that both the two models were in good agreement with the prototype experimental results, and the difference between them was very small.In comparison with the quasi-static model, the unsteady-BPM model established by area-averaged acceleration improves the performance very limitedly.
Li (2007) hypothesized that the force of the mixing mode MREA can be considered as a linear superposition of the MREA force in the flow mode and the shear mode.Different unsteady models were established for the mixed modes based on the BP and the Herschel-Bulkley (HB) constitutive models, respectively.Using the control expressions and the boundary conditions in the parallel plate model, the MREA force in the mixed modes was calculated.While the total MREA force was the sum of the two results above, and the prototype experiments were carried out to verify the proposed model.The unsteady model based on the BP constitutive model can be expressed as: where C * B1 , C * B2 , and C * B3 are the correction coefficients of the MREA force, the Coulomb friction force, and the inertia force, respectively, which can be obtained by correcting the experimental data; ẋc is the velocity of the cylinder in MREAs, during the prototype experiment, the piston rod is fixed, and the cylinder moves relative to the piston at the velocity, ẋc ; b is the width of the parallel plates, which is corresponding to the circumference of the inner wall of the cylinder in MREAs.The unsteady model based on the HB constitutive model can be expressed as: where C * H1 and C * H2 are the correction coefficients of the MREA force and the inertia force, respectively; δ = δ/h; δ is a parameter used to represent the width of the post-yield region; δ is the width of the post-yield region; n HB is the flow coefficient of the HB constitutive model; K is the scale parameter of the HB constitutive model.Shou et al (2019) proposed an unsteady model using areaaveraged acceleration to predict the force of the MREA under high-speed impact conditions, and the calculation of the flow field in the annular damping channel was simplified into a parallel plate model.As shown in figure 15, x and y are the abscissa and ordinate, respectively; y 1 is the upper boundary of the pre-yielding region; y 2 is the lower boundary of the preyielding region.
The MR fluid is assumed to be incompressible and flow along the axial direction only, the governing equation can be obtained based on the Navier-Stokes equation and the continuity equation: where u f is the flow velocity of the MR fluid in the parallel plate model and is a function of the ordinate, y, and the time, t; τ is the shear stress; ∆P g (t) is the pressure drop produced by the annular flow channel.The value of ∂u f (y,t) ∂t is assumed to be equal to the average acceleration of the MR fluid in the MR fluid gap, a f : where a p is the acceleration of the piston rod.
The force of the MREA can be expressed as: where The governing equation of MR fluid flow can be expressed as: Also, ∆P ml can be expressed by equation ( 11) and F gas can be expressed as: where P 0 is the initial pressure of the accumulator; V 0 is the initial volume of the accumulator; β is the thermal expansion coefficient of gas (it is a value within the range of 1.4-1.7);x p is the displacement of the piston.

Unsteady model using non-averaged acceleration.
As shown in figure 15, the flow of MR fluid can be divided into three regions.In order to further improve the prediction accuracy of the MREA force under high-speed impact conditions, it is necessary to respectively analyze the three flow regions of MR fluids in the damping channel.Region 1 and Region 3 are post-yield regions, in which the shear stress of MR fluids is greater than the yield stress followed by a plastic flow.While in Region 2, the shear stress is less than the yield stress resulting in a rigid flow.The partial differential equations of the three areas can be expressed as (Shou et al 2019): Region 1: Region 2: Region 3: where υ is the kinematic viscosity; u f1 (y, t), u f2 (y, t), u f3 (y, t) are the flow velocities in regions 1, 2, and 3, respectively; u f1t (y, t) and u f3t (y, t) are the first-order partial derivatives of the flow velocities in regions 1 and 3 with respect to time, respectively; u f1y (y, t) and u f3y (y, t) are the first-order partial derivatives of the flow velocities in regions 1 and 3 with respect to the vertical coordinate y, respectively; u f1yy (y, t) and u f3yy (y, t) are the second-order partial derivatives of the flow velocities in regions 1 and 3 with respect to y, respectively.The velocity distribution of MR fluid in the three regions can be obtained by using equations ( 32)-( 34), and the flow rate of MR fluid in each region can be further calculated.Finally, the pressure drop, ∆P g (t), can be solved by the numerical calculation method, and the total MREA force can be expressed as: To verify the superiority of the unsteady model, a prototype experiment was carried out on a shock test rig.The MREA force predicted by the two unsteady models and the quasistatic model was compared with the measured experimental data.The performance of the model was evaluated from three perspectives: (i) the accuracy of predicting the peak MREA force, the evaluation index is the relative error of the peak force of MREAs; (ii) the agreement between the MREA force curve predicted by the model and the experimental results, the evaluation index is the mean absolute percentage error of the MREA force curve; and (iii) the accuracy of the dynamic range predicted by the model, the evaluation index is the relative error of the dynamic range.
Shou (2020) showed that the unsteady model with nonaveraged acceleration had the best prediction performance under the three evaluation indicators.Therefore, when designing and optimizing the structure of the MREA used in highspeed impact conditions, the unsteady model using nonaverage acceleration should be considered.shown in figures 16(a) and (b), respectively.The model in figure 16(a) can be expressed as: where c 0BH is the damping coefficient; F 0 is the offset in the force included to account for the nonzero mean observed in the measured force due to the presence of the accumulator; F c is the frictional force related to the field-dependent yield stress, which can reflect the controllable force generated due to the MR effect in the MREA.The frictional force, F c , can be taken as a quadratic function of the applied current, I: where a 1 , a 2 , and a 3 are the constants for the controllable MREA force.The parameters, a 1 , a 2 , a 3 , c 0BH , and F 0 need to be identified in this model, and the main drawback is that the characteristics of the MREA cannot be described well at small deformations and low speeds by this model (Bai and Tang 2021).
where c 0BW and c 1BW are the viscous damping observed at low and high velocities; x p0 is the initial displacement of the piston accounting for the effect of the accumulator; k 0BW is the stiffness at high velocities; k 1BW is the accumulator stiffness; y BW is the internal displacement of the MREA; z is the hysteresis output item of the model; α BW , β BW , γ BW , n BW , and A BW are the basic hysteresis parameters, which can be used to adjust the scale and shape of the hysteresis loop.Among these parameters, α BW , c 1BW , and c 0BW , have linear function relationships with the applied current of the MREA.

Restructured model.
Although the modified Bouc-Wen model can effectively describe the nonlinear hysteresis characteristics of MREAs, the structure of this model is relatively complex and there are quite a few unknown parameters.
To solve this problem, Bai et al (2015) proposed a restructured model using the Bouc-Wen operator.The restructured model was built by removing k 0BW and changing the position of c 0BW .
The schematic diagram of the restructured model is shown in figure 18, which can be expressed as: where α BW z represents the hysteresis element; ρ PH and σ PH are the factors related to the hysteresis characteristic.In  comparison with the modified Bouc-Wen model, the number of model parameters was reduced by two with this restructured model.The eight parameters, α BW , c 1BW , c 0BW , f 0 , ρ PH , n BW , σ PH , and k 0BW , need to be identified in this model.Compared with the modified Bouc-Wen model, the prototype experimental results demonstrate that the restructured model can better describe the nonlinear characteristics of MREAs while avoiding the model being too complex, so as to achieve the accurate, fast, and effective control of systems.The RC operator-based hysteresis model can be expressed as: p , when (ẋ p (t) − ẏBW (t)) < 0 ż (t) = 0, when ((ẋ p (t) − ẏBW (t)) = 0 (46) where p > 0 is the hysteresis factor; S is the virtual displacement; g 1 and g 2 , are the tuning functions that are the monotonic functions of the virtual displacement, S; q is the time-varying hysteresis factor; g −1 1 and g −1 2 are the inverse functions of g 1 and g 2 , respectively; S 0 is the reference point of the virtual displacement, S, and x p0 − y BW0 is the displacement when the direction changes at the time, t * .
In equation ( 44), α BW and c 0BW are the functions of the applied current, which can be expressed as: The Bouc-Wen model, the restructured model, and the hysteresis model based on the RC operator were identified by using the multi-island genetic algorithm.The results demonstrated that the difficulty of parameter identification is greatly reduced in the hysteresis model based on the RC operator, and the identification time is only about half of those for the Bouc-Wen model and the restructured model.In comparison with the Bouc-Wen model and the reconstructed model, the hysteresis model based on the RC operator is proved to maintain better prediction accuracy.Shou (2020) showed that the parametric model performs well in predicting the force-velocity hysteresis characteristics of the MREA under sinusoidal excitations but cannot be used for impact excitations.Under high-speed impact loads, the force-displacement curve takes on an irregular shape and no longer approximates an ellipse.Under low-speed vibration conditions, the force-velocity curve of MREAs is similar to a hysteresis loop composed of two monotonous and smooth curves.One of the curves forming the hysteresis loop presents a monotonically increasing law, and the other presents a monotonically decreasing law.There are two inflection points on each of the two curves, and the area enclosed by the hysteresis loop increases with the increase of the applied current.In the case of high-speed impact conditions, the hysteresis loop of the force-velocity curve is composed of two strongly nonlinear non-monotonic curves, and there are an indeterminate number of inflection points on each of the two curves.The relationship between the area of the hysteresis loop and the magnitude of the applied current is not obvious, so the hysteresis loop is difficult to be described by using parametric models composed of elastic elements, damping elements, friction elements, and hysteresis operators.At present, there are not many studies on parametric models of MREAs under highspeed impact conditions.
In equation ( 54), ẋh (I) affects the width of the forcevelocity hysteresis curve.The fitting expressions of the damping coefficient, c 0BW (I), elastic coefficient, k 0BW (I), and ẋh can be respectively defined as: where c IM1 , c IM2 , k IM1 , k IM2 , k IM3 , x IM1 , and x IM2 are the parameters that need to be fitted according to the experimental data.
Xiang et al (2008) demonstrated the good prediction accuracy of the force-time curve of the MREA model, but the prediction accuracy of the force-displacement curve and forcevelocity curve has not been verified.In addition, the impact loads and impact speeds in the experiment were lower than the general high-speed impact conditions.
In order to solve the problem of predicting the force under impact excitations, a common solution is using a semiempirical model or a non-parametric model represented by the neural network model.The semi-empirical models consist of two parts, which are composed of parametric models and nonparametric models, respectively.Since the MR fluid flows fast, the influence of the inertia effect is not negligible, and the inertia effect has an important influence on the dynamic characteristics of the MREA under high-speed impact conditions.When the structure of the MREA is determined, the inertia effect is mainly related to the impact speed and the density of MR fluids.To improve the prediction accuracy at high speeds during the impact, the velocity of the piston should be also introduced as the input velocity in dynamic models, which is also beneficial to solving the problem that most current dynamic models lack the ability to accurately describe the force characteristics at low speeds in a wide range of frequency.

Semi-empirical models.
3.3.2.1.Dynamic RC operator-based hysteresis model.The advantage of parameter models is that they consume fewer computing resources, so they can be easily transformed into inverse models with desired damping force as input and command voltage as output to establish an MREA controller.Parameter models have been successfully used in MREA simulation of low-speed vibration control and the establishment of MREA controller under various working conditions.Semi-empirical and non-parametric models can solve the problem in cases where parametric models cannot describe the strong nonlinear hysteretic mechanical characteristics of MREA under high-speed impact conditions.
To accurately describe the dynamic characteristics of MREAs under impact loads and the force characteristics at low speeds in wide ranges of frequency, Bai and Tang (2021) proposed a dynamic RC operator-based hysteresis model.In the RC operator-based hysteresis model, the variable related to the velocity is introduced into the time-varying hysteresis factor, q, to construct the dynamic RC operator, which can be expressed as: where q 0 is the initial threshold of q and λ RC is the positive coefficient.By introducing a dynamic RC operator, the relationship between the input velocity and the force is established.
The parameter identification of the dynamic RC operator is carried out using the neural network.Firstly, a group of dynamic RC operators is given with hysteresis factors, p, set to different values, and the displacement and the velocity of the piston are input into these operators.Then the hysteresis values output from these operators together with the applied current input into MREAs are used as the input values of the neural network, the final force is output by the neural network as shown in figure 20, where z dRCOiRC with iRC = 1, 2 . . .n is the hysteresis output of the dynamic RC operator corresponding to different hysteresis parameters, a dRCOiRC .
The parameters of this model and the Bouc-Wen model were identified using the same genetic algorithm, and then the performances of the two models were compared with that of the RC operator-based hysteresis model.In all the prototype experiments, the dynamic RC operator-based hysteresis model has the best goodness of fit for the experimental data at any speeds, especially at low speeds.By testing the dynamic characteristics under excitations of different frequencies, the mean-squared error (MSE) of output force is the largest in the Bouc-Wen model, followed by the RC operator-based hysteresis model.As the excitation frequency increases, the MSE of output force in the three models all increase, but the dynamic RC operator-based hysteresis model increases the slowest, which proves the high prediction accuracy of the hysteresis model based on the dynamic RC operator in a wide frequency range.

Non-parametric models.
Non-parametric models commonly used in MREAs include polynomial models, neural network models, fuzzy logic models, and adaptive neuro-fuzzy inference systems.Most non-parametric models still focus on describing the force-velocity hysteresis characteristics of MREAs under low-speed vibration conditions, there are not many studies on non-parametric modeling of MREAs under high-speed impact conditions.Combining the self-learning ability of the neural network and the reasoning ability of the fuzzy system, the adaptive neural-fuzzy inference system (ANFIS) can be used to deal with nonlinear problems of complex systems, which is currently the main solution to establish the dynamic model of MREAs under high-speed impact conditions.

Time-delayed adaptive neuro-fuzzy inference system (TANFIS).
Arsava and Kim (2015) adopted the ANFIS to establish a non-parametric model of MREAs, which was compared with the modified Bouc-Wen model and the modified Bingham model proposed by Xiang et al (2008) and Fu et al (2012), respectively.The prediction results of the ANFIS model are found to be in good agreement with the experimental results.In order to enhance the prediction accuracy of the ANFIS model, the force output by the model at the previous moment is used as a feedback signal and added to the input terminal as an additional input variable to develop a TANFIS, as shown in figure 21.
The displacement, velocity, and acceleration of the piston, the applied current, and the output force at the previous moment are taken as the input of the model, and the membership function value of each input quantity is calculated through the membership function constructed for the fuzzy inference system in the first layer of the network.This is used as the input of the second layer, and the product of all input signals is output from the second layer.The data is regularized in the third layer, and then the node function is applied in the fourth layer to obtain the output of each fuzzy rule, which is summed in the fifth layer as the MREA force output.
The function of the first layer of the network can be expressed as: where µ TA is the appropriate parameterized membership function (MF); x a , y a , z a , and n a are the inputs; A iTA , B iTA , C iTA , D iTA , and E iTA are the MFs; f (t − l) is the output of the systems observed from the previous time steps; O 1 iTA is the membership function value of the input quantity output by different semantic variables.The second layer of the network outputs the product of all inputs, w iTA , which is expressed as: , where iTA = 1, 2, 3, 4. ( 62) The function of the third layer in the network is to regularize the output of the second layer, wiTA , which can be expressed as: , where iTA = 1, 2, 3, 4. (63) The function of the fourth layer of the network, O 4 iTA , can be expressed as: where the value of the node functions, f iTA , can be expressed as: , where iTA = 1, 2, 3, 4.
In equation ( 65), the output of the system at the previous moment, f t−l , is introduced, and l is the time delay which allows the model to use the observations at the previous time step, (t − l), of the system to calculate the current output value; in t iTA is the input at the current moment; e t is the error of time; e (t) is the error between the proposed TANFIS system's output, which is the predict value of the MREA force at the current time, and the measured actual MREA force at the current time; e t−l is the error between the system's output and the measured actual MREA force at previous time step, (t − l).In the first time step, f iTA , can be expressed as: where p iTA , q iTA , k iTA , l iTA , and r iTA are the parameters to be trained respectively.
The function of the fifth layer in the network can be expressed as: (67) Arsava and Kim (2015) demonstrated that the TANFIS model could better predict the nonlinear characteristics of MREAs compared with the common ANFIS models.

Wavelet-based time-delayed adaptive neuro-fuzzy inference model (W-TANFIS).
Arsava et al (2016) built a smart concrete beam equipped with MREAs and conducted drop-induced experiments.Due to the large amount of calculation of the TANFIS model, many design variables need to be considered.The wavelet transform can be introduced to realize data compression and noise reduction.Combining the ANFIS with the wavelet transform, a W-TANFIS model is established to predict the dynamic characteristics of the concrete structure, and the output of the previous moment is still used as a feedback signal.The impact force received by the concrete system and the applied current are taken as input signals, and the acceleration, the deformation, and the stress are taken as output signals to train the established W-TANFIS model, which is shown in figure 22.
In figure 22, the W-TANFIS retains the network architecture of the TANFIS, and only a discrete wavelet transform module is added between the input end and the first layer of the network.Using time and scale window functions, the wavelet transform provides the time-frequency representation of signals.The wavelet transform decomposes a given signal into sub-signals and reconstructs the original signal to compress data and reduce noise (Thuillard 2001).Then the actual test data can be stretched and filtered, and the output signal can be optimized, so as to achieve the purpose of reducing the time for model parameter identification.The continuous wavelet transform can be expressed as: where α walvet and β walvet are the scaling factor and the translation parameter, respectively; ψ is the wavelet function.The continuous wavelet transform can be expressed as: where x iWT (os) is defined as the original signal; X iWT , Y iWT , and Z iWT are the premise variables; s is the scale index; ϕ is the mother function.
Equation ( 69) can be transformed as: where N is the location index.The discrete wavelet transform isolates the high-frequency components from the original signal.In order to investigate both high-frequency and lowfrequency signals, multi-resolution analysis is used to divide the signal into segments.Multi-resolution analysis reduces the required data points by discretizing the function with the step size.The mother function and the corresponding wavelet function can be expressed as: The W-TANFIS model was compared with the ANFIS and the TANFIS models in the prototype experiment, and the results demonstrated that both the W-TANFIS and the TANFIS models predicted the nonlinear characteristics of smart concrete structures well, while the prediction accuracy of the ANFIS model was relatively low.However, the training time of the W-TANFIS model is much shorter than that of the TANFIS model, which is very conductive to the real-time control of the system.Therefore, the W-TANFIS model has the advantage of faster speed and higher accuracy in predicting the dynamic response of MREA under high-speed impact conditions.

Adaptive neuro-fuzzy inference model (FCM-ANFIS)
using fuzzy c-means clustering (FCM) algorithm.The two non-parametric models of MREAs developed based on the ANFIS model in sections 3.3.3.1 and 3.3.3.2 are established by the grid partition (GP) method, which is one of the typical input space partitioning methods for generating inference rules for fuzzy inference systems.However, when the dimensionality of the input space is large, the development of the ANFIS model is very time-consuming due to the need to determine a large number of parameters.In order to reduce the difficulty of parameter identification and improve the real-time response characteristics of the system, Shou et al (2022) developed the ANFIS model for MREAs by utilizing FCM algorithm to divide the input space and compared it with the ANFIS models that divide the input space using the GP method and the subtractive clustering algorithm, respectively.
According to the FCM algorithm, cluster centers are found by minimizing the objective function.Considering a dataset with nFC data points {x 1 , x 2 , …, x nFC } distributed in M-dimensional space is divided into cFC clusters, and each cluster center corresponds to a fuzzy inference rule of the fuzzy inference system.The objective function can be expressed as: where m FC is the weight index; x iFC is the iFCth data point; c jFC is the jFCth cluster center; u iFCjFC is the membership degree of the iFCth data point within the jFCth cluster center.In the execution of the FCM algorithm, it is necessary to normalize the input data, predefine the number of cluster centers, and assign an initial value to each cluster center to generate an initial cluster center matrix, C 0 .Then the membership degree matrix is initialized.The value of the iterative clustering center matrix, C k , and the membership degree matrix, U k , will be continuously updated until the absolute value of U k − U k−1 is less than the criterion value to stop clustering.Then the iterative process is terminated.The fuzzy inference system structure differs when a different cluster number is chosen, and the prediction results of established dynamic models will also be different.The dynamic model corresponding to the cluster number with the best performance is selected as the final training result.
The ANFIS model developed with the FCM algorithm is more suitable for the dynamic model construction of MREAs under high-speed impact conditions.

Control objectives
The control methods for MR shock mitigation control systems have been applied in various fields, including aircraft landing gears, automobile longitudinal impact energy absorbing devices, planetary landers, helicopter seat suspensions, and artillery/gun recoil systems (Han et al 2018).However, the control methods used in most MR shock mitigation control systems are mostly developed based on vibration control methods or general control methods, which cannot optimally satisfy the requirements of shock mitigation.
The research goals for the shock mitigation control are reducing the impact loads transmitted to the protective structure and protecting the structure and the occupants against the damage caused by the impact excitation.An ideal shock mitigation system aims at achieving soft-landing, which means that the piston velocity is reduced to zero when the piston stroke of MREAs is maximized during the process of shock mitigation control.In the coordinate system of figure 23, the soft-landing condition can be expressed as: where D is the maximum piston stroke; t sc is the duration of the compression process.
Wang et al (2021a) simulated a shock mitigation control system using the optimal Bi number control.They demonstrated that if the piston velocity decreases to zero before reaching the maximum stroke, the deceleration peak of the payload is greater than that of the soft-landing case.As a result, the deceleration of the payload may exceed the maximum allowable value, and the payload will be damaged.During the process of shock mitigation control, the shorter the piston stroke is, the greater the peak deceleration of the payload will be.If the payload deceleration exceeds the maximum allowable value during this process, the shorter the piston stroke, the longer the duration of the deceleration exceeding the maximum allowable value.The damage to the payload is not only related to the peak of the payload deceleration but also to the duration when the payload deceleration exceeds the maximum allowable value.The ideal soft-landing control method should ensure that the payload deceleration does not exceed the maximum allowable value, and at the same time minimize the duration of the payload deceleration sustained at the maximum allowable value.Experiments show that under the same impact excitation, the control method that satisfies the soft-landing condition can make full use of the piston stroke so that both the peak value of the payload deceleration and the duration for which the deceleration exceeds the maximum allowable value during the entire process of shock mitigation control is reduced, which improves the protection effect of the shock mitigation control system on the payload.
When the shock mitigation control system is inappropriately designed, there will be a probability that the velocity of the piston does not drop to zero yet even though it reaches the maximum stroke, and the piston will collide with the cylinder of MREAs at this time.As a result, the deceleration curve does not decline monotonically after reaching the maximum value, and the deceleration will suddenly become negative at the moment of impact.A negative deceleration means that the velocity of the payload starts to increase again, and then the deceleration increases to a positive value with a large rate of change.This phenomenon causes the most serious damage to the payload, and the control method is deemed ineffective.Therefore, the deceleration of the protective structure should be decreased according to the monotonically decreasing law after reaching the peak value as the basic control requirement, even when the control objective of the soft-landing cannot be met in the shock mitigation control system.
Current control methods to realize the soft-landing are mainly developed for the SDOF shock mitigation control system shown in figure 23 or the drop-induced shock mitigation control system shown in figure 24.In the SDOF shock mitigation control system, the MREA force control needs to be considered during the rebound stroke, because the protective structure will rebound under the action of the spring after the payload velocity drops to zero.
After the soft-landing was proposed as a control objective, the optimal Bi number control method was proposed and applied to the drop-induced shock mitigation control system and the SDOF shock mitigation control system (Wereley et al 2011, Singh and Wereley 2013).The time delay of MREAs was subsequently considered in the improvement of the optimal Bi number control method (Choi and Wereley 2015).While ensuring that the peak deceleration of the payload does not exceed the maximum allowable value, Wang et al (2021a) proposed the MDDE control method to achieve the control objective of the soft-landing, and the quadratic damping was also considered in the MDDE quadratic (MDDE-Q) control method based on the MDDE control method (Wang et al 2021b).
In addition to a series of control methods based on the optimal Bi number control method, considering the need for shock absorbers of the automobile seat suspension for both low-speed vibration attenuation and high-speed shock mitigation, Bai and Yang (2019) proposed a hybrid controller by combining the Skyhook control used for vibration with the constant force control used for shock mitigation, to achieve the soft-landing control of the payload under shock excitations that occur randomly during the vibration control process.
Compared with the SDOF system, the aircraft landing gears need to consider the absorption efficiency of the impact energy and the control performance of shock mitigation when the aircraft mass and sink velocity are different.The stability needs to be considered for planetary landers when landing on a slope.In shock mitigation control systems for the fixed-wing aircraft landing gears and planetary landers, the control method has not been developed to achieve the soft-landing yet.

Control methods with 'soft-landing' as shock mitigation control objective
4.2.1.Optimal Bi number control.In the drop-induced shock mitigation control system shown in figure 24, Wereley et al (2011) proposed the concept of optimal Bi number to achieve the control objective of the soft-landing.Bingham number, Bi, is defined as the ratio of the controllable MREA force to the passive MREA force.According to this control method, the initial and terminal conditions of the payload velocity and displacement are substituted into the non-dimensional governing equations of the payload displacement, velocity, and acceleration during the control process of shock mitigation, to solve the values of the unknown quantity, Bi, and record it as the optimal Bi number, Bi o .Based on the value of Bi o , the controllable MREA force required for the shock mitigation control process can be calculated as the desired controllable MREA force.In this case, the desired controllable MREA force will be constant during the MREA stroke.
The Bingham number, Bi, can be expressed as: where the passive MREA force is simplified as the product of the passive damping coefficient c and the initial drop velocity ẋe0 ; f MR is the desired controllable MREA force.According to the optimal Bi number control method, the solution process of f MR can be expressed as: where m is the payload mass; g is the acceleration of gravity; F d is the desired MREA force output by the system controller; x o0 is the total available MREA stroke before the impact.Firstly, by substituting the governing equation ( 76) for desired MREA force into equation ( 75), the payload displacement, x o , payload velocity, ẋo , and payload acceleration, ẍo , are obtained as shown in equations ( 77)-( 79).Here, The variables such as non-dimensional payload displacement, velocity, acceleration, and time, t, are given by: Using equations ( 77)-( 79) and the non-dimensional parameters mentioned in equation ( 81), the non-dimensional payload displacement, x o , the payload velocity, ẋo , and the payload acceleration, ẍo , can be expressed as: (84) Substituting the non-dimensional terminal conditions, ẋo ( t) = 0, and, x o = 0, into equations ( 82) and ( 83) to obtain the optimal Bingham number corresponding to the nondimensional velocity, Bi o,ẋo and the optimal Bingham number corresponding to the non-dimensional displacement, Bi o,x o : When Bi o,ẋo = Bi o,xo , the problem of solving Bi o can be transformed into the problem of solving a system of binary linear equations.The time required for the impact event, to , can be obtained by equating equations ( 85) and ( 86), and then equation ( 85) can be solved by replacing t with to to obtain Bi o : The optimal Bi number can also be obtained by using equations ( 85) and ( 87): In equations ( 87) and ( 88), W is the Lambert W function.Substituting the known initial drop velocity, ẋe0 , into Bi o = fMR cẋe0 , the desired controllable MREA force can be obtained in the control process of shock mitigation.

Optimal Bi number control considering rebound control.
Based on the research of Wereley et al (2011), Singh andWereley (2013) proposed an optimal Bi number control for the SDOF shock mitigation system equipped with Based on achieving the control objective of the soft-landing, this method extends the controlling idea of optimal Bi number to the process of rebounding stroke under the action of the spring after the payload compresses the MREA to the maximum stroke.The MREA force of the rebound stroke is controlled to ensure that the load can return smoothly to the position where it is before the impact while avoiding unnecessary oscillations.Li et al (2023a) conducted a prototype experiment on this control method, and a detailed schematic diagram of the control method is illustrated in figure 25.
Singh and Wereley (2013) adopted the opposite datum plane of Wereley et al (2011) to represent the displacement of the payload.The position of the payload when the piston was compressed to the maximum stroke was recorded as the displacement datum plane, and the vertical upward direction was taken as the positive direction.According to the control method, the initial drop velocity, −ẋ e0 , the governing equation of motion for the system, and the boundary conditions of the soft-landing are used to calculate the optimal Bi number in the compression stroke, Bi c o .The state of the payload movement can be determined by comparing the velocities of adjacent sampling points, and the optimal Bi number of the rebound process, Bi r o , can be calculated at the initial moment of the rebound process.
The desired controllable MREA force in the processes of compression and rebound can be determined by solving the optimal Bi c o and Bi r o , and the desired total MREA force can be subsequently calculated.During compression and rebound strokes, the force-tracking module generates an applied current to control the MREA to output the MREA force close to the desired total MREA force.
The desired controllable MREA force, f MRc , during the compression stroke and the desired controllable MREA force, f MRr , during the rebound stroke can be respectively expressed as: where k is the stiffness of the spring.
The core of the optimal Bi number control is to find the optimal Bi number, Bi o .When the optimal controllable MREA force is determined under the current excitation conditions, this optimal controllable MREA force will remain unchanged.When using the optimal Bi number control, the MREA's controllable MREA force remains constant during compression or rebound strokes.The compression process is taken as an example to illustrate the solution method of optimal Bi number, Bi o .
During the compression stroke, the governing expression for the payload displacement can be expressed as: where ζ is the damping ratio of the system; ω n is the undamped natural frequency, ω n = k m , ω d = ω n 1 − ζ 2 ; P 1 and P 2 are the constants.
In order to find the compression stroke, Bi c o , it is necessary to solve P 1 , P 2 , and t sc .The governing equations of the payload velocity and acceleration can be respectively expressed as: The initial conditions of payload displacement and velocity during the compression stroke can be expressed as: Substitute the initial conditions in equation ( 94) into the payload displacement equation ( 91) and payload velocity equation ( 92), P 1 and P 2 can be obtained as: The terminal conditions of the compression stroke can be expressed as: When equations ( 91)-( 93) are non-dimensionally rewritten in the way demonstrated in section 4.2.1, the termination condition of non-dimensional velocity during the compression stroke can be obtained, which will be substituted into the governing equation of non-dimensional payload velocity to solve the non-dimensional compression stroke duration: where ωd , P 2n , and P 1n can be expressed respectively as: 98) 99) where Bi c is the Bingham number corresponding to the MR yield force during the compression stroke, ẋe0 can be expressed as: By rewriting non-dimensionally, the displacement termination condition during the compression stroke in equation ( 96) can be substituted into the expression of non-dimensional payload displacement: where P 1n , P 2n , and tsc contain Bi c o thus equation ( 102) can be expressed as: The numerical method of fixed-point iteration is used to solve Bi c o .When γ Bi is the convergence criteria and the termination condition of iteration, is set as Bi c o , and then f MRc under the current impact excitation can be solved according to equation ( 74).
Using the solution method of Bi c o in the compression process, P 3n , P 4n , and the duration of the rebound process t sr can be solved and expressed as: The initial conditions of the rebound stroke are the terminal conditions of the compression stroke, x o (0) = x o (t sc ) ; ẋo (0) = 0.At the end of the rebound stroke, the payload returns to the equilibrium position, and the terminal conditions of the rebound stroke can be expressed as: x o (t sr ) = 0; ẋo (t sr ) = 0. Using the same non-dimensional analysis method to solve the optimal Bi number during the rebound stroke, the final solution, Bi r o , can be expressed as: (107)

Optimal Bi number control incorporating a time lag.
In sections 4.2.1 and 4.2.2, the optimal Bi number is solved on the assumption that MREAs can immediately provide the required controllable MREA force at the moment of the impact.In fact, the response time of the entire control system is composed of four units: the sensor, the system controller, the MREA controller, and the MREA.There is a time delay in each module between receiving the output from the previous module and the output of this module, and it takes some time for the output signal to reach a stable value from the generation.To improve the control performance of the shock mitigation control system, Choi and Wereley ( 2015) proposed an optimal Bi number control method incorporating a time lag into the drop-induced shock mitigation control system, as shown in figure 24.The control method solves the optimal Bi number, Bi o , under a specific impact condition in a similar way as illustrated in section 4.2.1, and the only difference is the time-delay effect of MREAs is considered in the desired controllable MREA force term in the system governing equation.Further improvement of control effectiveness can be achieved by considering the response time of each module of the shock mitigation system in the controller, but this will increase the complexity of the controller.The governing equation of the system motion can be expressed as: The desired controllable MREA force considering the time delay, f MRt , can be expressed as: where τ a is defined by the time for the step force response to reach 63.2% of its final value in a steady state.Substituting equation ( 109) into equation ( 108) with the initial conditions, the governing equations of payload displacement, velocity, and acceleration with respect to time can be obtained and then rewritten non-dimensionally, respectively.Substituting the termination conditions of displacement and velocity into the non-dimensional displacement-time and velocity-time governing equations, Bi o,ẋo and Bi o,xo can be obtained as: 1Rs −1) where the definitions of τ v , R s , and R v are the same as those in section 4.2.1;τ is the non-dimensional time constant and can be calculated by τ = τ a /τ v .The non-dimensional time, to , can be obtained by solving the following equation: In the process of solving to , it is necessary to approximate the exponential functions in equations ( 110) and ( 111).Using a Taylor expansion, the exponential functions can be approximated by discarding the high-order small quantities: Equations ( 113) and ( 114) are substituted into the equation, Bi o,ẋo ( t) = Bi o,xo ( t), to obtain to , and then to is substituted into equation ( 110) to obtain Bi oa .The desired controllable MREA force, f MR , calculated by the controller can be obtained from where In order to further improve the solution accuracy of optimal Bi number, Bi oa , needs to be corrected: where Bi oa is corrected by equation ( 118) to obtain Bi o .α Bi and β Bi coefficients determined as a function of the nondimensional time constant: where α Bi1 , α Bi2 , α Bi3 , α Bi4 , β Bi1 , β Bi2 , β Bi3 , and β Bi4 are parameters that need to be identified.
In the computer simulation experiment, benefiting from the fact that the simulation computer has more powerful computing power than the controller of the shock mitigation control system, the computer can utilize the graphical method with higher precision to solve equation ( 120) and determine the exact value of Bi o , which is denoted as Bi oexact : Simulation results show that when the impact speed is high, only Bi o calculated by the optimal Bi number control incorporating a time lag can follow the trend of accurate fit.
At this time, the control method that ignores the time lag cannot achieve the control objective of the soft-landing.However, using the control method described in this section, the calculated value of Bi o is basically consistent with the value of Bi oexact .The control performance of the optimal Bi number control can be improved in high-speed impact conditions by considering the time lag.Based on this control method, further consideration of the time lag in the entire shock mitigation control system will become a direction in the research field of MR shock mitigation control in the future.

MDDE control.
Optimal Bi number control methods mentioned in sections 4.2.1-4.2.3 have been successfully applied to shock mounts, automobile suspensions, sliding seat systems for ground vehicles, helicopter landing gears, and other fields.Wang et al (2021a) proposed the MDDE control method to make full use of the energy absorption capacity of MREA and achieve occupant protection at higher impact speeds than the optimal Bi number control method, which is composed of the MCD controller and the optimal Bi number controller.
The soft-landing control can be achieved by the optimal Bi number control method when the initial velocity of impacts is lower than the maximum allowable velocity.However, when the initial velocity of impacts exceeds the maximum allowable velocity, the optimal Bi number control method cannot achieve the soft-landing objective.The MCD controller is used to keep the deceleration of the payload at the maximum allowable value, and then the piston velocity is rapidly reduced.Afterwards, the optimal Bi number controller is used to achieve the soft-landing.
Switching between the two controllers enables the shock mitigation system to achieve a soft-landing at a higher impact speed while minimizing the duration of the payload deceleration at the maximum allowable value.In the control process, the total stroking time of the MCD control section, t 1 , and the output of the optimal Bi number controller, Bi o , are taken as variables of the controller to be solved.
Wang et al (2021a) took the equilibrium position of the payload before being impacted as the reference plane for the vertical displacement, and the vertical downward direction was taken as the positive direction.The values of the desired controllable MREA force, f MR , when the MCD controller works, can be solved: When the maximum allowable payload deceleration, ẍomax , is known, f MR can be obtained by solving equations ( 121) and ( 122) simultaneously.Supposing the total stroking time of the MCD control section is t 1 , when the optimal Bi number controller starts to work, the payload displacement, x o1 , and velocity, ẋo1 , in the vertical direction can be expressed as: Then the optimal Bi number controller starts working.Using the optimal Bi number control method in section 4.2.1 to solve the governing equations of the payload displacement, x o (t), and velocity, ẋo (t), in the vertical direction: where The stroking time of the optimal Bi number controller is supposed to be t 2 , and the soft-landing displacement and velocity termination conditions in equation ( 73) are substituted into equations ( 125) and ( 126), respectively: where t 1 and t 2 can be calculated with equations ( 128) and ( 129), and t 1 can be used to calculate the desired controllable MREA force, f MR1 , when the MCD controller finishes working: Substituting t 1 into equation ( 130), the desired controllable MREA force, f MR1 , is obtained when the MCD controller finishes working.The value of the desired controllable MREA force in the working process of the optimal Bi number controller is kept as f MR1 .The controller outputs the desired controllable MREA force according to equation ( 122) during the period from zero to t 1 , and maintains the value at f MR1 during the time period from t 1 to t 2 .When the initial velocity of impacts is low, t 1 is zero, and only the optimal Bi number controller is used in the MDDE controller to output the desired controllable MREA force.

MDDE-Q control.
In previous studies, the passive MREA force, F η (t) , was considered as a proportional function of the piston velocity.However, under high-speed impacts the Reynolds number of the MR fluid flowing in the MR fluid gap is greater than 2000.The MR fluid changes from laminar flow to turbulent flow, and the passive MREA force of MREAs has quadratic damping characteristics.Wang et al (2021b) regarded the passive MREA force, F η (t), as a quadratic function of the piston velocity, and proposed the MDDE-Q control method following the same way as demonstrated in section 4.2.4.
The passive MREA force can be expressed as: where c 1 and c 2 are the quadratic and viscous damping constants, respectively.Considering the quadratic damping characteristics, the desired controllable MREA force of the MCD controller can be expressed as: When the total stroking time of the MCD controller is supposed to be t 1 , and the optimal Bi number controller starts to work, the velocity, ẋo1 , and the displacement, x o1 , of the payload in vertical direction can be expressed as: Introducing the quadratic damping term, the governing equation of the system motion can be expressed as: where τ v = m c1 ; v1,2 can be expressed as: where v1, 2 will be solved in the subsequent calculation process.
Supposing the stroking time of the optimal Bi number control section is t 2 , the soft-landing termination conditions are substituted into the governing equations: The optimal value, vo 1,2 , of v1,2 can be calculated by equations ( 137) and ( 138).The value of f MR obtained by substituting vo 1,2 into the equation ( 136) is the desired controllable MREA force, f MR1 , when the MCD controller finishes working.
The critical yield force of the MREA is substituted into equation ( 136) to solve the value of v1,2 .Then, v1,2 and terminal conditions are substituted into the governing equations of the payload displacement and velocity of the optimal Bi number controller.The critical initial velocity is sequentially obtained.
When the initial velocity of impacts is smaller than the critical initial impact speed, the MCD controller does not work, while the optimal Bi number controller provides a constant desired controllable MREA force during the entire impact event.When the impact speed is higher than the critical initial impact speed, the MCD controller outputs the desired controllable MREA force until t = t 1 , and then switches to the optimal Bi number controller to output the desired controllable MREA force of f MR1 until the shock mitigation control process ends.
Compared with the optimal Bi number control, MDDE-Q control method can increase the maximum allowable initial drop velocity of the payload by keeping the deceleration of the payload and the total MREA force at the maximum allowable values.The initial drop velocity at which the payload reaches the maximum allowable deceleration under passive MREA force is the maximum allowable initial drop velocity of MDDE-Q control method.

Hybrid control for both shock and vibration mitigation.
For special vehicles, the seat suspension should not only control the vibration caused by uneven roads to improve the comfort of the occupants, but also mitigate the impact of the explosion of landmines and explosive devices to reduce the damage to the occupants, which requires a control method with the good ability of both vibration control and shock mitigation.Bai and Yang (2019) developed a hybrid controller for seat suspension, as shown in figure 26.The controller is composed of a Skyhook controller for vibration control, a soft-landing controller for impact conditions, and a switcher for judging the current state to switch the control.The controller takes the velocity and displacement of both the vehicle floor and the seat as input quantities, and the static equilibrium positions of the floor and the seat are used as the reference planes, respectively.The upward direction is taken as the positive direction of both velocity and displacement.The vehicle floor is considered as the excitation source, and the seat is taken as the payload in the system.
When it is judged that the system input excitation is vibration, the Skyhook controller outputs the desired controllable MREA force, F MRsky : where C sky is the Skyhook damping; F d-vibration is the desired MREA force of the vibration controller.
When it is judged that the system input excitation is impact, the shock mitigation control method is used to calculate the desired seat acceleration and the desired total MREA force, F d , at this stage.The entire process of shock mitigation control is divided into three stages: compression, rebound, and second compression.By solving the kinematic equations, the desired seat acceleration, a o , that the seat should bear to achieve the control objective of the soft-landing can be obtained at each stage.
Based on the solution of the desired seat acceleration, the difference, ∆, between the current external force, mẍ o , and the desired external force, ma o , needs to be solved for the seat.Then the actual total MREA force output by the tracking module of force is corrected using the difference, ∆, and this process can be expressed as: where F d-shock is the desired MREA force of the shock controller.
Since control objectives and corresponding control performances of vibration and shock controllers are different, a set of judging logic is introduced to select and switch between vibration and shock controllers.The judgment process is divided into four parts.(Bai and Yang, 2019).
Part 1: The states of systems are obtained including the displacement, x e , and velocity, ẋe , of the vehicle floor considered as the excitation, the displacement, x o , and velocity, ẋo , of the seat taken as the payload, the desired total MREA force calculated by the vibration controller, and the total MREA force calculated by the shock controller.The vehicle floor is considered as the excitation source of this system.In the initial state, the flag of the static variable used for the selection of vibration or shock controller is set to 0, and the controller works in the state of vibration control.
Part 2: The occupant and seat are regarded as the payload with unknown mass, m, and the displacement, velocity, and acceleration of the seat are regarded as the displacement, velocity, and acceleration of the payload.According to the velocities and displacements of the floor and the seat measured by the sensors, the solution process of m can be expressed as: Part 3: By using the kinematic equations to calculate circularly, the displacement of the seat is solved when the seat velocity is higher than or equal to the floor velocity under the action of the vibration controller.Denoting the to-be-calculated value of the seat displacement as Z, and the kinematic equation can be expressed as: where ∆h is the time step of program execution.The piston stroke required for the seat velocity to be increased to the floor velocity can be determined in the system as x e − Z.The seat will collide with the floor by using a vibration control strategy, once the piston stroke required x e − Z is greater than the currently available piston stroke, D tr .Then the flag is set to 1, and the controller will switch to the shock controller.
Part 4: When the sensors detect the floor velocity, ẋe , which is also the excitation velocity, is less than the minimum value, σ, set by the system, the floor velocity can be regarded as zero, and there is no impact excitation input in the system at this time.The default value of the flag is set to zero.When the shock mitigation control ends, the controller is switched to the vibration controller.
Simulation results show that, under vibration excitation, in comparison with the passive control and the constant current control, the vibration control performance of the hybrid control system improves significantly in the low frequency region, the resonance region, and the isolation region.Also, the system can achieve the control objective of soft-landing under impact excitation and can switch between vibration control and shock control on time according to working conditions.When multiobjective control is desired, such as the maximum deceleration threshold or vibration control and shock mitigation control effects, it can be considered to combine two or more control methods through reasonable switching algorithms.Although the system complexity as well as the cost are increased.

Skyhook control.
Typical controllers including Skyhook and PID controllers are widely used in various vibration controls based on MR fluids and MR elastomers and have been proven to be a simple and low-cost solution in the field of vibration control.At present, several research teams have applied the Skyhook controller to the MR shock mitigation control system.Han et al (2018) designed and manufactured an MREA for the aircraft landing gear, and a modified Skyhook controller was adopted in this MREA.The dynamic model of the MR aircraft landing gear system is shown in figure 27.The controllable MREA force of the Skyhook controller, F MRsky , can be expressed as: where ẋo is the sinking velocity of the airplane mass; ẋe is the descending velocity of the landing gear wheel.
In the shock mitigation control system of aircraft landing gears, the aircraft mass is used as the payload, and the impact excitation is input to the aircraft mass from the wheels.To ensure a smooth landing process, a larger total MREA force is required to absorb the impact energy when the aircraft mass compresses the MREA, while a smaller total MREA force is required to keep the tire grip in the rebound process by increasing the contact between the tire and the ground.When the aircraft mass and the wheels move in opposite directions, the MREA is in the rebound state, and the controllable MREA force provided by the MREA is zero.When the airplane mass and the wheels move in the same direction if the sinking velocity of the aircraft mass is greater than the velocity of the wheel, MREA needs to provide a larger MREA force to attenuate the vibration of the aircraft mass, and the output controllable MREA force is C sky ẋo ; if the sinking velocity of the airplane mass is lower than the velocity of the wheel, the smaller total MREA force helps to reduce the vibration transmitted from the wheel to the aircraft mass, and the controllable MREA force provided by the MREA is zero.
For the system of aircraft landing gears, the energy absorption efficiency is often used as an evaluation index of the shock mitigation control performance, and the peak acceleration transmitted to the aircraft mass can be used as another evaluation index to evaluate the shock mitigation control performance.Energy absorption efficiency refers to the ratio of the area enclosed by the MREA force-displacement curve to the product of the maximum MREA force and the maximum piston stroke, which reflects how much the control method can make full use of the energy absorption capacity of MREAs.Energy absorption efficiency, η efficiency , can be expressed as: where F totalmax is the peak of the MREA force.The higher the energy absorption efficiency, the lower the impact energy transmitted to the payload, and the better the protection effect on the payload.The 1/2 aircraft model is established for simulation analysis.The simulation results demonstrate that compared with the constant current control, the Skyhook control can improve energy absorption efficiency and reduce the peak acceleration of the aircraft mass peak by reducing the peak value of the MREA force.Bai et al (2017) established a SDOF shock and vibration control system by using the Skyhook control method, as shown in figure 28.The peak acceleration was taken as the evaluation index of the impact mitigation control performance, which also verified the effectiveness of the Skyhook control method.
For the aircraft landing gear shock mitigation control system and the SDOF shock and vibration control system, experimental results have shown that the Skyhook control method can effectively improve the shock control performance compared with the control method with a fixed applied current.The defect of the Skyhook control method is that there are different optimal gains for different impact excitations.Taking the aircraft landing gear as an example, due to the influence of remaining fuel, number of passengers, and weather conditions, the weight and the sink velocity of the aircraft lie in a certain range, and different optimal gains corresponding to different landing conditions will fluctuate accordingly.The conventional Skyhook controller cannot automatically adjust its gain to different impact excitations, thus new control methods are needed to enhance the shock mitigation control performance and finally achieve the control objective of the soft-landing.

PID control.
Similar to the Skyhook control method, as a simple and effective control method, the PID control is widely used in industrial automation.The PID control applied in the MR vibration and shock mitigation control system can be expressed as: where u PID is the desired MREA force after correcting the deviation (Choi et al 2016).Equation ( 147) indicated that the reasonable selection of K P , K I , and K D for the PID controller plays an important role in the control performance.The values of K P , K I , and K D are mainly determined by the trial-anderror method, and there is no method to prove mathematically whether the values of K P , K I , and K D values are the optimal selection.Hu et al (2012) applied PID control to an MR gun recoil system, and the maximum stroke of the piston and the peak MREA force were chosen as the evaluation index of the shock mitigation control performance.Four sets of K P , K I , and K D values were selected for prototype experiment, the significant difference in the stroke of the piston and peak values of the MREA force indicates that the reasonable selection of PID control parameters has a significant impact on control performance.They demonstrated that the PID controller can provide better control performance than the simple on-off control or the passive control.However, further improvement in control performance needs the application of advanced modern control methods, such as fuzzy control.

Hybrid control with Skyhook controller.
Han et al (2019) adopted a hybrid controller on the aircraft landing gear to enhance the effect of shock mitigation control and compared it with the Skyhook controller.The control schematic of the hybrid controller is shown in figure 29.
The governing equation of the desired controllable MREA force, F MRhy , can be expressed as: where F MRref is the reference force; F MRforcect is the controllable MREA force generated by the force controller and can be expressed by F MRref − F MRsky .Developing this control method is mainly to further improve the efficiency of energy absorption during the shock mitigation control based on the Skyhook control method.In the compression process, the desired controllable MREA force, F MRsky , will quickly reach the first peak with the increase of piston displacement, decrease with the further increase of displacement, and then rise sharply to the second peak.The sharp change of the desired controllable MREA force reduces the area enclosed by the MREA forcedisplacement curve, thereby reducing the efficiency of energy absorption.
When the piston displacement is small, the force controller does not work and the Skyhook control method is used and the desired controllable MREA force equals to F MRsky .As the piston displacement increases, the force, F MRsky , also increases.When F MRsky reaches the first peak point, the value of the reference force, F MRref , is set as the first peak value of F MRsky .Then, the force controller and the Skyhook controller start to work and the desired controllable MREA force is set and kept equal to the value of F MRhy .The desired controllable MREA force follows F MRhy , so that the desired controllable MREA force remains equal to the first peak point of F MRsky .This makes the MREA force-displacement curve more rectangular to improve the efficiency of energy absorption.
Han et al (2019) demonstrated that compared with the Skyhook controller, the hybrid controller could increase the efficiency of energy absorption of impact mitigation control while reducing the peak value of total MREA force.

Control method for 100% energy absorption efficiency.
In comparison with the Skyhook control method, Yoon et al (2020) proposed a control method based on the conservation of mechanical energy to improve the energy absorption efficiency and further reduce the rate of change of the MREA force during landing.This control method can achieve good performance of shock mitigation control under various landing conditions without adjusting the gain.
The sum of potential energy and kinetic energy at the moment of landing equals to the energy absorbed by the landing gear and tires during landing, and the equation can be expressed as: m e is the mass of the landing gear; v sink is the sink velocity of the aircraft each time; the difference between the displacement of the aircraft and the wheel, (x o − x e ), is the quantity to be solved and recorded as x end , which means that the piston displacement at the end of the first compression to fully absorb the impact energy input during landing at the current sink velocity; k t is equivalent stiffness of the wheel.
Using x end to solve the desired MREA force, F d , the solution process can be expressed as: Substituting x end into equation ( 150), the desired MREA force, F d , is obtained to completely absorb the impact energy input when the sink velocity is v sink , and the desired MREA force, F d , remains unchanged during the entire compression process.Because the spring force generated by the accumulator increases with the increase of piston displacement, the desired controllable MREA force decreases with the increase of piston displacement.Yoon et al (2020) showed that, under various simulation conditions, the control method maintained an energy absorption efficiency of more than 90%, and the energy absorption efficiency could reach 97% under specific conditions.Compared with the Skyhook control method, the efficiency of energy absorption during landing has been increased by 16%, and the derivative of the payload acceleration with respect to time has been reduced by 42%.Thus, the performance has been significantly improved in the shock mitigation control system of aircraft landing gears.

Adaptive sliding mode hybrid control.
Luong et al (2020) proposed an adaptive hybrid controller in sliding mode to achieve good performance of shock mitigation control for aircraft landing with different aircraft masses and sink velocities.The adaptive hybrid controller in sliding mode calculates the desired total MREA force according to the sink velocity and aircraft mass and reduces the deviation between the total MREA force provided by MREAs and the desired total MREA force.The adaptive hybrid controller in sliding mode consists of a sliding mode controller that outputs the desired controllable MREA force and a reference model that generates a reference value for the displacement of the aircraft mass, as shown in figures 30 and 31, respectively.
Estimating the mass, m, of the aircraft landings as: Combining with the kinematics control expression of the aircraft, equation ( 151) is rewritten as: where ṁ can be expressed by m and ẍo .The Lyapunov function is established using equation ( 152) to determine the convergence of the aircraft mass: The reference module consists of a Skyhook controller and a force control module, and the schematic diagram is shown in figure 31.The relationship between the passive force of MREAs and the piston stroke, aircraft mass, and sink speed is measured through experiments which can be used to calculate the target total MREA force, F target , and then the controller will make the total MREA force track the target total MREA force.The target total MREA force can be solved according to the relationship between the passive force of MREAs and the piston stroke: where x pmax,passive is the maximum piston stroke during the process of shock mitigation control with no current input.When the desired total MREA force is greater than the passive MREA force, the Skyhook controller determines the difference between the desired total MREA force and the passive MREA and the Skyhook gain can be expressed as: where ẋo | ẋp= ẋpistonmax is the maximum velocity of aircraft mass; F η | ẋp= ẋpistonmax is the passive MREA force at the maximum velocity of aircraft mass.
The desired controllable MREA force of the Skyhook controller can be expressed as: Figure 31 shows the reference module for calculating the history of the reference signal, x d o .The analytical model with an adaptive hybrid control algorithm is executed before the landing gear touches down, and this progress depends on the estimated mass and the sink velocity.The difference between the displacement, x o , measured by the sensor, and the reference signal, x d o , is calculated as the tracking error, xo , and the process is shown in figure 30.
Based on the desired controllable MREA force of the Skyhook controller, the sliding mode controller is used to correct the desired controllable MREA force.An estimate of the desired controllable MREA force, FdMR , is calculated with an estimate of the mass, m, and an estimate of the total MREA force, Ftotal , as: The corrected desired controllable MREA force can be expressed as: where Sat(•) is the saturation function; ξ f = 0.1 is the boundary thickness; k f is the discontinuous control gain, which is selected to satisfy the inequality in equation ( 159), the sliding surface, S f can be obtained according to equation ( 160): where λ f is the slope of the sliding surface line.Luong et al (2020) selected the peak actual total MREA force, the maximum piston stroke, and the energy absorption efficiency as evaluation indicators.By using the adaptive hybrid control method in sliding mode, computer simulation analysis proved that both the maximum stroke and the peak total MREA force are reduced under various working conditions and the energy absorption efficiency is improved.When the aircraft mass m o is 680 kg and the initial sink velocity is 3.5 m s −1 , the peak total MREA force is reduced from 27.9 kN of the Skyhook controller to 25.5 kN, the maximum stroke of the piston is reduced from 0.224 m to 0.220 m, and the energy absorption efficiency is significantly increased from 86.4% to 96.4%.

Control method based on Lyapunov function.
The lander used to detect extraterrestrial planets faces the risk of impact and overturning while landing on the surface of the planet.Therefore, it is necessary to improve the impact energy absorption capacity and the landing stability of the landing gear.Incorporating MREAs is a way to further improve the performance of the landing gear.Maeda et al (2019) developed a control method for the MREA semi-active landing gear system of a planetary lander to avoid capsizing.Aiming at minimizing the derivative of the Lyapunov function with respect to time, the derived control method based on the 2DOF model of the lander shown in figure 32 can be expressed as: where c R and c L are the damping coefficients of the left and right MREAs, respectively; ẋpR and ẋpL are the piston velocities of the left and right MREAs, respectively; ω is the roll rate of the lander and takes the counterclockwise direction as the positive direction.To minimize the derivative of the Lyapunov function with respect to time, the control method minimizes the force of the landing gear on the side that touches the ground first and maximizes the force of the landing gear on the other side that touches the ground later.
In figure 32, the force, f a , is exerted at the point of the lander via the landing gear from the terrain.Since the damping of the right landing gear that first touches the slope is reduced to the minimum, the force, f a , is lower than that of the passive honeycomb aluminum structured landing gear, and the contraction of the right landing gear is increased to suppress the change of the attitude of the lander during landing.When the left landing gear subsequently contacts the slope surface, because the damping is adjusted to the maximum value, the force, f a , increases in comparison with the passive landing gear and reaches the maximum value in the entire control process.At this time, increasing the force, f a , is beneficial to increasing the moment of suppressing rollover and preventing the lander from rolling over.
The landing process of the lander was simulated by computer, and the prototype experiment was carried out, compared with the passive MREA, the prototype experiment demonstrates that the semi-active landing gear system with the proposed control method can reduce the peak roll rate and the vertical acceleration of the lander when one side of the landing gear contacts for the first time, which solves the rollover problem of passive honeycomb aluminum structures.

Other control methods.
The controllers based on fuzzy logic do not require precise analytical models and have strong robustness to deal with nonlinearity and uncertain systems.The fuzzy control has been successfully applied in both low-speed vibration and high-speed shock mitigation control systems based on MREA.
On the basis of the research of Maeda et al (2019) on a lander equipped with MREA landing gear, Wang et al (2019) proposed a new method using the bypass MREA as the main structure of the landing gear and discussed the feasibility of the landing gear structure on a heavy-duty lander.Unlike Maeda et al (2019) who controlled the MREA independently, this fuzzy controller simultaneously and systematically controls the MREAs on four landing gears.The acceleration, the derivative of acceleration with respect to time, the pitch angle, and the roll angle of the cabin are defined as input values; and the voltage values of the four MREAs are defined as output values in this controller, which allows the landing gear to be adapted to different landing conditions.
When designing a fuzzy controller, it is necessary to first determine the respective control ranges of input and output, which is the discourse domain.Subsequently, set appropriate linguistic variables for each variable.Secondly, select the corresponding membership function for each linguistic variable.Finally, establish the fuzzy control rule to determine the corresponding relationship between input and output linguistic variables.
The fuzzy control rules selected by the fuzzy controller are: (i) when the acceleration is increasing and less than the set value, and the derivative of the acceleration with respect to time is positive, a low voltage is output; (ii) when the acceleration increases and is greater than the set value, and the derivative of acceleration with respect to time is positive, the output voltage should be zero; (iii) when the acceleration is decreasing and more than the set value, and the derivative of acceleration with respect to time is negative, the output voltage is zero; (iv) when the acceleration is less than the set value, and the derivative of the acceleration with respect to time is negative, the output voltage is high.The core of the fuzzy control rule is to ensure that the acceleration is as close as possible to the set value.In addition, the output voltage should also be affected by the pitch angle and the roll angle of landers to reduce the acceleration during landing, considering the precondition that the lander does not roll over.
Simulation results show that compared with the landing gear with a passive honeycomb aluminum structure, this landing gear can reduce the piston stroke while significantly reduce the peak deceleration under experimental conditions.The ability of the lander to land on terrains with large slopes and low friction coefficients has also been significantly improved.
In the MR shock mitigation control system, fuzzy control is often used in conjunction with other control methods.Li et al (2019) proposed an optimal control with a fuzzy compensation method for MREAs in full-scale gun recoil control.The optimal control needs to solve the desired MREA force in the process of shock mitigation control using the analytical model of MREAs.There are certain errors in the analytical model in describing the nonlinear characteristics of the MREA in the system at high speeds, for example, in sections 4.2.5 and 4.2.6 the items of passive MREA force in the total MREA force are regarded as a linear function and a quadratic function of the piston velocity, respectively.Thus, the accuracy of the optimal control will be affected by the accuracy of the analytical model.On the contrary, the output of fuzzy control is only related to the experimental data, so its accuracy is only affected by the measurement error of sensors.A reasonable combination of the two methods can achieve more accurate control of the total MREA force.
In this control scheme, the main part of the desired MREA force is output by the optimal controller and corrected by the fuzzy controller, which reduces the adverse effect on the control performance caused by the lack of precision in analytical models or the measurement error of sensors.Although the fuzzy controller cannot work when the sensor fails, the optimal controller can still provide the basic desired MREA force for the system.The control schematic diagram of the controller is demonstrated in figure 33.
According to figure 33, the governing equation of the desired controllable MREA force of the optimal controller, F MRoptimal , is: where F R is the ideal MREA force; F k0 is the initial spring force.The ideal recoil force can be obtained with the softlanding conditions: where ẋo0 is the initial velocity of the barrel.
The barrel displacement and the barrel velocity are taken as inputs for the fuzzy controller.The range of the barrel velocity is 0-4 m s −1 , and the range of the barrel displacement is 0-0.2 m.Accordingly, the discourse domain of the barrel velocity is set as [0, 4], and the discourse domain of barrel displacement is set as [0, 0.2].Using more linguistic variables of input or output can improve the accuracy of the controller which will lead to better buffering performance of the gun recoil system.The input linguistic variables are designed as seven subsets, which are denoted by zero small (ZS), zero medium (ZM), zero big (ZB), medium (M), big small (BS), big medium (BM), and big (B).The discourse domain of the output current is set as [0, 2].To obtain a more accurate output value, the output linguistic variables are designed as nine subsets, which are denoted by ZS, ZM, zero large (ZL), medium small (MS), medium (M), medium large (ML), large small (LS), large medium (LM), and large (L).The triangle-shaped function is chosen as the membership function of the inputs, and the Gaussian function is chosen as the membership function of the applied current.
The error of the desired total MREA force output by the optimal control gradually accumulates over time.The control rule can be expressed as: in the initial stage of the recoil process when the barrel displacement is small and the barrel velocity is high, the output compensation value, ∆F MRfuzzy , should be small; as time goes by, the barrel displacement increases, the barrel velocity decreases, and the output compensation value, ∆F MRfuzzy , should increase gradually.

Evaluation indexes for shock mitigation control performance
When designing the MR shock mitigation control system, it is necessary to evaluate the control performance.The commonly used evaluation indicators include energy absorption efficiency, the derivative of the payload acceleration with respect to time 'Jerk', V-ACR, AVCR, transient transmissibility, and the energy dissipation ratio.

Energy absorption efficiency.
The energy absorption efficiency refers to the ratio of the integral of the MREA force F total with respect to the piston displacement, x p , to the product of the maximum value of the MREA force, F totalmax , and the maximum piston stroke, D, during the compression process of the shock mitigation control.It is defined in equation ( 146).
The energy absorption efficiency reflects the ability of MREAs to absorb the impact energy (Yoon and Wereley 2020).The smaller the variations of the MREA force caused by different piston strokes, the more rectangular the MREA force-piston displacement curve.Then the higher the efficiency, η efficiency , the lower the impact energy transmitted to the payload, and the better the protection effect.

Derivative of payload acceleration with respect to time
Jerk.The derivative of the payload acceleration with respect to time is related to the riding comfort of passengers, so it is often used as an evaluation index for aircraft landing gears and automobile suspensions.It is defined as: ...
x o is the derivative of the payload acceleration with respect to time.The smaller the variation of the MREA force as time passes, the smaller the value of the Jerk.The control performance of the system can be evaluated by drawing the time-Jerk curve (Yoon et al 2020).
When optimizing the shock mitigation control system, both the peak and the average values of the Jerk should be reduced.Combined with section 4.3.1, the control system should minimize the variation of MREA force caused by different piston strokes as time passes and maximize the energy absorption efficiency, and at the same time the peak and the average values of Jerk should be reduced.

V-ACR and AVCR.
As two evaluation indexes, V-ACR and AVCR were proposed by Li et al (2023a) to characterize the stability and the efficiency of the shock mitigation control process, respectively.The V-ACR ratio is defined as the ratio of the payload acceleration to the initial velocity of the shock excitation: When using V-ACR to evaluate the control performance, the ratio of the maximum value to the minimum value during the shock mitigation process is usually considered.The smaller the ratio is, the smoother the shock mitigation process will be.
The AVCR is defined as: where ẋomax and ẋomin are the maximum and minimum payload velocities during the shock mitigation control process, respectively; x pu is the stroke utilized by the piston during the shock mitigation control process.The efficiency of the shock mitigation process is proportional to the AVCR.

Transient transmissibility and energy dissipation ratio.
Mao et al (2023) developed a sliding seat system based on MREA to reduce the impact loads transmitted to the occupant in the event of longitudinal collisions of the vehicle.The transient transmissibility and the energy dissipation ratio were induced to evaluate the protection effect of the sliding seat system on the occupant.The transient transmissibility is used to quantitatively describe the relationship between the deceleration transmitted to the payload and the deceleration of the input impact during the process of shock mitigation.Two forms of transient transmissibility were formulated, one is called the peak transient transmissibility (PTT) and the other is called the average transient transmissibility (ATT).The PTT is defined as the ratio of the peak seat deceleration to the peak floor deceleration: The ATT is defined as the ratio of the time averaged seat deceleration, ẍo-average , to the time averaged floor deceleration, ẍe-average , which is expressed as: The smaller the values of PTT and ATT, the smaller the proportion of the impact transmitted to the payload, and the better the control performance of the shock mitigation control system.
The energy dissipation ratio is defined as the ratio of the energy dissipated by the MREA during the shock mitigation process to the initial kinetic energy of the payload (including the occupant body mass and the seat): where t s and t e are the moments when the shock mitigation process begins and ends, respectively.The higher the energy dissipation rate, the stronger the ability of the system to absorb the impact energy, and the better the control performance of the shock mitigation control system.
4.4.5.Other evaluation indicators.In the shock mitigation control process, both the peak value of the MREA force and the piston stroke are used to compare the performances of different control methods in the same shock mitigation control system under the same excitation.For example, Han et al (2019) compared the proposed hybrid control method with the Skyhook control method in the simulation process.They showed that the energy absorption efficiency increased from 88.2% to 94.2%, the peak value of the MREA force decreased from 22.47 kN to 21.08 kN, and the piston stroke decreased from 0.2230 m to 0.2216 m, which are used as the basis for the effectiveness of the proposed hybrid control method based on the Skyhook control method.When adjusting the PID controller parameters of the gun recoil system, Hu et al (2012) also used the reductions in the peak MREA force and the piston stroke as the rule for adjusting the parameters.After the parameters of the PID controller were selected, the peak MREA force and the piston strokes were used to compare the PID controller with the on-off and the fuzzy controllers.
To evaluate the control performance of the gun recoil system, Li et al (2019) proposed the relative time from the start of the recoil process to the time when the MREA force first drops to 85% of the peak MREA force, χ 1 , the relative displacement from the start of the recoil process to the displacement where the MREA force first drops to 85% of the peak MREA force, χ 2 , and the relative variation of the MREA force from the peak MREA force to the end of the recoil process χ 3 .Among them, the larger the values of the relative time, χ 1 , and the relative displacement, χ 2 , the smaller the relative variation of the MREA force, χ 3 , and the better the shock mitigation control performance.
In addition, the peak deceleration of the payload is also an important index for evaluating the shock mitigation control performance.Choi et al (2005) proposed a control method for the shock mitigation of the helicopter seat shock absorber and achieved the reduction in peak deceleration as one of the control objectives.

Key issues
To achieve adaptive shock mitigation control using MREA, the following aspects can be focused on.For the structural design of MREAs, it is necessary to increase the pressure drop due to the MR yield stress while reduce the viscous pressure drop to improve the dynamic range and the controllable velocity range of MREA.This can be achieved by optimizing dimensions of components including the electromagnetic coil, the MR fluid gap, and the piston.Increasing the total active MR valve length or increasing the external magnetic field without magnetic saturation and reducing the fluid velocity in the annular duct or valve will be beneficial to improving the performance of MREAs.For high-speed impact conditions, typical structural configuration of MREAs used for low-speed vibration conditions may still not meet the performance requirements of shock mitigation control.Therefore, it is necessary to propose new MREA structures suitable for high-speed impact conditions like the internal bypass MREA shown in figure 12.In addition, the performance of the employed MR fluids in MREA for the high-speed impact applications should be considered to achieve high-performance MR control systems.
Due to the high demand for real-time controllability of the shock mitigation control system under high-speed impact conditions, the minimization of the response time of MREA should be considered when designing MREA.Electromagnetic simulation can be used to verify whether the response time of MREA meets the requirements.In the future, quantitative or qualitative research should be conducted on the relationship between the geometrical dimensions and parameters of MREAs and response time to optimize the response time when solving the geometrical dimensions and parameters of MREA.In order to make the actual mechanical performance of the designed MREA closer to the expected mechanical performance, the analytical model should take into account the minor losses and inertia effect of the MR fluid to make the relationship between geometric dimensions and parameters of MREAs and MREA force under high-speed impact conditions closer to the true value.Although there is currently limited research, unsteady models using non-averaged acceleration can more accurately describe the mechanical characteristics of MREA and will become the main research direction in the future.To more realistically describe the characteristics of MREA, the analytical model can consider the influence of the temperature of the MR fluid on the characteristics of the MR fluid, as well as the local losses that occur on passive fluid passage.
The research focus in the field of dynamic modeling is to accurately predict the strong nonlinear characteristics of MREA under high-speed impact conditions.The commonly used adaptive neural fuzzy inference system in the field of dynamic models can predict the nonlinear characteristics well under high-speed impact conditions, but it requires relatively high computational resources.It is difficult to convert into an inverse dynamic model used for the MREA controller.It has become a limiting factor for MR shock mitigation control performance that the MREA controller can accurately describe the strong nonlinear characteristics of MREA and has good real-time performance.
At present, the soft-landing control of SDOF systems is relatively mature.The soft-landing control of shock mitigation control systems such as artillery recoil buffer devices, aircraft landing gears, and landers based on MREA will become one of the main research directions in the field of MR impact buffer control in the future.The soft-landing control of artillery/gun recoil systems, aircraft landing gears, and planetary landers based on MREA will become one of the main research directions in the field of MR shock mitigation control.For actual shock mitigation control systems, there may be specific control performance requirements based on specific application scenarios.For example, the helicopter seat suspensions emphasize controlling the peak deceleration; the aircraft landing gears emphasize improving energy absorption efficiency; the planetary landers emphasize controlling the roll rate to avoid overturning when landing on a slope, and some artillery/gun recoil systems need to quickly reduce the gun body speed to zero to improve the artillery/gun's firing rate.How to achieve soft-landing control under the premise of meeting the control performance requirements of specific applications will become a research focus to further improve control effectiveness in the field of control.For highly complex and nonlinear systems, employing control methods based on artificial intelligence can achieve superior control performance.Artificial intelligence-based control methods, such as neural network models, can address challenges in system modeling, thereby facilitating the application of MREA in largerscale and more complex shock mitigation control systems with stronger nonlinear characteristics.
When evaluating the control performance of the MR shock mitigation control systems, multiple evaluation indicators such as the impact energy absorption, the stability and the efficiency of the shock mitigation control, and the transmission of impact loads to payloads should be selected to comprehensively analyze the characteristics of MR shock mitigation control systems.When optimizing one of the performance evaluation indicators according to the demand of a specific MR shock mitigation control system, it is necessary to avoid the other performance evaluation indicators being too poor.

Conclusions
This paper summarizes the state-of-the-art research in the field of MR shock mitigations, including the structural design and optimization of MREAs, the analytical and dynamic models of MREAs, the control methods of MREAs, and the performance evaluation.To achieve the best performance of shock mitigation controls, MREAs need to be optimized including the electromagnetic coil, the MR fluid gap, and the piston.Then, the passive MREA force can be reduced while increasing the controllable MREA force, and finally, the dynamic range of MREAs will be widened.The MREA with an internal bypass design has very good potential as the actuator of MR shock mitigation control systems because of its excellent properties.The MREA should be specialized based on the classical structure arrangement for different applications to enhance the control performance, and the control objective of the soft-landing can also be considered in the optimal design of MREAs.In addition, the inertial effect and minor losses of the MR fluid should be considered in the analytical model of MREAs at high flow velocities.To accurately describe the strong hysteresis nonlinear characteristics of the MREA force under high-speed impact conditions, it is appropriate to utilize semiempirical models or non-parametric models under high-speed impact conditions with respect to the prediction accuracy and the difficulty of parameter identification.The influence of temperature and other factors on the properties of MR fluids can be also considered to improve the accuracy of the analytical model.For MR shock mitigation control systems put into commercial applications, the characteristics of MREAs should be accurately predicted by the dynamic model under high-speed impact conditions, and the computational demands of the model should be minimized for the MREA controller.For the shock mitigation control system, the most ideal control target is the soft-landing in which the velocity of the control object is exactly reduced to zero when the available stroke is used up.Next, a series of control algorithms to achieve the soft-landing are emphasized in SDOF MR shock mitigation control systems.As the research on the soft-landing control method develops continuously in SDOF MR shock mitigation control systems, the soft-landing control of the artillery/gun recoil system, the aircraft landing gear and the planetary lander based on MREAs will become one of the main directions in the research field of MR shock mitigation control in the future.Finally, the evaluation indexes of performance are also reviewed for the current MR shock mitigation control.Due to the difference in control objectives, different evaluation indexes are selected for shock mitigation control systems applied on various occasions.
a f the instantaneous average acceleration of the fluid flow in the MR fluid gap ap the acceleration of the piston rod Ac the cross-sectional area of the electromagnetic coil gap Ag the cross-sectional area of the active MR fluid gap A h the cross-sectional area of the MR fluid gap Ap the effective piston area Apr the cross-sectional area of the piston rod c the passive damping coefficient D the maximum piston stroke D h the hydraulic diameter of the valve f the Darcy friction factor F d the desired MREA force F f the friction force Fgas the spring force generated by the accumulator Fη the passive MREA force, i.e. the uncontrollable MREA coefficient L the total MR gap length La the total active MR gap length L ni the length of the ith passive fluid gap Lc the length of all electromagnetic coils m the payload mass me the mass of the excitation source rp the radius of the piston rpr the radius of the piston rod v f the fluid velocity in the annular duct or valve v fc the average fluid velocity in the coil gap v i the average fluid velocity at the position where the ith minor losses occur V 0 the initial volume of the accumulator xp the displacement of the piston ẋp the piston velocity xo the payload displacement ẋo the payload velocity ẍo the payload acceleration ẋe the excitation velocity ẋc the velocity of MREA's cylinder ∆Pg the pressure drop produced by the annular damping channel ∆Pη the viscous pressure drop ∆P MR the pressure drop due to MR yield stress ∆P ml the pressure drop due to minor loss ∆P inertia the pressure drop through the MR valve due to inertia effect P 0 the initial pressure of the accumulator ε the average pipe wall roughness υ the kinematic viscosity η the post-yield viscosity ρ the density of the MR fluid τ the shear stress τy the yield stress ζ the damping ratio of the system 1.Introduction 1.1.Magnetorheological (MR) fluids As a kind of smart material, MR fluids consist of micronsized magnetic particles, carrier fluid, and various additives (de Vicente et al 2011).When there is no magnetic field, as shown in figure 1, the magnetic particles in MR fluids are uniformly distributed in the carrier fluid, and MR fluids behave like a Newtonian fluid.Under the application of an external magnetic field, as shown in figure 2, the magnetized particles in the carrier fluid tend to be arranged closely along the direction of magnetic flux paths, causing the apparent viscosity of MR fluids to increase and the fluidity to decline.The apparent viscosity of MR fluids can be controlled by adjusting the external magnetic field strength, and the response time of MR fluids is within 1 ms.The typical constitutive behavior of the MR fluids is shown in figure 3. The yield stress of MR fluids increases with the increasing magnetic field, and its relationship with shear strain initially exhibits an approximate linear relationship.After reaching the critical shear strain, the yield stress remains constant and the MR fluids exhibit non-Newtonian fluid characteristics, with viscosity decreasing with the shear rate (de Vicente et al 2011, Wang and Liao 2011, Pei and Peng 2022, Deng et al 2022a).Thanks to these properties, MR fluids are used in energy absorbers (EAs) (Yazid et al 2014, Elsaady et al 2020a, Abdalaziz et al 2023a), brakes (Nguyen et al 2013, Rossa et al 2013, Shiao and Nguyen 2013, Wellborn et al 2021), clutches (Rizzo 2016, Moghani and Kermani 2019, Fernández and Chang 2021), finishing devices (Kordonski and Gorodkin 2011, Srivastava et al 2021, Ghosh et al 2023, Singh and Jayant 2023), valves (Yoo and Wereley 2004), smart joints (Klingenberg 2001, Gao et al 2017, Khazoom et al 2020, Nordin et al 2022) and other actuators (Ahamed et al 2018, Phu and Choi 2019).In recent years, MR fluids have also been used for smart haptic devices (Park and Choi 2021, Liu et al 2023) and other applications (McDonald et al 2022, Aralikatti and Kumar 2023, Ntella et al 2023).

Figure 1 .
Figure 1.Microstructure of MR fluids in the field-off state: (a) SEM image, (b) schematic diagram.

Figure 2 .
Figure 2. Microstructure of MR fluids in the field-on state: (a) SEM image, (b) schematic diagram.

Figure 3 .
Figure 3.The typical constitutive behavior of MR fluids.

Figure 5 .
Figure 5. Schematic diagram of MREA shock mitigation control system.
, and the proportional-integral-derivative (PID) control (Hu et al 2012, Metered et al 2015), advanced modern control methods such as the adaptive control (Song et al 2005) and the sliding mode control (Dixit and Buckner 2005, Saleh et al 2021), hybrid control methods (Phu et al 2017, Turnip and Panggabean 2020) including the fuzzy PID control (Ahmed et al 2022, Gong et al 2022), the fuzzy sliding mode control (Phu et al 2014) and the fuzzy neural network control (Yu et al 2009) and other control methods (Choi et al 2016, Ma et al 2021) have been commonly used for MR vibration control systems.State-of-the-art control methods have also been systematically reviewed for MR vibration control systems (Masa'id et al 2023).

Figure 6 .
Figure 6.The functional block diagram (left) and full-text architecture diagram (right) of the MR shock mitigation control system to achieve the 'soft-landing' goal.
2.1.Force calculation of MREAsIn recent years, MR actuators including MR brakes (Jang et al 2022, Masjosthusmann et al 2022, Singh and Sarkar 2023), MR clutches (Bira et al 2022), and MREAs (Yan et al 2021, Tan et al 2022, Abdalaziz et al 2023b) have been designed according to the requirements of different applications, and a variety of MREAs dedicated to shock mitigation control have been developed (Jin et al 2021, Wang et al 2021c, Deng et al 2022b).When selecting the MREA for shock mitigation control systems, MREA force calculation is necessary for the structural configuration of MREAs based on the existing structures.

Figure 7 .
Figure 7. Schematic diagram of the single-ended MREA design.

Figure 8 .
Figure 8. Schematic diagram of the double-ended MREA design.

Figure 10 .
Figure 10.Schematic diagram of the MREA design applied to the aircraft landing gear.

Figure 11 .
Figure 11.Schematic diagram of the MREA design applied to the buffer device of the artillery recoil.

Figure 12 .
Figure 12.Schematic diagram of the internal bypass MREA.

Figure 14 .
Figure 14.Schematic diagram of MR fluid flow regions.

Figure 15 .
Figure 15.Schematic diagram of the fluid velocity profile in the parallel plate model.
3.3.DynamicMREA models   3.3.1.Parametric models   3.3.1.1.Bingham model.The Bingham model proposed byStanway et al (1987)  is widely used because of its computational simplicity.Using the friction force of the Coulomb friction element to represent the controllable force of the MREA is characteristic of Bingham-based models.The Bingham model and a modified Bingham model are

3. 3
.1.2.Modified Bouc-Wen model.Bouc-Wen hysteresis operator-based models have been widely used in the MREA modeling.A modified Bouc-Wen model is shown in figure 17, and the expression of this model is given by (Wang and Liao 2011):
3.3.1.4.RC operator-based hysteresis model.Ensuring the accuracy of the experimental data fitting, Bai et al (2019a) proposed a hysteresis model based on the RC operator to reduce the number of model parameters and improve computational efficiency.The model replaced the Bouc-Wen operator with the new RC operator based on the restructured model.This model has the advantages of no differential expression, few parameters, high calculation efficiency, and good storage characteristics.The schematic diagram is shown in figure 19.
3.3.1.5.Parametric models for shock response.In order to describe the dynamic characteristics of MREAs under shock loads, Xiang et al (2008) modified the expression of the hysteresis operator on the basis of a restructured Bouc-Wen model proposed byDominguez et al (2006).They proposed a new fitting expression for the elastic coefficient and the damping coefficient, and obtained a new dynamic model which can be expressed as: where α BW (I) and γ BW (I) are the shaping factors of the forcevelocity hysteresis curve.The fitting expressions of α BW (I) and γ BW (I) are close to the fitting expressions of a restructured Bouc-Wen model proposed by Dominguez et al.

Figure 20 .
Figure 20.Schematic diagram of the dynamic RC operator-based hysteresis model.

Figure 21 .
Figure 21.Topological structure diagram of the TANFIS model.

Figure 22 .
Figure 22.Topological structure diagram of the W-TANFIS model.

Figure 25 .
Figure 25.Schematic diagram of the optimal Bi number controller.

Figure 26 .
Figure 26.Schematic diagram of the hybrid controller for both shock mitigation and vibration attenuation(Bai and Yang, 2019).

Figure 27 .
Figure 27.Dynamic model of the MR landing gear system.

Figure 28 .
Figure 28.Dynamic model of the SDOF shock and vibration control system.

Figure 29 .
Figure 29.Schematic diagram the hybrid controller for MR landing gear.

Figure 30 .
Figure 30.Schematic diagram of the adaptive sliding mode hybrid controller.

Figure 31 .
Figure 31.Schematic diagram of the reference module.

Figure 32 .
Figure 32.Dynamic model of the 2DOF model of the lander.

Figure 33 .
Figure 33.Schematic diagram of the optimal control with fuzzy compensation.
The results of prototype experiments demonstrate that the proposed model can well describe the characteristics of MREAs in the low Reynolds number and high Reynolds number regions.But in the transition region from a low Reynolds number to a high Reynolds number, the simulation results of the model are not smooth, and there are sudden changes.To improve the predictive performance of the model in the transition region, Powers et al (2016) established a new rough wall model.In this model, the Darcy friction factor, f, can be calculated with a piecewise function according to the roughness Reynolds number, Re * : Bira N, Dhagat P and Davidson J R 2022 A low-power magnetorheological fluid clutch utilizing electropermanent magnet arrays Front.Mater.713 1039004 Bitaraf M, Ozbulut O E, Hurlebaus S and Barroso L 2010 Application of semi-active control strategies for seismic protection of buildings with MR dampers Eng.Struct.32 3040-7 Boada M J L, Calvo J A, Boada B L and Díaz V 2011 Modeling of a magnetorheological damper by recursive lazy learning Int.J. Non-Linear Mech.46 479-85 Bui Q D, Bai X X and Nguyen Q H 2021 Dynamic modeling of MR dampers based on quasi-static model and magic formula hysteresis multiplier Eng.Struct.245 112855 Carlson J D 2009 Magnetorheological fluids Smart Materials (CRC Press, Taylor & Francis Group, LLC) Caterino N, Spizzuoco M, Piccolo V and Magliulo G 2022 A semi-active control technique through MR fluid dampers for seismic protection of single-story rc precast buildings Materials 15 759 Chen P, Bai X X, Qian L J and Choi S B 2018 An approach for hysteresis modeling based on shape function and memory mechanism IEEE/ASME Trans.Mechatronics 23 1270-8 Cheng X, Bai Z, Zhu F, Chou C C, Jiang B and Xu S 2022 An optimized bio-inspired thin-walled structure with controllable crashworthiness based on magnetorheological fluid Mech.Adv.Mat.Struct.1-16 Choi S B, Lee S K and Park Y P 2001 A hysteresis model for the field-dependent damping force of a magnetorheological damper J. Sound Vib. 2 375-83 Choi S B, Li W, Yu M, Du H, Fu J and Do P X 2016 State of the art of control schemes for smart systems featuring magneto-rheological materials Smart Mater.Struct.25 043001 Choi Y T and Wereley N M 2005 Biodynamic response mitigation to shock loads using magnetorheological helicopter crew seat suspensions J. Aircr.42 1288-95 Choi Y T and Wereley N M 2015 Drop-induced shock mitigation using adaptive magnetorheological energy absorbers incorporating a time lag J. Vib.Acoust.137 1 Christie M D, Sun S, Deng L, Du H, Zhang S and Li W 2023 Shock absorption for legged locomotion through magnetorheological leg-stiffness control Machines 11 236 de Falco Manuel J G, Bombard A J F and Weeks E R 2023 Effect of polydispersity in concentrated magnetorheological fluids Smart Mater.Struct.32 045014 de Vicente J, Klingenberg D J and Hidalgo-Alvarez R 2011 Magnetorheological fluids: a review Soft Matter 7 3701-10 Deng H, Lian X and Gong X 2022a A brief review of variable stiffness and damping magnetorheological fluid dampers Front.Mater.9 1019426 Deng L, Sun S, Jin S, Li Z, Du H, Zhang S and Li W 2022b Development of a new magnetorheological impact damper with low velocity sensitivity Smart Mater.Struct.31 095042 Ding R, Wang R, Meng X and Chen L 2023 Research on time-delay-dependent H∞/H2 optimal control of magnetorheological semi-active suspension with response delay J. Vib.Control 29 1447-58 Dixit R K and Buckner G D 2005 Sliding mode observation and control for semiactive vehicle suspensions Veh.Syst.Dyn.43 83-105 Dominguez A, Sedaghati R and Stiharu I 2006 A new dynamic hysteresis model for magnetorheological dampers Smart Mater.Struct.15 1179 Dong X, Yu J, Wang W and Zhang Z 2017 Robust design of magneto-rheological (MR) shock absorber considering temperature effects Int.J. Adv.Manuf.Technol.90 1735-47 Ebrahimi B 2009 Development of hybrid electromagnetic dampers for vehicle suspension systems Doctoral Dissertation The University of Waterloo El Majdoub K, Giri F and Chaoui F Z 2020 Adaptive backstepping control design for semi-active suspension of half-vehicle with magnetorheological damper IEEE/CAA J. Autom.Sin. 8 582-96 Elsaady W, Oyadiji S O and Nasser A 2020a Magnetic circuit analysis and fluid flow modeling of an MR damper with enhanced magnetic characteristics IEEE Trans.Magn.56 1-20 Elsaady W, Oyadiji S O and Nasser A 2020b A review on multi-physics numerical modelling in different applications of magnetorheological fluids J. Intell.Mater.Syst.Struct.31 1855-97 Fernández M A and Chang J Y 2021 Design of an adjustable fail-safe MRF clutch with a novel field blocking mechanism for robotic applications IEEE Access 9 164912-27 Fu L, Lin L P and Xu X H 2012 Research on modeling and fuzzy control of magneto-rheological intelligent buffer system for impact load J. Shanghai Jiaotong Univ.17 567-72 Gao F, Liu Y N and Liao W H 2017 Optimal design of a magnetorheological damper used in smart prosthetic knees Smart Mater.Struct.26 035034 Gedik E, Kurt H, Recebli Z and Balan C 2012 Two-dimensional CFD simulation of magnetorheological fluid between two fixed parallel plates applied external magnetic field Comput.Fluids 63 128-34 Ghosh G, Sidpara A and Bandyopadhyay P P 2023 Performance improvement of magnetorheological finishing using chemical etchant and diamond-graphene based magnetic abrasives Precis.Eng.79 221-35 Gołdasz J and Sapiński B 2017 Magnetostatic analysis of a pinch mode magnetorheological valve Acta Mech.Autom.11 229-32 Goncalves F D 2005 Characterizing the behavior of magnetorheological fluids at high velocities and high shear rates Doctoral Dissertation Virginia Tech Gong F, Nie S, Ji H, Hong R, Yin F and Yan X 2022 Pipeline vibration control using magnetorheological damping clamps under fuzzy-PID control algorithm Micromachines 13 531 Gong X, Ruan X, Xuan S, Yan Q and Deng H 2014 Magnetorheological damper working in squeeze mode Adv.Mech.Eng. 6 410158 Gordaninejad F, Saiidi M, Hansen B C, Ericksen E O and Chang F K 2002 Magneto-rheological fluid dampers for control of bridges J. Intell.Mater.Syst.Struct.13 167-80 Grunwald A and Olabi A G 2008 Design of magneto-rheological (MR) valvel Sens. Actuators A 148 211-23 Gu X, Li Y and Li J 2016 Investigations on response time of magnetorheological elastomer isolator for real-time control implementation Smart Mater.Struct.25 11LT04 Han C, Kang B H, Choi S B, Tak J M and Hwang J H 2019 Control of landing efficiency of an aircraft landing gear system with magnetorheological dampers J. Aircr.56 1980-6 Han C, Kim B G and Choi S B 2018 Design of a new magnetorheological damper based on passive oleo-pneumatic landing gear J. Aircr.55 2510-20 Hitchcock G H, Wang X and Gordaninejad F 2007 A new bypass magnetorheological fluid damper J. Vib.Acoust.129 641-7 Hu H, Jiang X, Wang J and Li Y 2012 Design, modeling, and controlling of a large-scale magnetorheological shock absorber under high impact load J. Intell.Mater.Syst.Struct.23 635-45 Huang K and Li J 2022 Research on damping performance and parameter influence of new frame-MRD-core tube composite structural system Structures 41 215-21 Imaduddin F, Mazlan S A and Zamzuri H 2013 A design and modelling review of rotary magnetorheological damper Mater.Des.51 575-91 Ismail M, Ikhouane F and Rodellar J 2009 The hysteresis Bouc-Wen model, a survey Arch.Comput.Methods Eng.16 161-188