# Design and Multiobjective Optimization of Magnetorheological Damper considering the Consistency of Magnetic Flux Density

Authors: Zhizhen Dong; Zhimin Feng; Yuehua Chen; Kefan Yu; Gang Zhang
Journal: Shock and Vibration
 (2020)

---

## Abstract

The consistency of magnetic flux density of damping gap (CMDG) represents the balancing magnetic flux density in each damping gap of magnetorheological (MR) dampers. It can make influences on the performances of MR dampers and the accuracy of relevant objective functions. In order to improve the mechanical performances of the MR damper with a two-stage coil, the function for calculating CMDG needs to be found. By establishing an equivalent magnetic circuit model of the MR damper, the CMDG function is derived. Then, the multiobjective optimization function and the working flow of optimal design are presented by combining the parallel-plate model of the MR damper with the function posed before. Taking the damping force, the dynamic range, the response time, and the CMDG as the optimization objective, and the external geometric dimensions of the SG-MRD60 damper as the bound variable, this paper optimizes the internal geometric dimensions of MR damper by using a NSGA-III algorithm on the PlatEMO platform. The results show that the obtained scheme in Pareto-optimal solutions has existed with better performance than that of SG-MRD60 scheme. According to the results of the finite element analysis, the multiobjective optimization design including the CMDG function can improve the uniformity of magnetic flux density of the MR damper in damping gap, which meets the requirements of manufacture and application.


---

## Body

## 1. Introduction

As a semiactive controller for vibration attenuation, the magnetorheological (MR) dampers with excellent electromagnetic controllability have been widely used in various engineering applications [1–4]. A MR damper with high performance not only needs excellent MR fluid [5] but also should have a great damping force and a large dynamic range under certain constraints of geometric dimensions [6]. In order to improve the performance of MR dampers, many methods have been proposed to improve the internal structure of MR dampers, such as multistage piston [7], annular radial channel design [8], variable resistance gap [9], parallel double coil arrangement [10], and meandering magnetic circuit design [11]. However, the new structures resulted from these methods are usually complicated MR damper structures with expensive manufacturing costs and many maintenance difficulties, which significantly limit their engineering applications.Considering the problems arised from the change of the internal structure, many scholars have turned to the multiobjective optimization design method based on the original structure of MR dampers [12]. As we can see, they basically used the objective function optimization method and the finite element optimization method. Guan et al. [13] and Parlak et al. [14, 15] used the damping force and the dynamic range as the optimization objectives to optimize the performances of the MR damper, while Parlak et al. [16] used the damping force and the magnetic flux density as the optimization objectives based on the multiobjective genetic algorithm, Taguchi experimental method, finite element analysis, and CFD analysis technique. Compared to these studies that involve only two-objective optimization for MR dampers, Nguyen and Choi [17] used the specific volume of a vehicle MR damper as a geometric dimension limiting factor and transformed the three-objective optimization problem which involves the damping force, the dynamic range, and the response time into a single-objective optimization problem by means of a comprehensive evaluation function. By combining the results of the Bingham model with finite element analysis, they have identified the geometric dimension of the damper. As for the finite element optimization method, Ferdaus et al. [18] and Hu et al. [19] have explored the impact of different piston shapes on the performance of MR dampers by means of the finite element analysis method. It suggests optimization of the performances of MR dampers without changing its original structure on the basis of changing the internal geometric dimensions of that. Taken the temperature effects into consideration, besides, Dong et al. [20] used the Six Sigma robust method to optimize the main geometrical parameters of MR dampers. All the aforementioned optimization design methods for MR dampers only take two even less performance indexes into account at the same time. In addition, they usually take the middle piston thickness as twice of the flank piston thickness and do not specifically deduce and use any function to optimize the thickness of flank piston, which can influence the performance of MR dampers to some extent. Therefore, the optimization results of precise performance of MR damper are difficult to obtain.This paper takes the damping force, the dynamic range, the response time, and the function of the consistency of magnetic flux density of damping gap (CMDG) as the optimization objectives into account simultaneously and proposes a multiobjective optimization design method for MR dampers based on NSGA-III algorithm [21]. The SG-MRD60 damper applied in sea-crossing bridges is used as a subject, and the piston diameter, rod diameter, coil slot length, coil groove depth, cylinder body thickness, and flank piston thickness are used as variables of the geometric dimension. The formula for calculating CMDG performance of MR dampers is derived by establishing an equivalent magnetic circuit model of MR damper, which has reverse exciting current in each coil. After integrating this formula with the MR damper parallel-plate model, the multiobjective optimization function and the working flow of optimal design are proposed. Thereafter, the comprehensive performance of MR dampers is further optimized by using PlatEMO [22], which is a MATLAB Platform for evolutionary multiobjective optimization. The results show that the Pareto solution set formed by the optimization method provides better solutions for those with performance indexes (damping force, dynamic range, and response time) higher than those of SG-MRD60 scheme. And the performance index of CMDG decreases by 11.27%, compared to the corresponding result of SG-MRD60 scheme. According to the finite element analysis, these results demonstrate that the introduction of CMDG into multiobjective optimization design for MR damper is of significant importance.


## 2. Performance Calculation Model

The performance calculation function mainly reflects the functional relation between the design variables and performance indexes of MR dampers. In order to optimize and compare the design variables, we select SG-MRD60 two-stage coil-type MR dampers manufactured by NingBo ShanGong Intelligent Safety Technology Co., Ltd. as the subject (see in Figure1). The SG-MRD60 has a very high degree of representativeness in engineering, and its entity structure model is shown in Figure 2. As for materials, according to the SG-MRD60, the piston rod, piston, and cylinder block are chosen by 45# steel, electrician pure iron, and Q235, respectively. The diameter of the cylinder cavity is 60 mm. The maximum stroke of the damper is ±40 mm. The rated input current is 0–2A, and the maximum damping force is 5 kN. On the premise of maintaining the basic structure, material, and cavity dimension of the damper, the optimization of the internal geometric dimension parameters is studied. The design variables to be optimized include diameter of the piston DPis, diameter of the rod DRod, length of the coil slot LCoil, depth of the coil groove HCoil, thickness of the cylinder tCase, and thickness of the flank piston W1, which are shown in Figure 3.Figure 1
SG-MRD60 and its application.Figure 2
MR damper entity structure model.Figure 3
Schematic diagram of design variables.


### 2.1. Design Theory of Magnetic Circuit Performance

#### 2.1.1. Response Time Model

