A Comprehensive Study on the Optimal Design of Magnetorheological Dampers for Improved Damping Capacity and Dynamical Adjustability

Purpose: We aim to provide a systematic methodology for the optimal design of MRD for improved damping capacity and dynamical adjustability in performing its damping function. Methods: A modified Bingham model is employed to model and simulate the MRD considering the MR fluid’s compressibility. The parameters that describe the structure of MRD and the property of the fluid are systematically examined for their contributions to the damping capacity and dynamically adjustability. A response surface method is employed to optimize the damping force and dynamically adjustable coefficient for a more practical setting related to the parameters. Results: The simulation system effectively shows the hysteretic characteristics of MRDs and shows our common sense understanding that the damping gap width and yoke diameter have significant effects on the damping characteristics of MRD. By taking a typical MRD device setup, optimal design shows an increase of the damping force by 33% and an increase of the dynamically adjustable coefficient by 17%. It is also shown that the methodology is applicable to other types of MDR devices. Conclusion: The compressibility of MR fluid is one of the main reasons for the hysteretic characteristics of MRD. The proposed simulation and optimization methods can effectively improve the MRD’s damping performance in the design stage.


Introduction
Magnetorheological damper (MRD) is a semi-actively controlled intelligent vibration damping device. The rapid, continuous, and reversible rheological properties of the internal magnetorheological fluid make it a smart device [1,2]. Most studies on MRD have been devoted to vehicle vibration and control [3,4]. Recently, MRD has been used in the field of robotics [5]. For example, some work shows that MRD is used in rehabilitation medical devices [6,7]. Compared with the traditional gas damper, MRD can adjust its damping capability (measured by damping force) in real-time to accommodate a change in its environment, and thus improve the robot's controllability and stability as well as its resilience and robustness [8,9].
As a critical component of humanoid robot or interactive robot equipment, MRD has higher requirements on its size and power density than that of typical applications. Therefore, accurate prediction of the MRD's damping force is crucial for optimizing the device's performance. Regarding the optimization design of an MRD, the focus of the parameters to be optimized in the design is different for the different structures of the MRD [10]. Since the MRD is a multi-physics intersection involving electromagnetic field, flow field, structure, etc., to improve the performance of the MRD, the multi-objective optimization design of multi-physics coupling is a commonly used optimization design method [11,12]. Parlak et al. [13] performed a unidirectional coupling simulation analysis of MRD to study the effects of the current, damping gap, piston core radius, and flange length on the MRD's damping characteristics. Gurubasavaraju et al. [14] combined the finite element method with computational fluid dynamics to calculate the damper's maximum output force with a specified frequency, amplitude, and applied current. The effect of damping gap length and frequency on the damping force of MRD has also been investigated. Case et al. [15] performed multiphysics modeling on the dynamics of the small MRD and calculated the magnetic flux density, magnetic induction intensity, and the circulating fluid flow state in the MRD under sinusoidal excitation of the piston. The simulation result was consistent with the maximum damping force obtained from the experiment. Besides, in terms of power consumption and overall weight, aiming at reducing computational complexity, a comprehensive approach based on mathematical optimization and finite element analysis for the optimal design of MRDs was conducted [16]. The aforementioned studies could not accurately predict the MRD's hysteretic characteristics at low velocity, which is important to understanding the device's damping capability along with reason.
Indeed, maximization of the damping force and dynamically adjustable coefficient (which is a measure of the adjustability of MRD) are two important concerns from the applications point of view. Khan et al. [17] proposed an optimization method using a comparative analysis. The MRD's damping forces with nine different coils and different piston heads were simulated and analyzed. They found that the piston head can provide a higher damping force when the piston head is adjusted in the shape of rounded corners. These parameters are also factors that researchers often consider in the optimization design of an MRD [18]. Parlak et al. [19] proposed a DOE (design of experiments) optimization design method with the damping force and dynamically adjustable coefficient as two objectives, and the method was conducted by combining with the FEM method to achieve a multi-objective optimization goal [20]. Sung et al. [21] proposed a first-order optimization method using a general-purpose finite element package (ANSYS APDL). In addition, a multiphysics simulation method was performed using COMSOL Multiphysics to enhance an MRD's damping force [22]. The maximum output damping force was chosen as the optimization objective to optimize the MRD design, while the initial and optimized dampers were experimentally evaluated. The optimized one can provide superior performance in terms of the output damping force and power consumption.
An important shortcoming with the aforementioned studies on optimization of MRD design is the lack of comprehensive analysis of various parameters of the structure and property of MRD on how they contribute to the damping capability and dynamically adjustable range. This shortcoming further leads to ad-hoc choices of parameters for optimization, which either fail to get a truly optimized design or lead to a design with many side effects (although good damping force). The study presented in this paper attempted to overcome these shortcomings. The methodology taken in this study had three steps. In the first step, a simulation system for MRD was built. In the second step, the influence of the parameters of the structure and the material property was examined with the help of the simulation system, which led to a set of important parameters for optimization. In the third step, a response surface method was employed to optimize these parameters for two objectives: the damping force and the dynamically adjustable coefficient. The result of the optimization was validated with the simulation system. Throughout, an example MRD was taken; however, the generality of the methodology to other MRD devices was kept whenever possible.

