**3.4. Combination of modes**

Some of the applications of field responsive fluids take advantage of the combination of two modes for a greater strength and functionality. For instance, dampers can be constructed in three different modes [75]. In a general manner, shear mode exhibits Couette flow through the annular bypass, while a valve flow is characterized by Poiseiulle flow through the annular bypass. The combination of them often gives higher yield stress as compared to stress produced by individual operation modes. Kamath et al. [76] have shown in their analysis and testing of Bingham plastic behaviour that mixed (valve and shear) mode dashpot dampers have higher passive damping than flow mode dampers. The mixed mode damper has a secondary effect of viscous drag as a result of the motion of piston head, instead of relying on the pressure gradient developed by the piston head to push the fluid through the gap created by the fixed electrodes. Wereley and Pang [75] have developed nonlinear quasi-steady ER and MR damper models using idealized Bingham plastic shear flow mechanism to characterize the equivalent viscous damping constant of the dampers. Plug thickness is the strongest variable that contributed to the damper behaviour for both flow and mixed modes.

In another experimental study done by Kulkarni et al. [62], the performance ofhe combination of squeeze and shear modes of MRF in dynamic loading was investigated. Even though squeeze mode can produce the highest strength among all modes, the addition of squeeze mode to shear mode did not always give a better strength than the shear mode alone. However, Tang et al. [77] demonstrated that the yield shear tress can be significantly improved by compressing the MRF along the magnetic field direction before the shear process is performed.

### **4. Optimal design methodology of MRF-based mechanisms**

#### **4.1. Modeling of MRF based mechanisms**

It is well-known that modeling of the MRF based systems is a coupled analysis problem: electromagnetic analysis and fluid system analysis. The purpose of the modeling of an MRF based device is to find the relation between the applied electric power (usually the current applied to the coils) and the output mechanical power such as pressure drops for MR valves, damping force for MR damper, braking torque for MR brakes and transmitted torque for MR clutches. In order to deal with modeling of MRF based devices, firstly the magnetic circuit of the MRF based devices should be solved. In general, the magnetic circuit can be analyzed using the magnetic Kirchoff's law as follows:

$$\sum H\_k l\_k = N\_{turns} I \tag{14}$$

where *Hk* is the magnetic field intensity in the *kth* link of the circuit and *lk* is the effective length of that link. *Nturns* is the number of turns of the valve coil and *I* is the applied current. The magnetic flux conservation rule of the circuit is given by

$$
\Phi = \mathcal{B}\_k \mathcal{A}\_k \tag{15}
$$

where Φ is the magnetic flux of the circuit, *Ak* and *Bk* are the cross-sectional area and magnetic flux density of the *kth* link, respectively. It is noteworthy that the more links are used the more exact solution can be obtained. However, this increases computation load. At low magnetic field, the magnetic flux density, *Bk*, increases in proportion to the magnetic intensity *Hk* as follows:

$$B\_k = \mu\_0 \mu\_k H\_k \tag{16}$$

where μ*0* is the magnetic permeability of free space (μ*0*= 4π10-7*Tm/A*) and μ*<sup>k</sup>* is the relative permeability of the *kth* link material. As the magnetic field becomes large, its ability to polarize the magnetic material diminishes and the material is almost magnetically saturated. Generally, a nonlinear *B-H* curve is used to express the magnetic property of material. At low magnetic field, taking the linear relation (16) into consideration, the magnetic flux density and the field intensity of the *kth* link of magnetic circuit can be approximately calculated as follows:

Optimal Design Methodology of Magnetorheological Fluid Based Mechanisms 359

$$B\_k = \frac{\mu\_0 N\_{turns} I}{\frac{I\_k}{\mu\_k} + \sum\_{i=1, i \neq k}^n \frac{l\_i A\_k}{\mu\_i A\_i}}\tag{17}$$

$$H\_k = \frac{N\_{turns}I}{l\_k + \sum\_{i=1, i \neq k}^{n} \frac{\mu\_k A\_k}{\mu\_i A\_i} l\_i} \tag{18}$$