As shown in Figure3, the magnetic circuit area is divided into piston axial area (Area11), rod area (Area12), piston flank radial area (Area21), piston middle radial area (Area22), flank damping gap area (Areaf1), middle damping gap area (Areaf2), and cylinder area (Area3). The distribution of the corresponding magnetic circuits and areas is shown in Figure 4. In this figure, the red arrows represent the magnetic induction line direction of the magnetic circuit when the MR damper coil works with reverse exciting current.Figure 4
Distribution of magnetic circuits and areas.According to Gauss’s law [23], the sum of magnetic flux of the circuit ∅sum can be written as(1)∅sum=∫sBndS=0,where Bn is the magnetic flux density and varies with the change of area and S is the cross-sectional area formed by the passage of magnetic induction lines through various parts of the area.On the basis of Figure4, an equivalent magnetic circuit model (shown in Figure 5) is established. In this model, the magnetic reluctance of most regions is connected in series. However, since the piston rod and piston are made from different materials, the magnetic reluctance in the corresponding areas “Area11” and “Area12” is connected in parallel.Figure 5
Equivalent magnetic circuit model.According to the model in Figure5, neglecting the influence of magnetic flux leakage, the magnetic reluctances of each area can be written as(2a)Rf1=lnDPis+2h/DPis2πμoμMRW1,(2b)Rf2=lnDPis+2h/DPisπμoμMRW2,(2c)Rs11=2LPisπDPis−2HCoil2−DRod2μ0μDT4,(2d)Rs12=2LPisπDRod2μ0μ45,(2e)Rs21=lnDPis/DRod2πμ0μDT4W1,(2f)Rs22=lnDPis/DRodπμ0μDT4W2,(2g)Rs3=2LPisπDPis+tCase2−DPis2μ0μ45,where Rf1, Rf2, Rs11, Rs12, Rs21, Rs22, and Rs3 represent the magnetic reluctance of the flank damping gap, middle damping gap, piston axial, rod axial, flank piston radial, middle piston radial, and damper cylinder, respectively. The meaning of other parameters appearing in the formula is given in Table 1.Table 1
Temperature and wildlife count in the three areas covered by the study.
Symbol
Parameter
Value
D
c
Cavity diameter
60 mm
L
Pis
Piston length
120 mm
B
max
Magnetic flux density
0.34 T
η
Viscosity of the MR fluid
0.24 Pa s
v
0
Piston linear velocity
0.05 m/s
r
Wire radius
7.5
×
10
−
4 m
ρ
Resistivity
1.75
×
10
−
8 Ω·m
I
max
Rated current
2 A
μ
0
Permeability of vacuum
1.257
×
10
−
6 H/m
μ
MR
Relative permeability of MR fluid
4
μ
DT
4
Relative permeability of DT4
1600
μ
45
Relative permeability of 45# steel
1000
μ
Coil
Relative permeability of coil
1
α
Magnetic leakage coefficient
0.1538According to Ohm’s law and Ampere’s law of magnetic circuit [23], the number of coil turns can be expressed as(3)NI=Rm∅=Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅,(4)Rs1=Rs11Rs12Rs11+Rs12,where N is the number of single stage coil and I is the current of coil. Due to the fact that the piston is symmetrically arranged with double coils, Rm represents the total magnetic reluctance of one coil and Rs1 is the parallel magnetic reluctance. The magnetic flux of flank damping gap ∅ in (3) is given by(5)∅=BW1πDPis+h,where B is the magnetic flux density of flank damping gap.According to the rated current range of MR dampers, letI=2.0 A; the number of single-stage coil can be rewritten as(6)N=0.5Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅.According to the ratio of coil turns to coil slot length, the value of coil layern is expressed as(7)n=2rNLCoil,nεN+.The coil resistance of the MR damperR can be expressed as(8)R=ρlπr2,where ρ is the resistivity, r is the wire radius, and l is the wire length of coil and it is given by(9)l=2πnLCoilDPis−2HCoil+2rn.The response time of MR damperT can be expressed as(10)T=2N∅RImax.


#### 2.1.2. CMDG Model

According to the structure of MR dampers with a two-stage coil, the damping gap can be divided into three effective parts in magnetic field. And its magnetic flux equals the sum of each coil. When the current directions of the two exciting coils are opposite to each other, the middle part has the highest magnetic flux [24]. Therefore, the consistency of magnetic flux density between damping gaps in different areas (such as Areaf1 and Areaf2, as shown in Figure 4) should be taken into account during the design process.The key to control CMDG is to determine the geometric dimensions ofW1 and W2 (as shown in Figure 3). The traditional model [13] simply defines the size of W2 as twice as W1, which affects the magnetic circuit performance of MRD. Therefore, this study does not adopt the traditional way and deduces the optimization objective function of CMDG. The derivation process is as follows:According to Gauss’s law and Ohm’s law of magnetic circuit [23], taking the left piston exciting coil as the research subject, the formulas for calculating the magnetic flux generated by the coil in the middle of the piston and the right side of the piston are expressed, respectively, as(11a)∅=NIRm,(11b)∅mid=NIRm2Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(11c)∅right=NIRmRf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where ∅mid and ∅right are the magnetic flux of middle and right damping gap, respectively.Because of the symmetry of the exciting coils in both sides of the piston, the magnetic flux density of the damping gap in each side is the same. Therefore, the left excitation coil of the piston is taken as the subject for research. According to the superposition principle of magnetic field, ignoring the direction of the magnetic field, the formula for calculating the magnetic flux density of the flank and the middle damping gap can be written as(12)Bf=NIRmW1πDPis+h1−Rf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(13)Bm=NIRmW2πDPis+h4Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where Bf and Bm are the magnetic flux density of flank and middle of damping gap, respectively.Combining (12) and (13), the formula of CMDG C is given by(14)C=Bf−BmmaxBf,Bm+α=W2−2W1maxW2,2W1+α×100%,where α is the magnetic leakage coefficient, which is given by(15)α=Bf^−Bm^maxBf^,Bm^,where Bf^ and Bm^ are the simulation results corresponding to Bf and Bm in physical meaning, respectively.


### 2.2. Design Theory of Mechanical Performance

On the basis of the MR damper parallel-plate model [25, 26], damping force can be expressed as(16)F=Fη+Fτ+Ff,where Fη is the viscous force, Fτ is the controllable force, and Ff is the friction. Considering that friction force is hard to model and calculate, we assume that there is no friction force:(17)Fη=1+whv02Q12ηQLPisApwh3,(18)Fτ=12Qη12Qη+0.4wh2τ0+2.07τ0LaAphsgnv0,where h is the dimension of damping gap, w is the mean circumference of the damper’s annular flow path, La is the effective axial pole length, Ap is the effective extrusion surface on the side of the piston, and Q is the volumetric flow rate; these variables are given by(19)h=0.560−DPis,w=πDPis+h,La=LPis−2LCoil,Ap=0.25πDPis2−DRod2,Q=Ap×v0.Combining (16)–(19), the damping force of MR damper F can be rewritten as(20)F=1+whv02Q12ηQLPisApwh3+12Qη12Qη+0.4wh2τ0+2.07τ0LaAphsgnv0.The dynamic rangeK is defined as the ratio of the damping force to the viscous force, and it is given by(21)K=FFη=FtFη+1.


## 2.1. Design Theory of Magnetic Circuit Performance

### 2.1.1. Response Time Model

As shown in Figure3, the magnetic circuit area is divided into piston axial area (Area11), rod area (Area12), piston flank radial area (Area21), piston middle radial area (Area22), flank damping gap area (Areaf1), middle damping gap area (Areaf2), and cylinder area (Area3). The distribution of the corresponding magnetic circuits and areas is shown in Figure 4. In this figure, the red arrows represent the magnetic induction line direction of the magnetic circuit when the MR damper coil works with reverse exciting current.Figure 4
Distribution of magnetic circuits and areas.According to Gauss’s law [23], the sum of magnetic flux of the circuit ∅sum can be written as(1)∅sum=∫sBndS=0,where Bn is the magnetic flux density and varies with the change of area and S is the cross-sectional area formed by the passage of magnetic induction lines through various parts of the area.On the basis of Figure4, an equivalent magnetic circuit model (shown in Figure 5) is established. In this model, the magnetic reluctance of most regions is connected in series. However, since the piston rod and piston are made from different materials, the magnetic reluctance in the corresponding areas “Area11” and “Area12” is connected in parallel.Figure 5
Equivalent magnetic circuit model.According to the model in Figure5, neglecting the influence of magnetic flux leakage, the magnetic reluctances of each area can be written as(2a)Rf1=lnDPis+2h/DPis2πμoμMRW1,(2b)Rf2=lnDPis+2h/DPisπμoμMRW2,(2c)Rs11=2LPisπDPis−2HCoil2−DRod2μ0μDT4,(2d)Rs12=2LPisπDRod2μ0μ45,(2e)Rs21=lnDPis/DRod2πμ0μDT4W1,(2f)Rs22=lnDPis/DRodπμ0μDT4W2,(2g)Rs3=2LPisπDPis+tCase2−DPis2μ0μ45,where Rf1, Rf2, Rs11, Rs12, Rs21, Rs22, and Rs3 represent the magnetic reluctance of the flank damping gap, middle damping gap, piston axial, rod axial, flank piston radial, middle piston radial, and damper cylinder, respectively. The meaning of other parameters appearing in the formula is given in Table 1.Table 1
Temperature and wildlife count in the three areas covered by the study.
Symbol
Parameter
Value
D
c
Cavity diameter
60 mm
L
Pis
Piston length
120 mm
B
max
Magnetic flux density
0.34 T
η
Viscosity of the MR fluid
0.24 Pa s
v
0
Piston linear velocity
0.05 m/s
r
Wire radius
7.5
×
10
−
4 m
ρ
Resistivity
1.75
×
10
−
8 Ω·m
I
max
Rated current
2 A
μ
0
Permeability of vacuum
1.257
×
10
−
6 H/m
μ
MR
Relative permeability of MR fluid
4
μ
DT
4
Relative permeability of DT4
1600
μ
45
Relative permeability of 45# steel
1000
μ
Coil
Relative permeability of coil
1
α
Magnetic leakage coefficient
0.1538According to Ohm’s law and Ampere’s law of magnetic circuit [23], the number of coil turns can be expressed as(3)NI=Rm∅=Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅,(4)Rs1=Rs11Rs12Rs11+Rs12,where N is the number of single stage coil and I is the current of coil. Due to the fact that the piston is symmetrically arranged with double coils, Rm represents the total magnetic reluctance of one coil and Rs1 is the parallel magnetic reluctance. The magnetic flux of flank damping gap ∅ in (3) is given by(5)∅=BW1πDPis+h,where B is the magnetic flux density of flank damping gap.According to the rated current range of MR dampers, letI=2.0 A; the number of single-stage coil can be rewritten as(6)N=0.5Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅.According to the ratio of coil turns to coil slot length, the value of coil layern is expressed as(7)n=2rNLCoil,nεN+.The coil resistance of the MR damperR can be expressed as(8)R=ρlπr2,where ρ is the resistivity, r is the wire radius, and l is the wire length of coil and it is given by(9)l=2πnLCoilDPis−2HCoil+2rn.The response time of MR damperT can be expressed as(10)T=2N∅RImax.