Fundamental Principles of MRD
The studied MRD's compact volume makes it suitable for rehabilitation medical equipment or more generally human-computer interactive robots. A typical MRD device Actuators 2021, 10, 64 3 of 21 is shown in Figure 1, consisting of the piston (with a double coil inside), piston rod, cylinder, upper and lower end caps, floating piston, lifting lug, and other components. The high-pressure nitrogen compensation chamber is utilized to adjust the floating piston to compensate for the volume change caused by the piston rod entering and exiting.

Fundamental Principles of MRD
The studied MRD's compact volume makes it suitable for rehabilitation medical equipment or more generally human-computer interactive robots. A typical MRD device is shown in Figure 1, consisting of the piston (with a double coil inside), piston rod, cylinder, upper and lower end caps, floating piston, lifting lug, and other components. The high-pressure nitrogen compensation chamber is utilized to adjust the floating piston to compensate for the volume change caused by the piston rod entering and exiting. The magnetic flux generated by the excitation coil passes through the damping gap (the annular gap between the piston and the cylinder) in the middle and both ends of the piston. The magnetic flux passes through the cylinder via the damper gap while it passes through the other side's damper gap. Eventually, the magnetic flux flows back to the yoke to form a closed magnetic circuit. The applied magnetic field changes the viscosity of MR fluid in the damping gap as well as the damping force for the piston to move. Meanwhile, in the process of piston movement, as the current increases, the MR fluid viscosity increases, and the MR fluid's yield stress increases accordingly. Therefore, the MRD's damping force can be adjusted by the current applied to the excitation coil.

Mathematical Model for the Damping Force
According to the theory for the gap flow between two parallel plates, the MRD's damping force can be expressed as [23,24]: The magnetic flux generated by the excitation coil passes through the damping gap (the annular gap between the piston and the cylinder) in the middle and both ends of the piston. The magnetic flux passes through the cylinder via the damper gap while it passes through the other side's damper gap. Eventually, the magnetic flux flows back to the yoke to form a closed magnetic circuit. The applied magnetic field changes the viscosity of MR fluid in the damping gap as well as the damping force for the piston to move. Meanwhile, in the process of piston movement, as the current increases, the MR fluid viscosity increases, and the MR fluid's yield stress increases accordingly. Therefore, the MRD's damping force can be adjusted by the current applied to the excitation coil.