By assuming magnetic property of the structural materials of the MR devices is similar (μ*1*=μ*2*= …μ*n*=μ), the magnetic flux density and the field intensity across the active MRF volume can be approximately calculated as follows:

358 Smart Actuation and Sensing Systems – Recent Advances and Future Challenges

**4. Optimal design methodology of MRF-based mechanisms** 

**4.1. Modeling of MRF based mechanisms** 

can be analyzed using the magnetic Kirchoff's law as follows:

The magnetic flux conservation rule of the circuit is given by

*0* is the magnetic permeability of free space (

process is performed.

intensity *Hk* as follows:

calculated as follows:

where μ

In another experimental study done by Kulkarni et al. [62], the performance ofhe combination of squeeze and shear modes of MRF in dynamic loading was investigated. Even though squeeze mode can produce the highest strength among all modes, the addition of squeeze mode to shear mode did not always give a better strength than the shear mode alone. However, Tang et al. [77] demonstrated that the yield shear tress can be significantly improved by compressing the MRF along the magnetic field direction before the shear

It is well-known that modeling of the MRF based systems is a coupled analysis problem: electromagnetic analysis and fluid system analysis. The purpose of the modeling of an MRF based device is to find the relation between the applied electric power (usually the current applied to the coils) and the output mechanical power such as pressure drops for MR valves, damping force for MR damper, braking torque for MR brakes and transmitted torque for MR clutches. In order to deal with modeling of MRF based devices, firstly the magnetic circuit of the MRF based devices should be solved. In general, the magnetic circuit

where *Hk* is the magnetic field intensity in the *kth* link of the circuit and *lk* is the effective length of that link. *Nturns* is the number of turns of the valve coil and *I* is the applied current.

where Φ is the magnetic flux of the circuit, *Ak* and *Bk* are the cross-sectional area and magnetic flux density of the *kth* link, respectively. It is noteworthy that the more links are used the more exact solution can be obtained. However, this increases computation load. At low magnetic field, the magnetic flux density, *Bk*, increases in proportion to the magnetic

> *k kk* <sup>0</sup> *B H* = μ μ

permeability of the *kth* link material. As the magnetic field becomes large, its ability to polarize the magnetic material diminishes and the material is almost magnetically saturated. Generally, a nonlinear *B-H* curve is used to express the magnetic property of material. At low magnetic field, taking the linear relation (16) into consideration, the magnetic flux density and the field intensity of the *kth* link of magnetic circuit can be approximately

μ

*Hl N I k k turns* <sup>=</sup> (14)

*k k* Φ = *B A* (15)

*0*= 4π10-7*Tm/A*) and

(16)

μ

*<sup>k</sup>* is the relative

$$B\_{mr} = \frac{\mu\_0 N\_{turns} I}{\frac{l\_{mr}}{\mu\_{mr}} + \frac{1}{\mu} \sum\_{i} \frac{l\_i A\_{MR}}{A\_i}} \tag{19}$$

$$H\_{mr} = \frac{N\_{turns}I}{l\_{mr} + \frac{\mu\_{mr}A\_{mr}}{\mu} \sum\_{i} \frac{l\_i}{A\_i}}\tag{20}$$

where μ*mr* and μ are the relative permeability of MRF and the structural materials of the MR devices, respectively. It is noted that the permeability of the MRF is much smaller than that of the valve core material, therefore from Eq. (20) the magnetic field intensity of the MRF link can be approximated by

$$H\_{mr} = \mathcal{N}\_{turns} I \not\!\!/ \mathcal{l}\_{mr} \tag{21}$$

The inductive time constant (*Tin*) and the power consumption (*N*) of the MRF based devices can be calculated as follows:

$$T\_{in} = L\_{in} \not\subset \mathbb{R}\_w \tag{22}$$

$$N = I^2 R\_w \tag{23}$$

where *Lin* is the inductance of the coil given by *Lin=Nturns*Φ*/I*, *Rw* is the resistance of the coil wire which can be approximately calculated by