### 2.1.2. CMDG Model

According to the structure of MR dampers with a two-stage coil, the damping gap can be divided into three effective parts in magnetic field. And its magnetic flux equals the sum of each coil. When the current directions of the two exciting coils are opposite to each other, the middle part has the highest magnetic flux [24]. Therefore, the consistency of magnetic flux density between damping gaps in different areas (such as Areaf1 and Areaf2, as shown in Figure 4) should be taken into account during the design process.The key to control CMDG is to determine the geometric dimensions ofW1 and W2 (as shown in Figure 3). The traditional model [13] simply defines the size of W2 as twice as W1, which affects the magnetic circuit performance of MRD. Therefore, this study does not adopt the traditional way and deduces the optimization objective function of CMDG. The derivation process is as follows:According to Gauss’s law and Ohm’s law of magnetic circuit [23], taking the left piston exciting coil as the research subject, the formulas for calculating the magnetic flux generated by the coil in the middle of the piston and the right side of the piston are expressed, respectively, as(11a)∅=NIRm,(11b)∅mid=NIRm2Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(11c)∅right=NIRmRf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where ∅mid and ∅right are the magnetic flux of middle and right damping gap, respectively.Because of the symmetry of the exciting coils in both sides of the piston, the magnetic flux density of the damping gap in each side is the same. Therefore, the left excitation coil of the piston is taken as the subject for research. According to the superposition principle of magnetic field, ignoring the direction of the magnetic field, the formula for calculating the magnetic flux density of the flank and the middle damping gap can be written as(12)Bf=NIRmW1πDPis+h1−Rf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(13)Bm=NIRmW2πDPis+h4Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where Bf and Bm are the magnetic flux density of flank and middle of damping gap, respectively.Combining (12) and (13), the formula of CMDG C is given by(14)C=Bf−BmmaxBf,Bm+α=W2−2W1maxW2,2W1+α×100%,where α is the magnetic leakage coefficient, which is given by(15)α=Bf^−Bm^maxBf^,Bm^,where Bf^ and Bm^ are the simulation results corresponding to Bf and Bm in physical meaning, respectively.


## 2.1.1. Response Time Model

As shown in Figure3, the magnetic circuit area is divided into piston axial area (Area11), rod area (Area12), piston flank radial area (Area21), piston middle radial area (Area22), flank damping gap area (Areaf1), middle damping gap area (Areaf2), and cylinder area (Area3). The distribution of the corresponding magnetic circuits and areas is shown in Figure 4. In this figure, the red arrows represent the magnetic induction line direction of the magnetic circuit when the MR damper coil works with reverse exciting current.Figure 4
Distribution of magnetic circuits and areas.According to Gauss’s law [23], the sum of magnetic flux of the circuit ∅sum can be written as(1)∅sum=∫sBndS=0,where Bn is the magnetic flux density and varies with the change of area and S is the cross-sectional area formed by the passage of magnetic induction lines through various parts of the area.On the basis of Figure4, an equivalent magnetic circuit model (shown in Figure 5) is established. In this model, the magnetic reluctance of most regions is connected in series. However, since the piston rod and piston are made from different materials, the magnetic reluctance in the corresponding areas “Area11” and “Area12” is connected in parallel.Figure 5
Equivalent magnetic circuit model.According to the model in Figure5, neglecting the influence of magnetic flux leakage, the magnetic reluctances of each area can be written as(2a)Rf1=lnDPis+2h/DPis2πμoμMRW1,(2b)Rf2=lnDPis+2h/DPisπμoμMRW2,(2c)Rs11=2LPisπDPis−2HCoil2−DRod2μ0μDT4,(2d)Rs12=2LPisπDRod2μ0μ45,(2e)Rs21=lnDPis/DRod2πμ0μDT4W1,(2f)Rs22=lnDPis/DRodπμ0μDT4W2,(2g)Rs3=2LPisπDPis+tCase2−DPis2μ0μ45,where Rf1, Rf2, Rs11, Rs12, Rs21, Rs22, and Rs3 represent the magnetic reluctance of the flank damping gap, middle damping gap, piston axial, rod axial, flank piston radial, middle piston radial, and damper cylinder, respectively. The meaning of other parameters appearing in the formula is given in Table 1.Table 1
Temperature and wildlife count in the three areas covered by the study.
Symbol
Parameter
Value
D
c
Cavity diameter
60 mm
L
Pis
Piston length
120 mm
B
max
Magnetic flux density
0.34 T
η
Viscosity of the MR fluid
0.24 Pa s
v
0
Piston linear velocity
0.05 m/s
r
Wire radius
7.5
×
10
−
4 m
ρ
Resistivity
1.75
×
10
−
8 Ω·m
I
max
Rated current
2 A
μ
0
Permeability of vacuum
1.257
×
10
−
6 H/m
μ
MR
Relative permeability of MR fluid
4
μ
DT
4
Relative permeability of DT4
1600
μ
45
Relative permeability of 45# steel
1000
μ
Coil
Relative permeability of coil
1
α
Magnetic leakage coefficient
0.1538According to Ohm’s law and Ampere’s law of magnetic circuit [23], the number of coil turns can be expressed as(3)NI=Rm∅=Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅,(4)Rs1=Rs11Rs12Rs11+Rs12,where N is the number of single stage coil and I is the current of coil. Due to the fact that the piston is symmetrically arranged with double coils, Rm represents the total magnetic reluctance of one coil and Rs1 is the parallel magnetic reluctance. The magnetic flux of flank damping gap ∅ in (3) is given by(5)∅=BW1πDPis+h,where B is the magnetic flux density of flank damping gap.According to the rated current range of MR dampers, letI=2.0 A; the number of single-stage coil can be rewritten as(6)N=0.5Rf1+Rf2+Rs1+Rs21+Rs22+Rs3∅.According to the ratio of coil turns to coil slot length, the value of coil layern is expressed as(7)n=2rNLCoil,nεN+.The coil resistance of the MR damperR can be expressed as(8)R=ρlπr2,where ρ is the resistivity, r is the wire radius, and l is the wire length of coil and it is given by(9)l=2πnLCoilDPis−2HCoil+2rn.The response time of MR damperT can be expressed as(10)T=2N∅RImax.


## 2.1.2. CMDG Model

According to the structure of MR dampers with a two-stage coil, the damping gap can be divided into three effective parts in magnetic field. And its magnetic flux equals the sum of each coil. When the current directions of the two exciting coils are opposite to each other, the middle part has the highest magnetic flux [24]. Therefore, the consistency of magnetic flux density between damping gaps in different areas (such as Areaf1 and Areaf2, as shown in Figure 4) should be taken into account during the design process.The key to control CMDG is to determine the geometric dimensions ofW1 and W2 (as shown in Figure 3). The traditional model [13] simply defines the size of W2 as twice as W1, which affects the magnetic circuit performance of MRD. Therefore, this study does not adopt the traditional way and deduces the optimization objective function of CMDG. The derivation process is as follows:According to Gauss’s law and Ohm’s law of magnetic circuit [23], taking the left piston exciting coil as the research subject, the formulas for calculating the magnetic flux generated by the coil in the middle of the piston and the right side of the piston are expressed, respectively, as(11a)∅=NIRm,(11b)∅mid=NIRm2Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(11c)∅right=NIRmRf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where ∅mid and ∅right are the magnetic flux of middle and right damping gap, respectively.Because of the symmetry of the exciting coils in both sides of the piston, the magnetic flux density of the damping gap in each side is the same. Therefore, the left excitation coil of the piston is taken as the subject for research. According to the superposition principle of magnetic field, ignoring the direction of the magnetic field, the formula for calculating the magnetic flux density of the flank and the middle damping gap can be written as(12)Bf=NIRmW1πDPis+h1−Rf2+Rs222Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,(13)Bm=NIRmW2πDPis+h4Rf1+Rs1+Rs21+Rs32Rf1+Rs1+Rs21+Rs3+Rf2+Rs22,where Bf and Bm are the magnetic flux density of flank and middle of damping gap, respectively.Combining (12) and (13), the formula of CMDG C is given by(14)C=Bf−BmmaxBf,Bm+α=W2−2W1maxW2,2W1+α×100%,where α is the magnetic leakage coefficient, which is given by(15)α=Bf^−Bm^maxBf^,Bm^,where Bf^ and Bm^ are the simulation results corresponding to Bf and Bm in physical meaning, respectively.