Mathematical Model for the Damping Force
According to the theory for the gap flow between two parallel plates, the MRD's damping force can be expressed as [23,24]: where where F η , F τ , F g , and F f are the viscous damping force, the Coulomb damping force, the gas compensating force, and the friction force, respectively; η 0 is the MR fluid's viscosity in the absence of magnetic field; L denotes the piston length; L 1 and L 2 are the effective pole lengths in the middle and two ends, respectively; D and d describe the diameters of the piston, and the piston rod, respectively; D is the mean circumference of the damper's annular flow path; V 0 denotes the gas initial volume in the compensation chamber; x is the Actuators 2021, 10, 64 4 of 21 piston displacement; and v is the piston velocity. If the friction force is not considered, the dynamically adjustable coefficient λ of the MRD can be written as: 3. Simulation of MRD 3.1. The Fundamental Theory of MR Damping Force The governing equations of electromagnetism in the COMSOL solver, including Maxwell's equation and Ampere's law, are given by [15]: where H, J, D, and B are the magnetic field, current densities, electric displacement field, and magnetic flux densities, respectively. The flow of the magnetorheological fluid in the MRD can be described by the following weakly compressible Navier-Stokes equation to describe the velocity and pressure fields in the space motion coordinate system [25]: where ρ is the fluid density; u is the velocity vector; p is the pressure; I is the unit diagonal matrix; ∇ is the vector operator in Cartesian coordinates; η is the kinetic viscosity; and F is the volume force. The correlation between flow and electromagnetic fields in the MRD is described by the viscosity of MR fluid and magnetic induction intensity. Therefore, this study's solution strategy is to first obtain the magnetic induction intensity distribution by solving the electromagnetic field and then load the solved magnetic induction intensity distribution result as the initial values into the flow field problem for solution. To solve the nondifferentiable constitutive equation of the magnetorheological fluid without yield, the modified Bingham model [26] is utilized to model the magnetorheological fluid viscosity as: where η 0 is the zero-field viscosity; τ y (B) is the yield stress associated with the magnetic field; . γ is the shear rate; and m is a model parameter determining the viscosity growth of non-Newtonian fluids at low shear rates [13,27]. According to the experimental data in our research group, the shear stress and magnetic induction intensity of MRF450 was obtained as: For the flow field analysis of the MRD in the moving process, the relationship between the density and pressure of the MR fluid under the MR fluid compressibility is described using the Tait equation of state [28] as: where M v is the bulk elastic modulus, taking 3 × 10 8 Pa; ρ is the compressible fluid density; ρ 0 is the reference density; p 0 is the reference pressure; and n is the index regulating the fluid stiffness.

Material Properties and Meshing
The MRD's piston adopts DT4, which is a pure electrical iron with low coercivity, excellent processing performance, and a high magnetic flux performance; the piston rod and the cylinder adopt steel 45 and 20, respectively. The magnetization curves of each part of materials can be obtained from [29], while the magnetization curve of MRF450 is shown in Figure 2. Table 1 presents the MR fluid properties, and the main design parameters of the MRD are shown in Table 2.
tween the density and pressure of the MR fluid under the MR fluid compressibility is described using the Tait equation of state [28] as: where v M is the bulk elastic modulus, taking 8

10 Pa ×
; ρ is the compressible fluid density; 0 ρ is the reference density; 0 p is the reference pressure; and n is the index regulating the fluid stiffness.

Material Properties and Meshing
The MRD's piston adopts DT4, which is a pure electrical iron with low coercivity, excellent processing performance, and a high magnetic flux performance; the piston rod and the cylinder adopt steel 45 and 20, respectively. The magnetization curves of each part of materials can be obtained from [29], while the magnetization curve of MRF450 is shown in Figure 2. Table 1 presents the MR fluid properties, and the main design parameters of the MRD are shown in Table 2.

Model Number
Properties Value  The model should be reasonably simplified to reduce the calculation load. Therefore, the upper and lower end caps, seals, and fasteners are ignored in the modeling procedure. A two-dimensional axisymmetric model is established, and the properties of the materials above are incorporated into the corresponding modules of the solver. The meshing Actuators 2021, 10, 64 6 of 21 sequence controlled by the physical field, including flow model settings, specific feature processing, and grid cell size control, is utilized to divide the mesh (see Figure 3), and the mesh is refined at the damping gap and outlet of the damper. Boundary layer treatment is performed at the fluid boundary. The complete grid consists of 10,459 grid cells and 1699 edge cells.
number of turns per coil N 600 compensation chamber pressure a P

MPa
The model should be reasonably simplified to reduce the calculation load. Therefore, the upper and lower end caps, seals, and fasteners are ignored in the modeling procedure. A two-dimensional axisymmetric model is established, and the properties of the materials above are incorporated into the corresponding modules of the solver. The meshing sequence controlled by the physical field, including flow model settings, specific feature processing, and grid cell size control, is utilized to divide the mesh (see Figure 3), and the mesh is refined at the damping gap and outlet of the damper. Boundary layer treatment is performed at the fluid boundary. The complete grid consists of 10,459 grid cells and 1699 edge cells.