$$R\_w = L\_w r\_w = N\_{turns} \pi \overline{d}\_c \frac{r}{A\_w} \tag{24}$$

In the above, *Lw* is the length of the coil wire, *rw* is the resistance per unit length of the coil wire, *dc* is the average diameter of the coil, *Aw* is the cross sectional area of the coil wire, *r* is the resistivity of the coil wire, *r* = 0.01726E-6Ω*m* for copper wire, *Nturns* is the number of coil turns which can be approximated by *Nturns=Ac/Aw*, and *Ac* is the cross-sectional area of the coil.

In some applications, an electromagnet is used in combination with a permanent magnet to control the rheological properties of the MRF as shown in Figure 7. In this case, the permanent magnet is used to shift the off-state (no current in the coil) viscosity of the MRF to a selected value and the electromagnet is used to control the viscosity variations around this value. A frequent situation is that where the magnetic circuit is designed in such a way that the MRF viscosity is maximum when no current flows through the coil. This is particularly useful when the device based on such magnetic circuit has to be blocked the major part of its operation time (such as in release mechanisms, for instance). The magnetic intensity across the active volume of the MRF is determined by [78]

**Figure 7.** Magnetic coil with MRF filled gap and permanent magnet: a) magnetic circuit b) electric equivalence

$$\mathbf{H}\_{mr} = \frac{(\mu\_r/l\_m)\mathbf{N}\_{turn}\mathbf{I} - \mathbf{B}\_r/\mu\_0}{\mu\_r l\_{mr}/l\_m + \mu\_{mr}A\_{mr}/A\_m} \tag{25}$$