## 2.2. Design Theory of Mechanical Performance

On the basis of the MR damper parallel-plate model [25, 26], damping force can be expressed as(16)F=Fη+Fτ+Ff,where Fη is the viscous force, Fτ is the controllable force, and Ff is the friction. Considering that friction force is hard to model and calculate, we assume that there is no friction force:(17)Fη=1+whv02Q12ηQLPisApwh3,(18)Fτ=12Qη12Qη+0.4wh2τ0+2.07τ0LaAphsgnv0,where h is the dimension of damping gap, w is the mean circumference of the damper’s annular flow path, La is the effective axial pole length, Ap is the effective extrusion surface on the side of the piston, and Q is the volumetric flow rate; these variables are given by(19)h=0.560−DPis,w=πDPis+h,La=LPis−2LCoil,Ap=0.25πDPis2−DRod2,Q=Ap×v0.Combining (16)–(19), the damping force of MR damper F can be rewritten as(20)F=1+whv02Q12ηQLPisApwh3+12Qη12Qη+0.4wh2τ0+2.07τ0LaAphsgnv0.The dynamic rangeK is defined as the ratio of the damping force to the viscous force, and it is given by(21)K=FFη=FtFη+1.


## 3. Optimization Design Method

### 3.1. Background of Multiobject Optimization

A multiobjective optimization problem, as the name suggests, has many objective functions. Its general form [27] can be expressed as(22)min/maxfmxi,m=1,2,…,Ms.t.gjxi≥0,j=1,2,…,Js.t.hkxi=0,k=1,2,…,Ks.t.xil≤xi≤xiui=1,2,…,n,where fmxi is the objective function and gjxi and hkxi are the constraint function. M, J, and K are the numbers of each kind of function. xi is a vector of n decision variables, whose range of value is limited within a lower bound xil and an upper bound xiu.In this procedure of optimization design,fmxi represents four objective functions, the damping force, the dynamic range, the response time, and the CMDG. The vector xi represents all the parameters of geometric dimension. If we convert all the fmxi into a single-objective optimization problem, the objective function will be written as(23)Fxi=∑m=1Mwmfmxi,where Fxi is the objective function for the problem of single-objective optimization and wm is the weighting factor for each objective function.As shown in (23), the result of a single-objective optimization depends on the value of wm and there is only one optimal solution which can be obtained after optimization. It is impossible for the decision-maker to compare the solutions in the optimal solution set and select the required optimal solution from them. Therefore, the appropriate evolutionary algorithm needs to be selected and used to obtain the optimal solutions from the feasible solutions, which is the other focus in this paper.


### 3.2. Premise and Hypothesis

(1)
The cavity diameterDc (Dc=DPis+2h) is determined as 60 mm.
(2)
In order to keep the piston stroke constant when the length of the cavity is constant, let piston lengthLPis=120 mm.
(3)
According to the measured data of SG-MRF2035MR fluid, the relationship between yield stress and magnetic flux density is shown in Figure6, and the curve of the yield stress τ0 is fitted by the least-squares method and can be written as
(24)
τ
0
=
35.3
B
4
−
77.6
B
3
−
3.5
B
2
+
105.2
B
−
6.8.
According to the relevant parameters of SG-MRD60 under the maximum operating currentI = 2 A, we assume that the maximum magnetic flux density of damping gap Bmax=0.34 T. Then, the yield stress of MR fluid in damping gap can be calculated by (20) with a value of 26 kPa.
(4)
Substituting the finite element simulation results of SG-MRD60 into (15), we get the magnetic leakage coefficient α=0.1538.
(5)
Set the piston linear velocityv0=50 mm/s, the optimization parameters are listed in Table 1.Figure 6
Relationship betweenτ0 and B.


### 3.3. Relative Permeability of Materials

In fact, the permeability of the material tends to decrease with the increase of magnetic field strength, and different materials have different relative permeabilities at different temperatures [14]. Generally, the magnetization curve of DT4 pure iron and 45# steel can be seen in Figure 7, which shows the nonlinearity relationship between B and H. Since the relative permeability is the derivative of H to B, it can be approximated to a constant when the magnetic flux density B of the materials is less than 1.5 T.Figure 7
Relationship betweenB and H [26].To simplify the calculation, it is assumed that the relative permeability of each material is a constant which does not vary with temperature. And the values of relative permeability [26] of MR fluid (μMR), DT4 pure iron (μDT4), 45# steel (μ45), and Coil (μCoil) are listed in Table 1.


### 3.4. Constraint Condition

(1)
Diameter of the piston must be smaller than diameter of the inner cavity, that is,DPis<60mm.
(2)
According to the GB/T699-1999 standard, 45# steel with a diameter of 15 mm can withstand the static axial force of 106 kN, which is much larger than that of the SG-MRD60 design with a damping force 5 kN; therefore, letDRod>15 mm.
(3)
According to MR damper geometric relations, the constraint conditions of the geometric dimensionW2 and h can be expressed as
(25)
W
2
=
L
Pis
−
2
L
Coil
−
2
W
1
,
(26)
h
=
60
−
D
Pis
2
.
(4)
The remaining constraints are selected with reference to the SG-MRD60 geometric dimensions.Based on the geometric relationship between the constraint conditions and the geometric dimension parameters, six geometric dimension parameters are determined as design variables to be optimized. The corresponding meaning and range of values are shown in Table2.Table 2
Parameters of geometric dimension.
Symbol
Parameter
Range
Unit
D
Pis
Diameter of the piston
[54, 60]
mm
D
Rod
Diameter of the rod
[15, 20]
mm
L
Coil
Length of the coil slot
[35, 45]
mm
H
Coil
Depth of the coil groove
[14, 16]
mm
h
Width of the damping gap
—
mm
t
Case
Thickness of the cylinder body
[6, 9]
mm
W
1
Thickness of the flank piston
[8, 12]
mm
W
2
Thickness of the middle piston
—
mm


### 3.5. Algorithm Principle and Optimization Method

In multiobjective optimization problem, it is difficult to achieve the best optimal results simultaneously for all optimal objectives. The NSGA-III algorithm [21] is used in this study. This algorithm seeks to find all solutions with at least one target having better solution while making no target worse, namely, Pareto-optimal solutions. Each performance index is overall optimal in the front surface formed by the Pareto-optimal solution set, which is convenient for decision-makers to filter out the optimal solution to meet their own needs.


#### 3.5.1. Principle of NSGA-III Algorithm

In order to handle the multiobjective optimization problems, NSGA-III algorithm [21] is proposed as an evolutionary multiobjective optimization, which, compared with the NSGA-II algorithm [28], achieves further improvement of population selection mechanism. NSGA-III algorithm can ensure diversity of obtained solutions through a predefined set of reference points. As for the operating steps of the NSGA-III algorithm, it can be written as follows. Firstly, we use the genetic algorithm to initialize the population and layer the population at different levels according to the nondomination sorting algorithm, so a new population St can be sorted from each level of the population. Secondly, on the basis of the minimum value of each objective function and the reference point supplied by a systematic approach or user, we construct the reference line and associate the population St with it. According to the niching algorithm, niche counts that determine the next generation of populations are obtained by calculating the distance between population St and reference line. Finally, the iterative process is repeated and the algorithm does not stop until it reaches the maximum iterations. When the algorithm stops, it outputs the Pareto-optimal solutions.


#### 3.5.2. Procedure of Optimization Design