Boundary Conditions and Initial Values
Boundary conditions are necessary to find electromagnetic and flow fields and should have permissible values. Here, the magnetic field's air boundary condition is chosen as the magnetic insulation, while the boundary condition of the flow field in contact with the piston and cylinder is adjusted as a non-slip flow on the wall. To more accurately simulate the transient flow field variations of the magnetorheological fluid inside the damper under the magnetic field action, the dynamic grid technology is utilized to describe the piston displacement through the following sinusoidal relation: where A is the piston displacement amplitude, while its initial value is adjusted as 10 mm; f is the piston's vibration frequency, while its initial value is chosen as 1 Hz; and t s is the time step, which is considered as 0.01 s.

Simulation Results
The magnetic induction intensity and magnetic induction line distribution of the MRD at different currents (for I = 0.25, 0.5, 0.75, and 1 A) are shown in Figure 4a-d. The magnetic induction intensity at the midline of the damping gap under different currents is extracted (see Figure 5) to further understand the distribution of magnetic induction intensity in the MR region. It can be seen that the magnetic induction intensity gradually increases with the increase of current. In the damping gap, the magnetic induction intensity in the effective damping gap between the two coils reaches its maximum while rapidly reduces in the damping gap area adjacent to the coils. magnetic induction intensity at the midline of the damping gap under different currents is extracted (see Figure 5) to further understand the distribution of magnetic induction intensity in the MR region. It can be seen that the magnetic induction intensity gradually increases with the increase of current. In the damping gap, the magnetic induction intensity in the effective damping gap between the two coils reaches its maximum while rapidly reduces in the damping gap area adjacent to the coils.   The pressure at the middle line of the damping gap with different currents is shown in Figure 6. There is a significant pressure rise in the effective damping gap, causing the damping force.    The pressure at the middle line of the damping gap with different currents is shown in Figure 6. There is a significant pressure rise in the effective damping gap, causing the damping force. Figure 7a-d shows the pressure distribution of MR fluid at different currents (for I = 0.25, 0.5, 0.75, and 1 A) when the piston rebounds upward. When I = 1 A, the damping pressure rise at both ends of the piston reaches 1.25 MPa. The pressure at the middle line of the damping gap with different currents is shown in Figure 6. There is a significant pressure rise in the effective damping gap, causing the damping force. Figure Figure 9a-d, which indicates that the MR fluid viscosity adjacent to the effective magnetic poles gradually increases with the increase of current, while the MR fluid viscosity adjacent to the effective magnetic poles in the middle is higher than that at both ends. The MR fluid's maximum viscosity in the damping gap reaches 127 Pa · s( I = 1 A).
Two damping force-displacement and force-velocity curves are utilized to describe the MRD's damping characteristics under sinusoidal excitation. The damping forcedisplacement curve describes the relationship between the MRD's output damping force and piston displacement, reflecting the MRD's ability to attenuate vibration in one motion period, while the damping force-velocity curve describes the relationship between the damper's output damping force and the piston velocity, reflecting the damping force response to the piston velocity. The damping characteristics of dampers considering compressibility and incompressibility are shown in Figure 10 for I = 1 A, A = 10 mm, and f = 1 Hz. According to Figure 10a, the simulation model considering the MR fluid incompressibility and the simulation model considering the MR fluid compressibility have good agreement on the MRD's maximum output damping force. However, near the maximum displacement, when the piston reaches its maximum displacement, the damping force does not change immediately but gradually. It also should be noticed that the damping force is not symmetric about the x-axis due to the role of the gas compensation chamber. When the piston is in extrusion motion, more pressurized gas is applied, which applies a greater force than in rebound motion. Actuators 2021, 10, x FOR PEER REVIEW 10 of 22   According to Figure 10b, compared with the incompressible model, it demonstrates a strong nonlinear relationship between the output damping force predicted by the simulation model (considering the magnetorheological fluid compressibility) and the velocity. Moreover, an evident hysteretic phenomenon can be observed in the low-velocity region, which cannot be predicted by neglecting the MR fluid compressibility, as one of the essential factors for the MRD's hysteretic characteristics at low velocity in previous studies. Actuators 2021, 10, x FOR PEER REVIEW 11 of 22    Two damping force-displacement and force-velocity curves are utilized to describe the MRD's damping characteristics under sinusoidal excitation. The damping force-displacement curve describes the relationship between the MRD's output damping force and piston displacement, reflecting the MRD's ability to attenuate vibration in one motion period, while the damping force-velocity curve describes the relationship between the damper's output damping force and the piston velocity, reflecting the damping force response to the piston velocity. The damping characteristics of dampers considering compressibility and incompressibility are shown in Figure 10 for I = 1 A, A = 10 mm, and f = 1 Hz. According to Figure 10a, the simulation model considering the MR fluid incompressibility and the simulation model considering the MR fluid compressibility have good agreement on the MRD's maximum output damping force. However, near the maximum displacement, when the piston reaches its maximum displacement, the damping force does not change immediately but gradually. It also should be noticed that the damping force is not symmetric about the x-axis due to the role of the gas compensation chamber. When the piston is in extrusion motion, more pressurized gas is applied, which applies a greater force than in rebound motion. According to Figure 10b, compared with the incompressible model, it demonstrates a strong nonlinear relationship between the output damping force predicted by the simulation model (considering the magnetorheological fluid compressibility) and the velocity. Moreover, an evident hysteretic phenomenon can be observed in the low-velocity region,  Two damping force-displacement and force-velocity curves are utilized to describe the MRD's damping characteristics under sinusoidal excitation. The damping force-displacement curve describes the relationship between the MRD's output damping force and piston displacement, reflecting the MRD's ability to attenuate vibration in one motion period, while the damping force-velocity curve describes the relationship between the damper's output damping force and the piston velocity, reflecting the damping force response to the piston velocity. The damping characteristics of dampers considering compressibility and incompressibility are shown in Figure 10 for I = 1 A, A = 10 mm, and f = 1 Hz. According to Figure 10a, the simulation model considering the MR fluid incompressibility and the simulation model considering the MR fluid compressibility have good agreement on the MRD's maximum output damping force. However, near the maximum displacement, when the piston reaches its maximum displacement, the damping force does not change immediately but gradually. It also should be noticed that the damping force is not symmetric about the x-axis due to the role of the gas compensation chamber. When the piston is in extrusion motion, more pressurized gas is applied, which applies a greater force than in rebound motion. According to Figure 10b, compared with the incompressible model, it demonstrates a strong nonlinear relationship between the output damping force predicted by the simulation model (considering the magnetorheological fluid compressibility) and the velocity. Moreover, an evident hysteretic phenomenon can be observed in the low-velocity region,