where *lm* and *Am* are the length and cross-sectional area of the permanent magnet, *Br* and μ*r* are remanent flux density and relative permeability of the magnet (μ*<sup>r</sup>*≅1). In order to cancel the flux inside the MRF, *NI* has to be equal to the magnetomotive force, *Fm=Brlm/*μ*<sup>0</sup>*, which is highly influenced by the magnet length. It should however be noticed that such a magnetic circuit will not be used in practice since it might lead to demagnetization of the permanent magnet. A solution to this problem is to include a secondary path inside the magnetic circuit as shown in Figure 8 [78]

**Figure 8.** Magnetic coil with MRF filled gap and permanent magnet with secondary path: a) magnetic circuit b) electric equivalence

The electromagnet will thus not be used to completely cancel the flux produced by the permanent magnet but will only redirect it to the secondary path. This secondary path comprises a higher reluctance air gap in order to concentrate the major part of the flux generated by the permanent magnet in the primary path (comprising the MRF gap) when no current is flowing through the coil. In this case, the magnetic intensity across the active volume of the MRF is determined by [78]

360 Smart Actuation and Sensing Systems – Recent Advances and Future Challenges

intensity across the active volume of the MRF is determined by [78]

*mr*

are remanent flux density and relative permeability of the magnet (

μ

the flux inside the MRF, *NI* has to be equal to the magnetomotive force, *Fm=Brlm/*

equivalence

as shown in Figure 8 [78]

circuit b) electric equivalence

In some applications, an electromagnet is used in combination with a permanent magnet to control the rheological properties of the MRF as shown in Figure 7. In this case, the permanent magnet is used to shift the off-state (no current in the coil) viscosity of the MRF to a selected value and the electromagnet is used to control the viscosity variations around this value. A frequent situation is that where the magnetic circuit is designed in such a way that the MRF viscosity is maximum when no current flows through the coil. This is particularly useful when the device based on such magnetic circuit has to be blocked the major part of its operation time (such as in release mechanisms, for instance). The magnetic

**Figure 7.** Magnetic coil with MRF filled gap and permanent magnet: a) magnetic circuit b) electric

*lN IB <sup>H</sup>*

<sup>−</sup> <sup>=</sup> <sup>+</sup> μ

where *lm* and *Am* are the length and cross-sectional area of the permanent magnet, *Br* and

highly influenced by the magnet length. It should however be noticed that such a magnetic circuit will not be used in practice since it might lead to demagnetization of the permanent magnet. A solution to this problem is to include a secondary path inside the magnetic circuit

**Figure 8.** Magnetic coil with MRF filled gap and permanent magnet with secondary path: a) magnetic

<sup>0</sup> () / / / *r m turn r*

*r mr m mr mr m*

 μ

(25)

*<sup>r</sup>*≅1). In order to cancel

μ

μ

μ*r*

*<sup>0</sup>*, which is

*ll AA*

 μ

$$\mathbf{H}\_{mr} = \frac{(\mu\_r/\mathbf{l}\_m + \mathbf{A}\_a \mu\_a \,/\, \mathbf{A}\_m \mathbf{l}\_a) \mathbf{N}\_{\rm tur} \mathbf{I} - \mathbf{B}\_r \,/\, \mu\_0}{\mu\_r \mathbf{l}\_{mr} \,/\, \mathbf{l}\_m + \mathbf{A}\_a \mu\_a \mathbf{l}\_{mr} \,/\, \mathbf{A}\_m \mathbf{l}\_a + \mu\_{mr} \mathbf{A}\_{mr} \,/\, \mathbf{A}\_m} \tag{26}$$

It is interesting to note that, if *ga*→∞, we come back to Eq. (25). In order to cancel the magnetic flux inside the MRF gap, we need:

$$N\_{turn}I = \frac{B\_r l\_m}{\mu\_r \mu\_{m0} \left[ (\frac{A\_a}{A\_m})(\frac{l\_m}{l\_a})(\frac{\mu\_a}{A\_r}) + 1 \right]} \tag{27}$$

This value may seem smaller than what was obtained in the previous case; however, to obtain the same magnetic field inside the MRF gap, the magnet has to be more powerful since it has to compensate for the loss of magnetic flux in the secondary circuit.

In the above, magnetic circuit of the MRF based devices is solved based on the approximation of the analytical analysis. This approach can only used in case of simple geometry. In case of complex geometry or several coils are used, the approach becomes very complicated. Therefore, practically, the magnetic circuit of the MRF based devices is solved by finite element method (FEM). Once the magnetic solution is obtained, the magnetic intensity and magnetic flux density across the active MRF volume can be calculated. The rheological properties of MRF in the active volume are then determined based on the behaviour characteristics of the employed MRF. The behavior characteristics of MRF are usually obtained from experimental results with a curve-fitting algorithm. The most important parameter of MRF is the field-dependent yield stress. There have been several approximate functions have been used to express the dependence of the induced yield stress of MRF on the applied magnetic field. The two most widely used functions are the exponential function and the polynomial function. The former can well expressed the saturation of MRF yield stress as a function of the applied magnetic intensity. However, it exhibits large error at the small value of the applied magnetic intensity. In general, the approximate exponential function of induced yield stress is expressed as following

$$
\pi\_y(H) = \pi\_0 + \alpha H^\beta \tag{28}
$$

where τ*<sup>y</sup>*(*H*) is the induced yield stress of MRF as a function of the applied magnetic intensity (*H*), α and β are the curve parameters determined from experimental results using a curve-fitting algorithm, and τ*<sup>0</sup>* is the zero-field yield stress of the MRF.

The latter, the approximate polynomial function, can well predict the MRF yield stress at small value of the applied magnetic intensity. The higher order of the polynomial is the more accurate value of the yield stress can be predicted. In practice, the third order polynomial is often used. However, the polynomial function can not express the saturation of the induced yield stress. Therefore, a saturation condition should be added. The 3rd order approximate polynomial function of MRF yield stress can be expressed by

$$
\pi\_y(H) = \pi\_0 + c\_1 H + c\_2 H^2 + c\_3 H^3 \tag{29}
$$

where *c1*, *c2*, and *c3* are the curve parameters determined from experimental results using a curve-fitting algorithm, and τ*<sup>0</sup>* is the zero-field yield stress of the MRF.

In many researches, other characteristics of MRF such as the post yield viscosity (η) in Bingham model, the consistency parameter (*K*) and the fluid behavior index (*m*) are determined from experimental results on rheological properties of MRF and assumed to be independent of the applied magnetic. However, in practice, these parameters are slightly affected by the applied magnetic field. In order to take this into account, Zubieta et al. [44] have proposed a field-dependent plastic model for MRF based on original Bingham plastic and Herschel-Bulkley plastic models as mention in Section 2.5. Once obtaining the yield stress and other rheological parameters of the MRF, the output mechanical power such as pressure drops, damping force for MR damper, braking torque and transmitted torque can be determined on governing equations of the MRF based devices.

#### **4.2. Optimization problems in design of MRF based devices**

As aforementioned modeling of MRF based systems is a coupled analysis problem. Therefore, output mechanical power of these systems depends not only on their mechanics behaviors but also on their magnetic circuits. It is obvious that in order to improve performance of the MRF based systems, the optimal design should be taken into account. Generally, the objective of the optimal design is to find significant geometric dimensions of the MRF based devices that maximize an objective function considering typical characteristics such as pressure drop, damping force, dynamic range, braking torque, transmitted force, mass, time response constant and power consumption. Some constraints such as available space, allowable operating temperature, uncontrollable torque etc. may be also considered in the optimal design. There have been several researches focusing on optimal design of MRF devices. Rosenfield and Wereley [48] proposed analytical optimization design method for MR valves and dampers based on the assumption of constant magnetic flux density throughout the magnetic circuit to ensure that one region of the magnetic circuit does not saturate prematurely and cause a bottleneck problem. Nguyen et al. [79] proposed a FEM based optimal design of MR valves (single-coil, two-coil, threecoil and radial-annular types) constrained in a specified volume. This work considered the effects of all geometric variables of MR valves by minimizing the valve ratio calculated from the FE analysis. Later on Nguyen el al. [80] have developed an optimization procedure based on the finite element method in order to find the optimal geometry of MR valves constrained to a specific volume, satisfying a required pressure drop with minimal power consumption. The time response of the valves was also taken into account by considering the inductive time constant as a state variable. The optimization results showed the significance of the optimal design of the MR valves in order to minimize the power consumption. It was also shown that the wire diameter does not significantly affect the optimization solution and can be neglected. The optimal design of MR damper was also performed by Nguyen et al [81], in which the objective function was proposed by a linear combination of the ratios of the damping force, dynamic range and the inductive time constant and their reference values using corresponding weighting factors. Recently, there have been several researches on the optimal design of MR brakes and clutches. Park et al. [82] have performed multidisciplinary design optimization of an automotive MR brake, in which a multi-objective function considering both braking torque and mass of the brake was considered. Nguyen et al. [45] have performed a thorough research on optimal design of MR brake for middle-sized vehicle considering the available space, mass, braking torque and steady heat generated by a zero-field friction torque of the MR brake on cruising at a speed of 100*km/h.* Furthermore, different configurations of MR brake and different types of MRF are taken into account in that research. More recently, Nguyen et al. [83] have performed the optimal design of common types of MR brakes such as disc-type, drum-type, inverted drum-type, single-coil hybrid-types, inverted single-coil hybrid-types, two-coil hybridtypes, inverted two-coil hybrid-types and T-type. The objective of the optimization was to maximize the braking torque while torque ratio (the ratio of maximum braking torque and the zero field friction torque) is constrained not to exceed a certain value. Based on the optimal solutions, the advices on optimal selection of MR brakes type were addressed. It was showed that the guide on optimal selection of MR brake types can be applied for different types of MRF and different constrains of torque ratio.

362 Smart Actuation and Sensing Systems – Recent Advances and Future Challenges

τ

be determined on governing equations of the MRF based devices.

**4.2. Optimization problems in design of MRF based devices** 

τ

curve-fitting algorithm, and

approximate polynomial function of MRF yield stress can be expressed by

 τ

The latter, the approximate polynomial function, can well predict the MRF yield stress at small value of the applied magnetic intensity. The higher order of the polynomial is the more accurate value of the yield stress can be predicted. In practice, the third order polynomial is often used. However, the polynomial function can not express the saturation of the induced yield stress. Therefore, a saturation condition should be added. The 3rd order

01 2 3 ( )

where *c1*, *c2*, and *c3* are the curve parameters determined from experimental results using a

Bingham model, the consistency parameter (*K*) and the fluid behavior index (*m*) are determined from experimental results on rheological properties of MRF and assumed to be independent of the applied magnetic. However, in practice, these parameters are slightly affected by the applied magnetic field. In order to take this into account, Zubieta et al. [44] have proposed a field-dependent plastic model for MRF based on original Bingham plastic and Herschel-Bulkley plastic models as mention in Section 2.5. Once obtaining the yield stress and other rheological parameters of the MRF, the output mechanical power such as pressure drops, damping force for MR damper, braking torque and transmitted torque can

As aforementioned modeling of MRF based systems is a coupled analysis problem. Therefore, output mechanical power of these systems depends not only on their mechanics behaviors but also on their magnetic circuits. It is obvious that in order to improve performance of the MRF based systems, the optimal design should be taken into account. Generally, the objective of the optimal design is to find significant geometric dimensions of the MRF based devices that maximize an objective function considering typical characteristics such as pressure drop, damping force, dynamic range, braking torque, transmitted force, mass, time response constant and power consumption. Some constraints such as available space, allowable operating temperature, uncontrollable torque etc. may be also considered in the optimal design. There have been several researches focusing on optimal design of MRF devices. Rosenfield and Wereley [48] proposed analytical optimization design method for MR valves and dampers based on the assumption of constant magnetic flux density throughout the magnetic circuit to ensure that one region of the magnetic circuit does not saturate prematurely and cause a bottleneck problem. Nguyen et al. [79] proposed a FEM based optimal design of MR valves (single-coil, two-coil, threecoil and radial-annular types) constrained in a specified volume. This work considered the effects of all geometric variables of MR valves by minimizing the valve ratio calculated from the FE analysis. Later on Nguyen el al. [80] have developed an optimization procedure based on the finite element method in order to find the optimal geometry of MR valves

In many researches, other characteristics of MRF such as the post yield viscosity (

*<sup>0</sup>* is the zero-field yield stress of the MRF.

2 3

*<sup>y</sup> H cH cH cH* =+ + + (29)

η) in

### **4.3. Optimal design of MRF devices based on finite element analysis**

As abovementioned, the magnetic circuit of the MRF based devices can be solved by an approximation of analytical solution or by FEM. Therefore, the optimal design of these devices can be performed based on either the analytical analysis or finite element analysis (FEA). The former is used only for simple devices such as single coil MR damper [84]. In this section, the optimal design of MRF devices based on FEA is introduced. First of all, an objective function should be proposed depending on the purpose of the optimal design and the application of the devices. It is noted that in the optimization problem the objective function is always minimized. Therefore, if the purpose of the optimization is to maximize a performance function of the devices, that function should be transform to an equivalent objective function. The equivalent objective function is the function that when it is minimized, the corresponding performance function is maximized. After the objective function is constructed, the design parameters of the optimization problem should be identified. In addition, the constraints of the optimization problem should be determined if there any. In the next step, an algorithm to obtain the optimal solution should be chosen. It is well-known that there have been numerous methods to find the optimal solution of an

optimization problem. They may be non-derivative, first-derivative or second-derivative methods. The non-derivative methods that do not require any derivative of the function are not usually applied to MRF based systems. Although they are generally easy to implement, their convergence properties are rather poor. They may work well in special cases when the function is quite random in character or the variables are essentially uncorrelated. Some typical non-derivative algorithms are the Simplex, Genetic Algorithms and Neural Networks. The second-derivative optimization methods are characterized by fast convergence and affine invariance. However, they require second derivatives and the solution of linear equation can be too expensive for large scale applications. The most popular optimization method, which is widely used in optimal design of MRF based devices, is the first order (derivative) method. Although the convergence rate of the first derivative method is somewhat slower than that of the second-derivative one, the first derivative method is still preferred in many applications because of its inexpensive cost for computation and programming. A typical first derivative optimization algorithm is the conjugate gradient method. The flow chart in Figure 9 shows how to find the optimal solution of MRF based devices based on ANSYS finite element software using the first order method. The procedures from the flow chart can be easily extended to other finite element software.

First of all, initial value of the design variables (DV) should be decided. Computation time of the optimization process significantly depends on the initial value of the DVs. Therefore, the initial value of design variables should be calculated based on a draft calculation or based on practical experience. Then, an analysis file for solving the magnetic circuit and calculating performance characteristics of the devices such as control energy, the inductive time constant, pressure drops, damping force, braking torque and transmitted torque is built. In ANSYS, the analysis file is built using parametric design language (APDL). It is noted that this analysis file can be created from a graphic user interface (GUI) model of ANSYS by using the *list>log file* submenu from the *file* menu of the ANSYS software. In the analysis file, the DVs must be used as symbolic variables and initial value is assigned to them. Generally, in order to calculate performance characteristics of the devices, the magnetic flux density across the active volume of MRF should be calculated. The magnetic flux density (*B*) and magnetic intensity (*H*) are not constant along the MRF duct, so an average should be used. The average magnetic flux density and intensity across the MR ducts was calculated by integrating flux density along a predefined path, then divided by the path length [79, 80]. In order to calculate the inductive time constant, firstly the magnetic flux is determined as follows:

$$\Phi = 2\pi R\_d \int\_{L\_p} B(s)ds\tag{30}$$

where *B(s)* is the magnetic flux density at each nodal point on the path, *s* is a dummy variable for the integration. The integration was performed along the path length, *Lp*. It is noteworthy that geometric dimensions of the MRF devices change during the optimization process, so that the meshing size of FE model should be specified by the number of elements per line rather than element size.

**Figure 9.** Flow chart for optimal design of MRF based devices using FEM.

extended to other finite element software.

flux is determined as follows:

per line rather than element size.

optimization problem. They may be non-derivative, first-derivative or second-derivative methods. The non-derivative methods that do not require any derivative of the function are not usually applied to MRF based systems. Although they are generally easy to implement, their convergence properties are rather poor. They may work well in special cases when the function is quite random in character or the variables are essentially uncorrelated. Some typical non-derivative algorithms are the Simplex, Genetic Algorithms and Neural Networks. The second-derivative optimization methods are characterized by fast convergence and affine invariance. However, they require second derivatives and the solution of linear equation can be too expensive for large scale applications. The most popular optimization method, which is widely used in optimal design of MRF based devices, is the first order (derivative) method. Although the convergence rate of the first derivative method is somewhat slower than that of the second-derivative one, the first derivative method is still preferred in many applications because of its inexpensive cost for computation and programming. A typical first derivative optimization algorithm is the conjugate gradient method. The flow chart in Figure 9 shows how to find the optimal solution of MRF based devices based on ANSYS finite element software using the first order method. The procedures from the flow chart can be easily

First of all, initial value of the design variables (DV) should be decided. Computation time of the optimization process significantly depends on the initial value of the DVs. Therefore, the initial value of design variables should be calculated based on a draft calculation or based on practical experience. Then, an analysis file for solving the magnetic circuit and calculating performance characteristics of the devices such as control energy, the inductive time constant, pressure drops, damping force, braking torque and transmitted torque is built. In ANSYS, the analysis file is built using parametric design language (APDL). It is noted that this analysis file can be created from a graphic user interface (GUI) model of ANSYS by using the *list>log file* submenu from the *file* menu of the ANSYS software. In the analysis file, the DVs must be used as symbolic variables and initial value is assigned to them. Generally, in order to calculate performance characteristics of the devices, the magnetic flux density across the active volume of MRF should be calculated. The magnetic flux density (*B*) and magnetic intensity (*H*) are not constant along the MRF duct, so an average should be used. The average magnetic flux density and intensity across the MR ducts was calculated by integrating flux density along a predefined path, then divided by the path length [79, 80]. In order to calculate the inductive time constant, firstly the magnetic

> 2 () *P d L* Φ = *R B s ds* π

where *B(s)* is the magnetic flux density at each nodal point on the path, *s* is a dummy variable for the integration. The integration was performed along the path length, *Lp*. It is noteworthy that geometric dimensions of the MRF devices change during the optimization process, so that the meshing size of FE model should be specified by the number of elements

(30)

After the analysis file is prepared, the procedures to achieve optimal design parameters of the MRF devices using the first order method of ANSYS optimization tool are performed as shown in Figure 9. Starting with initial value of DVs, by executing the analysis file, the initial value of the performance characteristics of the devices such as control energy, the inductive time constant, pressure drops, damping force, braking torque and transmitted torque are obtained. The ANSYS optimization tool then transforms the constrained optimization problem to an unconstrained one via penalty functions. The dimensionless, unconstrained objective function is formulated as follows:

$$\text{O}\left(\mathbf{x}, \boldsymbol{\eta}\right) = \frac{\text{OBJ}}{\text{OBJ}\_0} + \sum\_{i=1}^{n} P\_{\mathbf{x}}(\mathbf{x}\_i) + \eta \sum\_{i=1}^{m} P\_{\mathbf{g}}(\mathbf{g}\_i) \tag{31}$$

where *OBJ0* is the reference objective function value that is selected from the current group of design sets, *q* is the response surface parameter which controls constraint satisfaction. *Px* is the exterior penalty function applied to the design variable *x*. *Pg* is extended-interior penalty function applied to state variable (the constraint) *g*. For the initial iteration (*j = 0*), the search direction of DVs is assumed to be the negative of the gradient of the unconstrained objective function. Thus, the direction vector is calculated by

$$d^{(0)} = -\nabla Q(x^{(0)}, 1) \tag{32}$$

The values of DVs in next iteration (*j+1*) is obtained from the following equation,

$$\mathbf{x}^{(j+1)} = \mathbf{x}^{(j)} + s\_j d^{(j)} \tag{33}$$

where the line search parameter *sj* is calculated by using a combination of a golden-section algorithm and a local quadratic fitting technique. The analysis file is then executed with the new values of DVs and the convergence of the objective function, *OBJ,* is checked. If the convergence occurs, the values of DVs at the *jth* iteration are the optimum. If not, the subsequent iterations will be performed. In the subsequent iterations, the procedures are similar to those of the initial iteration with the exception of the direction vectors which are calculated according to Polak-Ribiere recursion formula as follows:

$$d^{(j)} = -\nabla \mathbb{Q}(\mathbf{x}^{(j)}, q\_k) + r\_{j-1} d^{(j-1)} \tag{34}$$

$$r\_{j-1} = \frac{\left[\nabla\mathbb{Q}(\mathbf{x}^{(j)}, \boldsymbol{\eta}) - \nabla\mathbb{Q}(\mathbf{x}^{(j-1)}, \boldsymbol{\eta})\right]^T \nabla\mathbb{Q}(\mathbf{x}^{(j)}, \boldsymbol{\eta})}{\left|\nabla\mathbb{Q}(\mathbf{x}^{(j-1)}, \boldsymbol{\eta})\right|^2} \tag{35}$$

Thus, each iteration is composed of a number of sub-iterations that include search direction and gradient computations.

It is noted that ANSYS software supports optimal design problems by integrating an optimization tool. Therefore, in most cases the optimal solution of the MRF based devices can be solved directly by the ANSYS software without interfacing with any programming software. In order to use the ANSYS optimization tool, it is necessary to set up optimization parameters. To do this, firstly the analysis file should be manually executed once to load all parameters in the analysis file into software buffer memory. After that, from the *Design Opt* menu, we specify the analysis file which will be used during optimization process, the DVs with their limits and tolerances, the state variables (if there are any) with limits and tolerances, the objective function with a convergence criteria, the method for solving the optimal solution, and the optimal output control option if necessary. In some cases, it is expected to employ some advanced optimization algorithms such as Genetic Algorithms, Neural Network, or user defined algorithms, the interfacing between the ANSYS and other software to perform the optimization such as Matlab, FORTRAN, C languages is required.