According to the needs of engineering, the objectives of the damping force and the dynamic range are necessary to be maximized. However, the NSGA-III algorithm is used to minimize the objective function in the PlatEMO [22]. Therefore, we transfer the optimization problem of max fxi into the problem of min 1/fxi. Taking the damping force, the dynamic range, the response time, and the CMDG as the optimization objectives simultaneously, as is mentioned in Section 3.1, the multiobjective optimization problem can be rewritten as(27)f1xi=1Fxi,f2xi=1Kxi,f3xi=Txi,f4xi=Cxi,(28)minimizeFxi=f1xi,f2xi,f3xi,f4xiT,xi=DPis,i,DRod,i,LCoil,i,HCoil,i,tCase,i,W1,iT,where f1xi, f2xi, f3xi, and f4xi are the optimization objective function about the damping force, the dynamic range, the response time, and the CMDG, respectively.Equation (28) is written to a MRD file, and the above optimization algorithm is implemented by setting parameters (shown as Table 3) in the PlatEMO [22], a MATLAB-based open source platform with more than 50 popular evolutionary multiobjective optimization algorithms including the NSGA-III algorithm, is easily used by the MATLAB GUI.Table 3
Acceptable parameters for PlatEMO.
Parameter name
Description
Value
-algorithm
MOEA function
@NSGA-III
-problem
MOP function
@MRD
-N
Population size
6000
-M
Number of objectives
3
-D
Number of variables
6
-evaluation
Number of evaluations
30,000At the end of each iteration in the running NSGA-III algorithm, it is necessary to judge all the design variables whether the depth of the coil grooveHCoil is equal to or greater than the product of the number of coil layers n and the wire radius r. If a set of design variables satisfies HCoil<n×2r, the set of design variables will be defined as an invalid variable and eliminated at the next iteration of NSGA-III algorithm.After the calculation of NSGA-III algorithm, the variablePopulation including both the Pareto-optimal solutions and its geometric dimensions is outputted by the PlatEMO. And then, the scheme index which is required by engineering will be filtered from Population by using Excel. Besides, the corresponding geometric dimension vectors are found from Population. The working flow of optimal design is shown in Figure 8.Figure 8
Working flow of optimal design.


## 3.1. Background of Multiobject Optimization

A multiobjective optimization problem, as the name suggests, has many objective functions. Its general form [27] can be expressed as(22)min/maxfmxi,m=1,2,…,Ms.t.gjxi≥0,j=1,2,…,Js.t.hkxi=0,k=1,2,…,Ks.t.xil≤xi≤xiui=1,2,…,n,where fmxi is the objective function and gjxi and hkxi are the constraint function. M, J, and K are the numbers of each kind of function. xi is a vector of n decision variables, whose range of value is limited within a lower bound xil and an upper bound xiu.In this procedure of optimization design,fmxi represents four objective functions, the damping force, the dynamic range, the response time, and the CMDG. The vector xi represents all the parameters of geometric dimension. If we convert all the fmxi into a single-objective optimization problem, the objective function will be written as(23)Fxi=∑m=1Mwmfmxi,where Fxi is the objective function for the problem of single-objective optimization and wm is the weighting factor for each objective function.As shown in (23), the result of a single-objective optimization depends on the value of wm and there is only one optimal solution which can be obtained after optimization. It is impossible for the decision-maker to compare the solutions in the optimal solution set and select the required optimal solution from them. Therefore, the appropriate evolutionary algorithm needs to be selected and used to obtain the optimal solutions from the feasible solutions, which is the other focus in this paper.


## 3.2. Premise and Hypothesis

(1)
The cavity diameterDc (Dc=DPis+2h) is determined as 60 mm.
(2)
In order to keep the piston stroke constant when the length of the cavity is constant, let piston lengthLPis=120 mm.
(3)
According to the measured data of SG-MRF2035MR fluid, the relationship between yield stress and magnetic flux density is shown in Figure6, and the curve of the yield stress τ0 is fitted by the least-squares method and can be written as
(24)
τ
0
=
35.3
B
4
−
77.6
B
3
−
3.5
B
2
+
105.2
B
−
6.8.
According to the relevant parameters of SG-MRD60 under the maximum operating currentI = 2 A, we assume that the maximum magnetic flux density of damping gap Bmax=0.34 T. Then, the yield stress of MR fluid in damping gap can be calculated by (20) with a value of 26 kPa.
(4)
Substituting the finite element simulation results of SG-MRD60 into (15), we get the magnetic leakage coefficient α=0.1538.
(5)
Set the piston linear velocityv0=50 mm/s, the optimization parameters are listed in Table 1.Figure 6
Relationship betweenτ0 and B.


## 3.3. Relative Permeability of Materials

In fact, the permeability of the material tends to decrease with the increase of magnetic field strength, and different materials have different relative permeabilities at different temperatures [14]. Generally, the magnetization curve of DT4 pure iron and 45# steel can be seen in Figure 7, which shows the nonlinearity relationship between B and H. Since the relative permeability is the derivative of H to B, it can be approximated to a constant when the magnetic flux density B of the materials is less than 1.5 T.Figure 7
Relationship betweenB and H [26].To simplify the calculation, it is assumed that the relative permeability of each material is a constant which does not vary with temperature. And the values of relative permeability [26] of MR fluid (μMR), DT4 pure iron (μDT4), 45# steel (μ45), and Coil (μCoil) are listed in Table 1.


## 3.4. Constraint Condition

(1)
Diameter of the piston must be smaller than diameter of the inner cavity, that is,DPis<60mm.
(2)
According to the GB/T699-1999 standard, 45# steel with a diameter of 15 mm can withstand the static axial force of 106 kN, which is much larger than that of the SG-MRD60 design with a damping force 5 kN; therefore, letDRod>15 mm.
(3)
According to MR damper geometric relations, the constraint conditions of the geometric dimensionW2 and h can be expressed as
(25)
W
2
=
L
Pis
−
2
L
Coil
−
2
W
1
,
(26)
h
=
60
−
D
Pis
2
.
(4)
The remaining constraints are selected with reference to the SG-MRD60 geometric dimensions.Based on the geometric relationship between the constraint conditions and the geometric dimension parameters, six geometric dimension parameters are determined as design variables to be optimized. The corresponding meaning and range of values are shown in Table2.Table 2
Parameters of geometric dimension.
Symbol
Parameter
Range
Unit
D
Pis
Diameter of the piston
[54, 60]
mm
D
Rod
Diameter of the rod
[15, 20]
mm
L
Coil
Length of the coil slot
[35, 45]
mm
H
Coil
Depth of the coil groove
[14, 16]
mm
h
Width of the damping gap
—
mm
t
Case
Thickness of the cylinder body
[6, 9]
mm
W
1
Thickness of the flank piston
[8, 12]
mm
W
2
Thickness of the middle piston
—
mm


## 3.5. Algorithm Principle and Optimization Method

In multiobjective optimization problem, it is difficult to achieve the best optimal results simultaneously for all optimal objectives. The NSGA-III algorithm [21] is used in this study. This algorithm seeks to find all solutions with at least one target having better solution while making no target worse, namely, Pareto-optimal solutions. Each performance index is overall optimal in the front surface formed by the Pareto-optimal solution set, which is convenient for decision-makers to filter out the optimal solution to meet their own needs.


### 3.5.1. Principle of NSGA-III Algorithm

In order to handle the multiobjective optimization problems, NSGA-III algorithm [21] is proposed as an evolutionary multiobjective optimization, which, compared with the NSGA-II algorithm [28], achieves further improvement of population selection mechanism. NSGA-III algorithm can ensure diversity of obtained solutions through a predefined set of reference points. As for the operating steps of the NSGA-III algorithm, it can be written as follows. Firstly, we use the genetic algorithm to initialize the population and layer the population at different levels according to the nondomination sorting algorithm, so a new population St can be sorted from each level of the population. Secondly, on the basis of the minimum value of each objective function and the reference point supplied by a systematic approach or user, we construct the reference line and associate the population St with it. According to the niching algorithm, niche counts that determine the next generation of populations are obtained by calculating the distance between population St and reference line. Finally, the iterative process is repeated and the algorithm does not stop until it reaches the maximum iterations. When the algorithm stops, it outputs the Pareto-optimal solutions.


### 3.5.2. Procedure of Optimization Design