Structural Parameters Influence
In this section, the established finite element simulation model is employed to study the influence of several structural parameters, including the damping gap, effective magnetic poles at both ends, effective magnetic poles in the middle, and yoke diameter.

Damping Gap
The damping gaps are chosen as g = 0.8, 1, and 1.2 mm under different currents, while the other parameters remain unchanged. The damping characteristics predicted by the simulation model are shown in Figure 11. The output damping force decreases with the increase of the damping gap width. The area enclosed by the damping force-displacement in Figure 11a decreases with the increase of the damping gap width, indicating that the MRD's energy dissipation performance decreases. It is not difficult to be explained that, with the increase of the damping gap width, the magnetic induction intensity in the damping gap, the yield stress of the MR fluid, and the Coulomb damping force decrease accordingly. Figure 11b shows that the slope of the damping force and velocity curve decreases, while the hysteretic width of the velocity curve reduces. That may be due to the increase of damping gap, which alleviates the pressure inside the damper and changes the MR fluid compressibility accordingly. In summary, although the damping gap width variation is small, it has a significant impact on the MRD's output performance. When the damping gap is relatively large, the damping force output is slightly affected by the amplitude and velocity variations, demonstrating the superior controllability of the MRD. damping gap, the yield stress of the MR fluid, and the Coulomb damping force decrease accordingly. Figure 11b shows that the slope of the damping force and velocity curve decreases, while the hysteretic width of the velocity curve reduces. That may be due to the increase of damping gap, which alleviates the pressure inside the damper and changes the MR fluid compressibility accordingly. In summary, although the damping gap width variation is small, it has a significant impact on the MRD's output performance. When the damping gap is relatively large, the damping force output is slightly affected by the amplitude and velocity variations, demonstrating the superior controllability of the MRD.