According to the needs of engineering, the objectives of the damping force and the dynamic range are necessary to be maximized. However, the NSGA-III algorithm is used to minimize the objective function in the PlatEMO [22]. Therefore, we transfer the optimization problem of max fxi into the problem of min 1/fxi. Taking the damping force, the dynamic range, the response time, and the CMDG as the optimization objectives simultaneously, as is mentioned in Section 3.1, the multiobjective optimization problem can be rewritten as(27)f1xi=1Fxi,f2xi=1Kxi,f3xi=Txi,f4xi=Cxi,(28)minimizeFxi=f1xi,f2xi,f3xi,f4xiT,xi=DPis,i,DRod,i,LCoil,i,HCoil,i,tCase,i,W1,iT,where f1xi, f2xi, f3xi, and f4xi are the optimization objective function about the damping force, the dynamic range, the response time, and the CMDG, respectively.Equation (28) is written to a MRD file, and the above optimization algorithm is implemented by setting parameters (shown as Table 3) in the PlatEMO [22], a MATLAB-based open source platform with more than 50 popular evolutionary multiobjective optimization algorithms including the NSGA-III algorithm, is easily used by the MATLAB GUI.Table 3
Acceptable parameters for PlatEMO.
Parameter name
Description
Value
-algorithm
MOEA function
@NSGA-III
-problem
MOP function
@MRD
-N
Population size
6000
-M
Number of objectives
3
-D
Number of variables
6
-evaluation
Number of evaluations
30,000At the end of each iteration in the running NSGA-III algorithm, it is necessary to judge all the design variables whether the depth of the coil grooveHCoil is equal to or greater than the product of the number of coil layers n and the wire radius r. If a set of design variables satisfies HCoil<n×2r, the set of design variables will be defined as an invalid variable and eliminated at the next iteration of NSGA-III algorithm.After the calculation of NSGA-III algorithm, the variablePopulation including both the Pareto-optimal solutions and its geometric dimensions is outputted by the PlatEMO. And then, the scheme index which is required by engineering will be filtered from Population by using Excel. Besides, the corresponding geometric dimension vectors are found from Population. The working flow of optimal design is shown in Figure 8.Figure 8
Working flow of optimal design.


## 3.5.1. Principle of NSGA-III Algorithm

In order to handle the multiobjective optimization problems, NSGA-III algorithm [21] is proposed as an evolutionary multiobjective optimization, which, compared with the NSGA-II algorithm [28], achieves further improvement of population selection mechanism. NSGA-III algorithm can ensure diversity of obtained solutions through a predefined set of reference points. As for the operating steps of the NSGA-III algorithm, it can be written as follows. Firstly, we use the genetic algorithm to initialize the population and layer the population at different levels according to the nondomination sorting algorithm, so a new population St can be sorted from each level of the population. Secondly, on the basis of the minimum value of each objective function and the reference point supplied by a systematic approach or user, we construct the reference line and associate the population St with it. According to the niching algorithm, niche counts that determine the next generation of populations are obtained by calculating the distance between population St and reference line. Finally, the iterative process is repeated and the algorithm does not stop until it reaches the maximum iterations. When the algorithm stops, it outputs the Pareto-optimal solutions.


## 3.5.2. Procedure of Optimization Design

According to the needs of engineering, the objectives of the damping force and the dynamic range are necessary to be maximized. However, the NSGA-III algorithm is used to minimize the objective function in the PlatEMO [22]. Therefore, we transfer the optimization problem of max fxi into the problem of min 1/fxi. Taking the damping force, the dynamic range, the response time, and the CMDG as the optimization objectives simultaneously, as is mentioned in Section 3.1, the multiobjective optimization problem can be rewritten as(27)f1xi=1Fxi,f2xi=1Kxi,f3xi=Txi,f4xi=Cxi,(28)minimizeFxi=f1xi,f2xi,f3xi,f4xiT,xi=DPis,i,DRod,i,LCoil,i,HCoil,i,tCase,i,W1,iT,where f1xi, f2xi, f3xi, and f4xi are the optimization objective function about the damping force, the dynamic range, the response time, and the CMDG, respectively.Equation (28) is written to a MRD file, and the above optimization algorithm is implemented by setting parameters (shown as Table 3) in the PlatEMO [22], a MATLAB-based open source platform with more than 50 popular evolutionary multiobjective optimization algorithms including the NSGA-III algorithm, is easily used by the MATLAB GUI.Table 3
Acceptable parameters for PlatEMO.
Parameter name
Description
Value
-algorithm
MOEA function
@NSGA-III
-problem
MOP function
@MRD
-N
Population size
6000
-M
Number of objectives
3
-D
Number of variables
6
-evaluation
Number of evaluations
30,000At the end of each iteration in the running NSGA-III algorithm, it is necessary to judge all the design variables whether the depth of the coil grooveHCoil is equal to or greater than the product of the number of coil layers n and the wire radius r. If a set of design variables satisfies HCoil<n×2r, the set of design variables will be defined as an invalid variable and eliminated at the next iteration of NSGA-III algorithm.After the calculation of NSGA-III algorithm, the variablePopulation including both the Pareto-optimal solutions and its geometric dimensions is outputted by the PlatEMO. And then, the scheme index which is required by engineering will be filtered from Population by using Excel. Besides, the corresponding geometric dimension vectors are found from Population. The working flow of optimal design is shown in Figure 8.Figure 8
Working flow of optimal design.


## 4. Results and Discussion

### 4.1. Results of Optimization

#### 4.1.1. Results of Algorithm Calculation

The NSGA-III algorithm generates a total of 5986 Pareto-optimal solutions after 30000 iterations. The variables of geometric dimension corresponding to the Pareto-optimal solutions are taken into (10), (14), (20), and (21) respectively, and then, the abnormal solutions in the course of preserving the operation precision are removed. With damping force less than 35 kN, 4844 sets of data points are retained. The distribution of 3D Pareto-optimal front is drawn and shown as Figure 9. Then, we divide the Pareto-optimal solutions into 3 layers according to the range of response time and CMDG from the 4844 sets of data points. And the distribution of points is shown in Figures 10 and 11.Figure 9
3D Pareto-optimal solutions.Figure 10
Layer of Pareto-optimal solutions for response time.Figure 11
Layer of Pareto-optimal solutions for CMDG.As shown in Figure9, the damping force is negatively correlated with the dynamic range, which is in accordance with results of similar research [13]. According to Figure 10, the Pareto-optimal solutions are concentrated in a region with low damping force and high dynamic range when response time is less than 160 ms. With the increase of response time, the Pareto-optimal solutions tend to move to the region with high damping force and low dynamic range. And the comprehensive mechanical performance which refers to damping force and dynamic range will be improved. As shown in Figure 11, the design scheme with CMDG less than 0.2% has a large distribution area of Pareto-optimal solutions and the frontier surface formed by the Pareto solution set is uniformly distributed, which meets the needs of design. Besides, with the significant increase of CMDG, the improvement of mechanical performance of MR damper is not obvious.


#### 4.1.2. Analysis of Parametric Correlation

In order to analyze the uncertainty in MR damper parameters and to find which parameter will influence the performances of MR dampers most, we have done some research studies to obtain the correlation coefficient between each value of the optimal objective and the geometric dimension parameter in the Pareto-optimal solutions.The linear correlation coefficient is calculated by using the linear correlation analysis formula (29). The calculating results are listed in Table 4 and shown in Figure 12:(29)ρxy=n∑xiyi−∑xi∑yin∑xi2−∑xi2n∑yi2−∑yi2×100%,where ρxy is the linear correlation coefficient and xi and yi are the ith parameter of the optimal objective and the structural parameter, respectively.Table 4
Correlation coefficient matrix (%).
F
K
T
C
D
Pis
79.66
−96.13
21.70
1.25
D
Rod
−0.03
1.23
−4.60
−2.91
L
Coil
−13.87
−7.64
−54.55
−51.13
H
Coil
12.79
−5.28
47.52
5.37
t
Case
1.67
0.99
0.06
−0.57
w
1
6.68
4.15
94.53
−42.28Figure 12
Correlation coefficient between the optimal objective and the geometric dimension parameters. (a) The correlation betweenF and geometric dimensions. (b) The correlation between K and geometric dimensions. (c) The correlation between (T) and geometric dimensions. (d) The correlation between C and geometric dimensions.
(a)
(b)
(c)
(d)As shown in Figure12, the piston diameter DPis makes a great influence on the performances of the damping force and the dynamic range in all the geometric dimension parameters. And the flank piston thickness W1 has higher influence on the performances of the response time, and the CMDG also performs apparently than other parameters. Besides, the DRod and tCase are useless to improve the optimization objective of the MR dampers, as shown in Table 4. It may be partly because of the small values of geometric dimension in the rod diameter and the cylinder body thickness compared to other dimensions.


#### 4.1.3. Scheme Comparison

To compare the influence of CMDG difference, we filter out solutions where the performance indexes of the damping forceF, the dynamic range K, and the response time T are better than those of SG-MRD60 scheme from the Pareto-optimal solutions set by using Excel. According to the optimization results of CMDG, the solution whose CMDG is approximately 0%, 15%, and 30% are selected and named as scheme A, B, and C, respectively. All schemes including SG-MRD60 scheme are shown in Table 5.Table 5
Optimization results of partial geometry dimensions.
Parameter
SG-MRD60
Scheme A
Scheme B
Scheme C
F (kN)
5.03
5.14
5.05
5.27
K
32.24
33.69
43.03
41.73
T (ms)
215
190
182
166
C (%)
15.3
0.6
15.1
29.9
D
Pis (mm)
57.00
56.99
56.72
56.73
D
Rod (mm)
18.00
18.37
18.15
15.26
L
Coil (mm)
40.00
39.34
37.55
37.42
H
Coil (mm)
15.00
14.38
14.10
14.72
t
Case (mm)
7.00
8.08
6.92
7.79
W
1 (mm)
10.00
9.44
9.22
8.00
N
108
134
136
121According to Table5, compared with the SG-MRD60 scheme, the damping force, the dynamic range, and the response time in each scheme are increased by up to 2.2%, 4.5%, and 11.6%, respectively.


### 4.2. Results of Simulation

In order to compare the optimization results with the simulation results of CMDG, we import the parameters of geometric dimension (see in Table5) into ANSYS finite element analysis software. The Nephogram magnetic flux density distribution for each scheme is obtained and shown in Figure 13. And the simulation results of magnetic flux density of the damping gap are listed in Table 6. Besides, the average magnetic flux density of the damping gap Bf−m¯ is given by(30)Bf−m¯=Bf×2+Bm3.Figure 13
Nephogram magnetic flux density distribution obtained by FEA. (a) SG-MRD60 scheme. (b) Scheme A. (c) Scheme B. (d) scheme C.
(a)
(b)
(c)
(d)Table 6
Magnetic flux density result obtained by FEA.
Parameter
SG-MRD60
Scheme A
Scheme B
Scheme C
B
f (T)
0.242
0.326
0.333
0.312
B
m (T)
0.286
0.340
0.299
0.242
B
f
−
m
¯ (T)
0.257
0.331
0.332
0.289
C
ANSYS (%)
15.38
4.11
11.37
28.93As can be seen in Figure13(a)–13(c), the original scheme from SG-MRD60, whose magnetic flux density of middle damping gap is significantly greater than that of in the flank damping gap, has poor performance in CMDG. The CMDG of schemes A and B are better than that of SG-MRD60. It is in accordance with the result in Table 6. Besides, the trend of CMDG in Table 6 is the same as that in Table 5. The average deviation between the optimization result and the simulation result is 0.4%, which reflects that the CMDG function can express the distribution effect of the magnetic flux density in damping gap well.As shown in Table6, the CMDG in scheme C is 28.93%, which is the largest of all schemes. It leads to the consequence that the simulation result of CMDG in scheme C is the worst as shown in Figure 13(d), which means CMDG with maximum value has the worst performance. The results show that the worse the CMDG performance is, the more unbalanced the magnetic flux density of flank damping gap and middle damping gap will be. It will affect the calculation result of the optimization objective function of damping force and make the optimization scheme invalid. Therefore, the introduction of CMDG into multiobjective optimization design for MR dampers is very important.


## 4.1. Results of Optimization

### 4.1.1. Results of Algorithm Calculation

The NSGA-III algorithm generates a total of 5986 Pareto-optimal solutions after 30000 iterations. The variables of geometric dimension corresponding to the Pareto-optimal solutions are taken into (10), (14), (20), and (21) respectively, and then, the abnormal solutions in the course of preserving the operation precision are removed. With damping force less than 35 kN, 4844 sets of data points are retained. The distribution of 3D Pareto-optimal front is drawn and shown as Figure 9. Then, we divide the Pareto-optimal solutions into 3 layers according to the range of response time and CMDG from the 4844 sets of data points. And the distribution of points is shown in Figures 10 and 11.Figure 9
3D Pareto-optimal solutions.Figure 10
Layer of Pareto-optimal solutions for response time.Figure 11
Layer of Pareto-optimal solutions for CMDG.As shown in Figure9, the damping force is negatively correlated with the dynamic range, which is in accordance with results of similar research [13]. According to Figure 10, the Pareto-optimal solutions are concentrated in a region with low damping force and high dynamic range when response time is less than 160 ms. With the increase of response time, the Pareto-optimal solutions tend to move to the region with high damping force and low dynamic range. And the comprehensive mechanical performance which refers to damping force and dynamic range will be improved. As shown in Figure 11, the design scheme with CMDG less than 0.2% has a large distribution area of Pareto-optimal solutions and the frontier surface formed by the Pareto solution set is uniformly distributed, which meets the needs of design. Besides, with the significant increase of CMDG, the improvement of mechanical performance of MR damper is not obvious.


### 4.1.2. Analysis of Parametric Correlation

In order to analyze the uncertainty in MR damper parameters and to find which parameter will influence the performances of MR dampers most, we have done some research studies to obtain the correlation coefficient between each value of the optimal objective and the geometric dimension parameter in the Pareto-optimal solutions.The linear correlation coefficient is calculated by using the linear correlation analysis formula (29). The calculating results are listed in Table 4 and shown in Figure 12:(29)ρxy=n∑xiyi−∑xi∑yin∑xi2−∑xi2n∑yi2−∑yi2×100%,where ρxy is the linear correlation coefficient and xi and yi are the ith parameter of the optimal objective and the structural parameter, respectively.Table 4
Correlation coefficient matrix (%).
F
K
T
C
D
Pis
79.66
−96.13
21.70
1.25
D
Rod
−0.03
1.23
−4.60
−2.91
L
Coil
−13.87
−7.64
−54.55
−51.13
H
Coil
12.79
−5.28
47.52
5.37
t
Case
1.67
0.99
0.06
−0.57
w
1
6.68
4.15
94.53
−42.28Figure 12
Correlation coefficient between the optimal objective and the geometric dimension parameters. (a) The correlation betweenF and geometric dimensions. (b) The correlation between K and geometric dimensions. (c) The correlation between (T) and geometric dimensions. (d) The correlation between C and geometric dimensions.
(a)
(b)
(c)
(d)As shown in Figure12, the piston diameter DPis makes a great influence on the performances of the damping force and the dynamic range in all the geometric dimension parameters. And the flank piston thickness W1 has higher influence on the performances of the response time, and the CMDG also performs apparently than other parameters. Besides, the DRod and tCase are useless to improve the optimization objective of the MR dampers, as shown in Table 4. It may be partly because of the small values of geometric dimension in the rod diameter and the cylinder body thickness compared to other dimensions.


### 4.1.3. Scheme Comparison

To compare the influence of CMDG difference, we filter out solutions where the performance indexes of the damping forceF, the dynamic range K, and the response time T are better than those of SG-MRD60 scheme from the Pareto-optimal solutions set by using Excel. According to the optimization results of CMDG, the solution whose CMDG is approximately 0%, 15%, and 30% are selected and named as scheme A, B, and C, respectively. All schemes including SG-MRD60 scheme are shown in Table 5.Table 5
Optimization results of partial geometry dimensions.
Parameter
SG-MRD60
Scheme A
Scheme B
Scheme C
F (kN)
5.03
5.14
5.05
5.27
K
32.24
33.69
43.03
41.73
T (ms)
215
190
182
166
C (%)
15.3
0.6
15.1
29.9
D
Pis (mm)
57.00
56.99
56.72
56.73
D
Rod (mm)
18.00
18.37
18.15
15.26
L
Coil (mm)
40.00
39.34
37.55
37.42
H
Coil (mm)
15.00
14.38
14.10
14.72
t
Case (mm)
7.00
8.08
6.92
7.79
W
1 (mm)
10.00
9.44
9.22
8.00
N
108
134
136
121According to Table5, compared with the SG-MRD60 scheme, the damping force, the dynamic range, and the response time in each scheme are increased by up to 2.2%, 4.5%, and 11.6%, respectively.