Magnetic Pole Length at Both Ends
The MRD's damping characteristics are investigated when effective magnetic poles at both ends are considered as L 1 = 3, 5, and 9 mm under different currents, while the other parameters are constant. The damping characteristics predicted by the simulation model are shown in Figure 12. As the effective pole length increases, the damping force decreases. That is because, when the effective pole length increases, the magnetic flux density via the damping gap declines, and the shear yield stress of the MR fluid and the damper's Coulomb damping force descend accordingly. However, the viscous damping force increases with the increase of the total damping gap length. Due to these two factors' combined action, the damper's output damping force increases first and decreases in a small range. The area enclosed by the damping force-displacement in Figure 12a increases first and then drops, while the energy dissipation performance is approximately constant. Since the output damping force and the MR fluid compressibility change slowly, as shown in Figure 12b, the damping force has the same trend with the slope and hysteresis width of the velocity curve. It can be drawn from above that the effective magnetic pole lengths at both ends have little influence on the MRD's output performance.
first and then drops, while the energy dissipation performance is approximately constant. Since the output damping force and the MR fluid compressibility change slowly, as shown in Figure 12b, the damping force has the same trend with the slope and hysteresis width of the velocity curve. It can be drawn from above that the effective magnetic pole lengths at both ends have little influence on the MRD's output performance.

Magnetic Pole Length in the Middle
The MRD's damping characteristics are evaluated when effective magnetic pole length in the middle is selected as 2 L = 5, 8, 10, and 15 mm under various currents, while the remaining parameters are kept fixed. The damping characteristics estimated by the simulation model are presented in Figure 13. The effect of effective magnetic poles in the middle on MRD characteristics is similar to that of effective magnetic poles at both ends, Velocity (mm/s) Velocity (mm/s) Velocity (mm/s) Velocity (mm/s)

Magnetic Pole Length in the Middle
The MRD's damping characteristics are evaluated when effective magnetic pole length in the middle is selected as L 2 = 5, 8, 10, and 15 mm under various currents, while the remaining parameters are kept fixed. The damping characteristics estimated by the simulation model are presented in Figure 13. The effect of effective magnetic poles in the middle on MRD characteristics is similar to that of effective magnetic poles at both ends, which is not described. Compared with the two effective magnetic poles, the middle's effective magnetic pole has a more significant effect on the MRD's output performance. However, the damping force output performance is affected by the amplitude and velocity variations, and the trend is nearly unchanged.
which is not described. Compared with the two effective magnetic poles, the middle's effective magnetic pole has a more significant effect on the MRD's output performance. However, the damping force output performance is affected by the amplitude and velocity variations, and the trend is nearly unchanged.

Yoke Diameter
The MRD's damping characteristics are verified when the yoke diameters are 1 d = 12, 14, and 16 mm under several currents while considering constant values for the other parameters. The damping characteristics of various currents predicted by the simulation model are depicted in Figure 14. The yoke diameter is different from the above parameters since it does not change the damping gap structure. Accordingly, it does not affect the viscous damping force. The increase in the yoke diameter increases the magnetic induction in the damping gap, accordingly the MR fluid's yield stress and the damper's output damping force. The area enclosed by the damping force-displacement in Figure 14a increases with the increase of yoke diameter, indicating that the MRD's energy dissipation performance increases. Figure 14b illustrates that, when the damping force and velocity curve slope increase, the velocity curve's hysteretic width increases for the fluid's compressibility improvement resulting from improved damping force output performance. The analysis states that, when the damping gap structure is constant, the yoke diameter change significantly influences the MRD's damping force output. However, the damping force output performance is influenced by amplitude, and velocity does not change much.

Yoke Diameter
The MRD's damping characteristics are verified when the yoke diameters are d 1 = 12, 14, and 16 mm under several currents while considering constant values for the other parameters. The damping characteristics of various currents predicted by the simulation model are depicted in Figure 14. The yoke diameter is different from the above parameters since it does not change the damping gap structure. Accordingly, it does not affect the viscous damping force. The increase in the yoke diameter increases the magnetic induction in the damping gap, accordingly the MR fluid's yield stress and the damper's output damping force. The area enclosed by the damping force-displacement in Figure 14a increases with the increase of yoke diameter, indicating that the MRD's energy dissipation performance increases. Figure 14b illustrates that, when the damping force and velocity curve slope increase, the velocity curve's hysteretic width increases for the fluid's compressibility improvement resulting from improved damping force output performance. The analysis states that, when the damping gap structure is constant, the yoke diameter change significantly influences the MRD's damping force output. However, the damping force output performance is influenced by amplitude, and velocity does not change much.