## 4.1.1. Results of Algorithm Calculation

The NSGA-III algorithm generates a total of 5986 Pareto-optimal solutions after 30000 iterations. The variables of geometric dimension corresponding to the Pareto-optimal solutions are taken into (10), (14), (20), and (21) respectively, and then, the abnormal solutions in the course of preserving the operation precision are removed. With damping force less than 35 kN, 4844 sets of data points are retained. The distribution of 3D Pareto-optimal front is drawn and shown as Figure 9. Then, we divide the Pareto-optimal solutions into 3 layers according to the range of response time and CMDG from the 4844 sets of data points. And the distribution of points is shown in Figures 10 and 11.Figure 9
3D Pareto-optimal solutions.Figure 10
Layer of Pareto-optimal solutions for response time.Figure 11
Layer of Pareto-optimal solutions for CMDG.As shown in Figure9, the damping force is negatively correlated with the dynamic range, which is in accordance with results of similar research [13]. According to Figure 10, the Pareto-optimal solutions are concentrated in a region with low damping force and high dynamic range when response time is less than 160 ms. With the increase of response time, the Pareto-optimal solutions tend to move to the region with high damping force and low dynamic range. And the comprehensive mechanical performance which refers to damping force and dynamic range will be improved. As shown in Figure 11, the design scheme with CMDG less than 0.2% has a large distribution area of Pareto-optimal solutions and the frontier surface formed by the Pareto solution set is uniformly distributed, which meets the needs of design. Besides, with the significant increase of CMDG, the improvement of mechanical performance of MR damper is not obvious.


## 4.1.2. Analysis of Parametric Correlation

In order to analyze the uncertainty in MR damper parameters and to find which parameter will influence the performances of MR dampers most, we have done some research studies to obtain the correlation coefficient between each value of the optimal objective and the geometric dimension parameter in the Pareto-optimal solutions.The linear correlation coefficient is calculated by using the linear correlation analysis formula (29). The calculating results are listed in Table 4 and shown in Figure 12:(29)ρxy=n∑xiyi−∑xi∑yin∑xi2−∑xi2n∑yi2−∑yi2×100%,where ρxy is the linear correlation coefficient and xi and yi are the ith parameter of the optimal objective and the structural parameter, respectively.Table 4
Correlation coefficient matrix (%).
F
K
T
C
D
Pis
79.66
−96.13
21.70
1.25
D
Rod
−0.03
1.23
−4.60
−2.91
L
Coil
−13.87
−7.64
−54.55
−51.13
H
Coil
12.79
−5.28
47.52
5.37
t
Case
1.67
0.99
0.06
−0.57
w
1
6.68
4.15
94.53
−42.28Figure 12
Correlation coefficient between the optimal objective and the geometric dimension parameters. (a) The correlation betweenF and geometric dimensions. (b) The correlation between K and geometric dimensions. (c) The correlation between (T) and geometric dimensions. (d) The correlation between C and geometric dimensions.
(a)
(b)
(c)
(d)As shown in Figure12, the piston diameter DPis makes a great influence on the performances of the damping force and the dynamic range in all the geometric dimension parameters. And the flank piston thickness W1 has higher influence on the performances of the response time, and the CMDG also performs apparently than other parameters. Besides, the DRod and tCase are useless to improve the optimization objective of the MR dampers, as shown in Table 4. It may be partly because of the small values of geometric dimension in the rod diameter and the cylinder body thickness compared to other dimensions.


## 4.1.3. Scheme Comparison

To compare the influence of CMDG difference, we filter out solutions where the performance indexes of the damping forceF, the dynamic range K, and the response time T are better than those of SG-MRD60 scheme from the Pareto-optimal solutions set by using Excel. According to the optimization results of CMDG, the solution whose CMDG is approximately 0%, 15%, and 30% are selected and named as scheme A, B, and C, respectively. All schemes including SG-MRD60 scheme are shown in Table 5.Table 5
Optimization results of partial geometry dimensions.
Parameter
SG-MRD60
Scheme A
Scheme B
Scheme C
F (kN)
5.03
5.14
5.05
5.27
K
32.24
33.69
43.03
41.73
T (ms)
215
190
182
166
C (%)
15.3
0.6
15.1
29.9
D
Pis (mm)
57.00
56.99
56.72
56.73
D
Rod (mm)
18.00
18.37
18.15
15.26
L
Coil (mm)
40.00
39.34
37.55
37.42
H
Coil (mm)
15.00
14.38
14.10
14.72
t
Case (mm)
7.00
8.08
6.92
7.79
W
1 (mm)
10.00
9.44
9.22
8.00
N
108
134
136
121According to Table5, compared with the SG-MRD60 scheme, the damping force, the dynamic range, and the response time in each scheme are increased by up to 2.2%, 4.5%, and 11.6%, respectively.


## 4.2. Results of Simulation

In order to compare the optimization results with the simulation results of CMDG, we import the parameters of geometric dimension (see in Table5) into ANSYS finite element analysis software. The Nephogram magnetic flux density distribution for each scheme is obtained and shown in Figure 13. And the simulation results of magnetic flux density of the damping gap are listed in Table 6. Besides, the average magnetic flux density of the damping gap Bf−m¯ is given by(30)Bf−m¯=Bf×2+Bm3.Figure 13
Nephogram magnetic flux density distribution obtained by FEA. (a) SG-MRD60 scheme. (b) Scheme A. (c) Scheme B. (d) scheme C.
(a)
(b)
(c)
(d)Table 6
Magnetic flux density result obtained by FEA.
Parameter
SG-MRD60
Scheme A
Scheme B
Scheme C
B
f (T)
0.242
0.326
0.333
0.312
B
m (T)
0.286
0.340
0.299
0.242
B
f
−
m
¯ (T)
0.257
0.331
0.332
0.289
C
ANSYS (%)
15.38
4.11
11.37
28.93As can be seen in Figure13(a)–13(c), the original scheme from SG-MRD60, whose magnetic flux density of middle damping gap is significantly greater than that of in the flank damping gap, has poor performance in CMDG. The CMDG of schemes A and B are better than that of SG-MRD60. It is in accordance with the result in Table 6. Besides, the trend of CMDG in Table 6 is the same as that in Table 5. The average deviation between the optimization result and the simulation result is 0.4%, which reflects that the CMDG function can express the distribution effect of the magnetic flux density in damping gap well.As shown in Table6, the CMDG in scheme C is 28.93%, which is the largest of all schemes. It leads to the consequence that the simulation result of CMDG in scheme C is the worst as shown in Figure 13(d), which means CMDG with maximum value has the worst performance. The results show that the worse the CMDG performance is, the more unbalanced the magnetic flux density of flank damping gap and middle damping gap will be. It will affect the calculation result of the optimization objective function of damping force and make the optimization scheme invalid. Therefore, the introduction of CMDG into multiobjective optimization design for MR dampers is very important.


## 5. Conclusions

This study optimizes the geometric dimensions of the selected MR damper using NSGA-III algorithm in the PlatEMO platform. The objectives of the optimization are maximizing the damping force and the dynamic range and minimizing the response time and the CMDG of the damper. This study also proposes a CMDG formula and introduces the CMDG into the process of optimization design. The main conclusions can be listed as follows.Optimized by NSGA-III algorithm, the Pareto-optimal front shows that the damping force is negatively correlated with the dynamic range and the response time in the Pareto-optimal solutions. And the CMDG index is greatly influenced by the geometric dimensions of theLCoil and W1. The Pareto-optimal front also shows that the damping force and the dynamic range have no correlation with CMDG, which indicates that during the process of final design scheme selection from the Pareto-optimal solutions, the performance of CMDG does not need to be ensured by reducing other performances of the MR dampers.Compared with the SG-MRD60 scheme, the optimized design scheme achieves higher performance of the MR damper. The finite element analysis results of CMDG are close to the optimization results with an average deviation of 0.4%. Besides, the simulation results of CMDG decrease by 11.27%, which is superior to the simulation results using the SG-MRD60 scheme. It indicates that the introduction of the CMDG index can make the procedure of multiobjective optimization design more valuable for engineering application.In addition, the optimization results of actual damping force, dynamic range, and response time need to be further tested and verified after the MR damper is manufactured.



---

*Source: 7050356-2020-01-17.xml*