Optimization Based on RSM
As an optimal design method based on statistics and engineering [30], RSM (Response surface methodology) is an alternative model which obtains relevant data through a series of experimental schemes determined in the design domain and employs them to extract the explicit relationship between design variables and response values. When a response or a set of responses is affected by multiple variables, the level of these design variables can be determined to obtain an optimal response. Here, the damper structure is optimized based on response surface methodology (RSM) to achieve higher damping force and damping performance of the adjustable damping coefficient.

Response Surface Test Design
The output damping and dynamically adjustable coefficient are the two most important parameters when evaluating the MRD's overall performance [31]. Therefore, the maximum damping force and dynamically adjustable coefficient are chosen as the optimization objectives. According to the optimization goals, the optimization design process is designed (see Figure 15).

Optimization Based on RSM
As an optimal design method based on statistics and engineering [30], RSM (Response surface methodology) is an alternative model which obtains relevant data through a series of experimental schemes determined in the design domain and employs them to extract the explicit relationship between design variables and response values. When a response or a set of responses is affected by multiple variables, the level of these design variables can be determined to obtain an optimal response. Here, the damper structure is optimized based on response surface methodology (RSM) to achieve higher damping force and damping performance of the adjustable damping coefficient.

Response Surface Test Design
The output damping and dynamically adjustable coefficient are the two most important parameters when evaluating the MRD's overall performance [31]. Therefore, the maximum damping force and dynamically adjustable coefficient are chosen as the optimization objectives. According to the optimization goals, the optimization design process is designed (see Figure 15). The main structural parameters studied in the previous section are considered as the optimization variables, while their upper and lower limits are presented in Table 3. By considering the Box-Behnken scheme for experimental design, 25 groups of schemes are generated, and the corresponding response values for each group are shown in Table 4.   The main structural parameters studied in the previous section are considered as the optimization variables, while their upper and lower limits are presented in Table 3. By considering the Box-Behnken scheme for experimental design, 25 groups of schemes are generated, and the corresponding response values for each group are shown in Table 4.

Response Surface Model
Regression analysis is commonly utilized in RSM to represent polynomial functions. The RSM-based approximate response can be described as: where x 1 ,x 2 , . . . x n are design variables and ε is the error between the actual response value and the approximate one. In this paper, the RSM-based approximate function fitting is directly described by the following quadratic model: From left to right, on the right side of the equation are the intercept, linear, quadratic interaction, and square terms of the model. Regression analysis is applied to the design variables and their response values through Design Expert. Table 5 gives the finallyestablished response functions.
The relationship between the other two design variables obtained by fixing two design variables, the damping force and adjustable damping coefficient, is shown in Figures 16 and 17. It indicates the change and distribution of the damping force and damping adjustable coefficient while changing the design variables. Meanwhile, the relationship between the two response values and different design variables can be obtained through the control variable method, which is not listed here.  The variance analysis of damping force d F and adjustable damping coefficient λ obtained from the RSM is shown in Table 6. The variance analysis can evaluate the response surface function's applicability through the two variance indices F, P, and the determination coefficient R 2 . According to the authors of [32], as both responses' p values are lower than 0.001, it indicates that design variables can significantly affect the response values. R 2 of both responses is higher than the standard value of 0.9, demonstrating the prediction model validity for both responses.

Optimization Results
The MRD's maximum output damping force determines the robot's impact on energy absorption and stability during walking, while the dynamically adjustable coefficient determines the adjustment range of the damping force. The following multi-objective optimization model [33,34] containing two responses is established to increase the output damping force and dynamically adjustable coefficient of the MRD:  The variance analysis of damping force d F and adjustable damping coefficient λ obtained from the RSM is shown in Table 6. The variance analysis can evaluate the response surface function's applicability through the two variance indices F, P, and the determination coefficient R 2 . According to the authors of [32], as both responses' p values are lower than 0.001, it indicates that design variables can significantly affect the response values. R 2 of both responses is higher than the standard value of 0.9, demonstrating the prediction model validity for both responses.

Optimization Results
The MRD's maximum output damping force determines the robot's impact on energy absorption and stability during walking, while the dynamically adjustable coefficient determines the adjustment range of the damping force. The following multi-objective optimization model [33,34] containing two responses is established to increase the output damping force and dynamically adjustable coefficient of the MRD: The variance analysis of damping force F d and adjustable damping coefficient λ obtained from the RSM is shown in Table 6. The variance analysis can evaluate the response surface function's applicability through the two variance indices F, p, and the determination coefficient R 2 . According to the authors of [32], as both responses' p values are lower than 0.001, it indicates that design variables can significantly affect the response values. R 2 of both responses is higher than the standard value of 0.9, demonstrating the prediction model validity for both responses.

Optimization Results
The MRD's maximum output damping force determines the robot's impact on energy absorption and stability during walking, while the dynamically adjustable coefficient determines the adjustment range of the damping force. The following multi-objective optimization model [33,34] containing two responses is established to increase the output damping force and dynamically adjustable coefficient of the MRD: where D is an expectation function, reflecting the appropriate scope for each response.
The expectation function combines several responses into a dimensionless performance measure, while the expectation function can be changed for each response to specify different importance and weights. This study employs the expectation method to solve multi-goals optimization problems.
The structural parameters, maximum damping force, and dynamically adjustable coefficient of the initial and optimization models of the MRD are shown in Table 7. Compared with the initial model, the optimized model's damping force is increased by 33% and the dynamically adjustable coefficient is increased by 17%. The optimized model provides superior performance and higher values for the damping force and dynamically adjustable coefficient. The structural parameters of the optimization model are substituted into the simulation model to solve the problem. The damping characteristics before and after optimization are shown in Figure 18. The indicator area enclosed by the optimal damping force and displacement curve is increased, demonstrating the optimal damper's superior energy dissipation performance and more robustness of the damping force's output performance under the amplitude and velocity variations. Simultaneously, the damping force and velocity curve slope increase, which means that the piston's damping force significantly changes and responds quickly at the stroke's two ends. The optimal MRD provides superior stability and controllability for the structural system.
where D is an expectation function, reflecting the appropriate scope for each response.
The expectation function combines several responses into a dimensionless performance measure, while the expectation function can be changed for each response to specify different importance and weights. This study employs the expectation method to solve multi-goals optimization problems.
The structural parameters, maximum damping force, and dynamically adjustable coefficient of the initial and optimization models of the MRD are shown in Table 7. Compared with the initial model, the optimized model's damping force is increased by 33% and the dynamically adjustable coefficient is increased by 17%. The optimized model provides superior performance and higher values for the damping force and dynamically adjustable coefficient.  The structural parameters of the optimization model are substituted into the simulation model to solve the problem. The damping characteristics before and after optimization are shown in Figure 18. The indicator area enclosed by the optimal damping force and displacement curve is increased, demonstrating the optimal damper's superior energy dissipation performance and more robustness of the damping force's output performance under the amplitude and velocity variations. Simultaneously, the damping force and velocity curve slope increase, which means that the piston's damping force significantly changes and responds quickly at the stroke's two ends. The optimal MRD provides superior stability and controllability for the structural system.

Conclusions
In this study, the modified Bingham constitutive relation for MR fluid was utilized to establish a simulation model considering fluid compressibility. The MRD's damping

Conclusions
In this study, the modified Bingham constitutive relation for MR fluid was utilized to establish a simulation model considering fluid compressibility. The MRD's damping characteristics under sinusoidal excitation were predicted, and a feasible and effective optimization design method was proposed. The damping force and dynamically adjustable coefficient of the MRD are significantly improved. The main conclusions can be drawn as follows: (1) Hysteresis phenomenon for the MRD's damping characteristics is exhibited from the simulation system, which shows the simulation's validity. The MR fluid compressibility is one of the essential factors affecting the MRD's low-velocity hysteretic characteristics. As such, the enabler to simulate and predict the hysteresis characteristics of MRD is an important contribution to the optimal design of MDR devices.
(2) Compared with the effective pole length, the damping gap width and yoke diameter significantly influence the MRD's damping characteristics. When the damping gap is relatively large, the damping force is slightly affected by the amplitude and velocity variation, while the MRD controllability is improved. The yoke diameter has a significant effect on the MRD's damping force output, while the changing trend of damping force output performance affected by amplitude and speed does not change much.
(3) The maximal damping force output performance and dynamically adjustable coefficient of an MRD device can be increased by 33% and 17%, respectively, with the optimization procedure developed in this study.
A physical test-bed will be built to verify the simulation system further and improve the optimal design quality in the near future. Another future work is a more careful examination of the damping characteristics considering the friction behavior in the damping force equation. The friction behavior may not be governed by Coulomb's law, as conventionally thought, but by other laws such as the LuGre model (see a comprehensive discussion of friction in micro-motion systems in [35]).