Modeling the Response of Magnetorheological Fluid Dampers under Seismic Conditions

Magnetorheological (MR) fluid is a smart material fabricated by mixing magnetic-responsive particles with non-magnetic-responsive carrier fluids. MR fluid dampers are able to provide rapid and reversible changes to their damping coefficient. To optimize the efficiency and effectiveness of such devices, a computational model is developed and presented where the flow field is simulated using the computational fluid dynamics approach, coupled with the magnetohydrodynamics model. Three different inlet pressure profiles were designed to replicate real loading conditions are examined, namely a constant pressure, a sinusoidal pressure profile, and a pressure profile mimicking the 1994 Northbridge earthquake. When the MR fluid damper was in its off-state, a linear pressure drop between the inlet and the outlet was observed. When a uniform perpendicular external magnetic field was applied to the annular orifice of the MR damper, a significantly larger pressure drop was observed across the annular orifice for all three inlet pressure profiles. It was shown that the fluid velocity within the magnetized annular orifice decreased proportionally with respect to the strength of the applied magnetic field until saturation was reached. Therefore, it was clearly demonstrated that the present model was capable of accurately capturing the damping characteristics of MR fluid dampers.


Introduction
Advancements in nanotechnology have allowed for the fabrication of micrometer-and nanometer-sized magnetic-responsive particles, and when mixed with non-magnetic-responsive carrier fluids and at least one type of thiocarbamate and/or thiophosphrous additives, a type of smart fluid material is developed, known as the magnetorheological (MR) fluid [1]. Under the application of external magnetic fields, magnetic-responsive particles interact with each other through dipole moments to form chain-like structures along the direction of the applied magnetic field [2]. These chain structures contribute to resisting shear stress in the perpendicular direction to the magnetic field. Therefore, by varying the magnitude of external magnetic fields, MR fluids can exhibit a rapid, reversible, and tuneable transition between a liquid state (off-state) and a more viscous semisolid state (on-state) [3].
The initial development of MR fluid devices can be traced back to the late 1940s by Rabinow [4]. Carbonyl iron particles obtained from the thermal decomposition of iron pentacarbonyl are the most 2 of 16 commonly used magnetic-responsive particles in MR fluids due to their large saturation magnetization. Typical base fluids used in MR fluids include mineral and silicone oils, synthetic hydrocarbons, polyethers, polyesters, and water [5].
Bell et al. [6] studied the influences of particle shape on the properties of MR fluids. Experimental results showed that the off-state viscosity of MR fluids containing only wires was substantially higher than MR fluids containing only spherical particles. The effects of particle size on the properties of MR fluids were examined in References [7,8]. It was shown that larger magnetic-responsive particles had a larger yield stress compared to particles with smaller sizes. However, MR fluids with smaller particle sizes showed higher stability against sedimentation. Lemaire et al. [7] explained these effects as the results of fluctuations of particle positions in the particle chain. Figure 1 shows the three basic modes of operations of MR fluid devices [9]. The valve mode, or the pressure-driven flow mode, is often observed in dampers, shock absorbers, and servo-valves [5]. The shear mode is often utilized in clutches and brakes. The squeeze mode has not been as thoroughly studied as the previous two operating modes but can be used as low amplitude vibration dampers [9].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 17 commonly used magnetic-responsive particles in MR fluids due to their large saturation magnetization. Typical base fluids used in MR fluids include mineral and silicone oils, synthetic hydrocarbons, polyethers, polyesters, and water [5]. Bell et al. [6] studied the influences of particle shape on the properties of MR fluids. Experimental results showed that the off-state viscosity of MR fluids containing only wires was substantially higher than MR fluids containing only spherical particles. The effects of particle size on the properties of MR fluids were examined in References [7,8]. It was shown that larger magnetic-responsive particles had a larger yield stress compared to particles with smaller sizes. However, MR fluids with smaller particle sizes showed higher stability against sedimentation. Lemaire et al. [7] explained these effects as the results of fluctuations of particle positions in the particle chain. Figure 1 shows the three basic modes of operations of MR fluid devices [9]. The valve mode, or the pressure-driven flow mode, is often observed in dampers, shock absorbers, and servo-valves [5]. The shear mode is often utilized in clutches and brakes. The squeeze mode has not been as thoroughly studied as the previous two operating modes but can be used as low amplitude vibration dampers [9]. Among these applications, MR fluid dampers have received the most extensive research interest. MR fluid dampers, first introduced in 1992 by Carlson and Chrzan [10], utilize the controllable properties of MR fluids, including viscosity and shear stress, while only requiring a low power input. The peak power requirement in some commercially available dampers is less than 10 W, where the power from a small camera battery would allow for more than one hour of continuous operation [11,12]. The schematic diagram of a typical MR fluid damper is shown in Figure 2. Compared to conventional dampers, the damping coefficient of MR fluid dampers can be continuously and quickly adjusted by modifying the strength of the applied external magnetic field. Hence, they can be effectively used as shock absorbers to increase passenger comfort and safety in automobile ground vehicles [14][15][16][17][18], railway vehicles [19][20][21][22][23], and aircraft [24][25][26][27]. They can also be built into various civil structures, such as buildings [28][29][30], bridges [31,32], and piping systems [33], Among these applications, MR fluid dampers have received the most extensive research interest. MR fluid dampers, first introduced in 1992 by Carlson and Chrzan [10], utilize the controllable properties of MR fluids, including viscosity and shear stress, while only requiring a low power input. The peak power requirement in some commercially available dampers is less than 10 W, where the power from a small camera battery would allow for more than one hour of continuous operation [11,12]. The schematic diagram of a typical MR fluid damper is shown in Figure 2.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 17 commonly used magnetic-responsive particles in MR fluids due to their large saturation magnetization. Typical base fluids used in MR fluids include mineral and silicone oils, synthetic hydrocarbons, polyethers, polyesters, and water [5]. Bell et al. [6] studied the influences of particle shape on the properties of MR fluids. Experimental results showed that the off-state viscosity of MR fluids containing only wires was substantially higher than MR fluids containing only spherical particles. The effects of particle size on the properties of MR fluids were examined in References [7,8]. It was shown that larger magnetic-responsive particles had a larger yield stress compared to particles with smaller sizes. However, MR fluids with smaller particle sizes showed higher stability against sedimentation. Lemaire et al. [7] explained these effects as the results of fluctuations of particle positions in the particle chain. Figure 1 shows the three basic modes of operations of MR fluid devices [9]. The valve mode, or the pressure-driven flow mode, is often observed in dampers, shock absorbers, and servo-valves [5]. The shear mode is often utilized in clutches and brakes. The squeeze mode has not been as thoroughly studied as the previous two operating modes but can be used as low amplitude vibration dampers [9]. Among these applications, MR fluid dampers have received the most extensive research interest. MR fluid dampers, first introduced in 1992 by Carlson and Chrzan [10], utilize the controllable properties of MR fluids, including viscosity and shear stress, while only requiring a low power input. The peak power requirement in some commercially available dampers is less than 10 W, where the power from a small camera battery would allow for more than one hour of continuous operation [11,12]. The schematic diagram of a typical MR fluid damper is shown in Figure 2. Compared to conventional dampers, the damping coefficient of MR fluid dampers can be continuously and quickly adjusted by modifying the strength of the applied external magnetic field. Hence, they can be effectively used as shock absorbers to increase passenger comfort and safety in automobile ground vehicles [14][15][16][17][18], railway vehicles [19][20][21][22][23], and aircraft [24][25][26][27]. They can also be built into various civil structures, such as buildings [28][29][30], bridges [31,32], and piping systems [33], Compared to conventional dampers, the damping coefficient of MR fluid dampers can be continuously and quickly adjusted by modifying the strength of the applied external magnetic field. Hence, they can be effectively used as shock absorbers to increase passenger comfort and safety in Appl. Sci. 2019, 9, 4189 3 of 16 automobile ground vehicles [14][15][16][17][18], railway vehicles [19][20][21][22][23], and aircraft [24][25][26][27]. They can also be built into various civil structures, such as buildings [28][29][30], bridges [31,32], and piping systems [33], in order to minimize their responses to vibrations caused by strong winds and earthquakes. Studies on these systems were reviewed by Housner et al. [34], Spencer Jr. and Nagarajaiah [35], and Jung et al. [36].
In order to optimize the performance of MR fluid dampers in these applications, it is of fundamental importance to quantitatively examine the response (output) of MR fluid dampers under various input pressure profiles. In some modern studies [2,[37][38][39], MR fluids and other particulate mixtures are modelled using coupled computational fluid dynamics and discrete phase modelling approaches. However, this method is very computationally intensive, and can only be applied in modelling MR fluids in microscale configurations. The macroscale behavior of MR fluid dampers has been experimentally investigated in various studies [40][41][42]. However, effective approaches in numerically modelling such systems at the device level remain underdeveloped.
In this study, a macroscopic computational model was developed to simulate the response of MR fluid dampers under various input pressure profiles, including constant pressure, a sinusoidal pressure profile, and a seismic pressure based on the Northridge earthquake in 1994, utilizing the computational fluid dynamics and magnetohydrodynamics models. The effectiveness of damping was investigated by examining the static pressure at different locations along the annular orifice of the MR fluid damper under a range of applied external magnetic fields. This novel computational approach enabled the macroscale simulation of MR fluids by significantly reducing the computational resources required compared to other computational approaches, where only microscale simulations are possible. It can also significantly reduce the human and economic resources required in the optimization of MR fluid dampers when compared to experimental approaches.

Computational Fluid Dynamics
In computational fluid dynamics, the continuity and Navier-Stokes governing equations are used to simulate the flow behavior. Given the small size of the geometry, and hence the small hydrodynamic length scale, the flow can be assumed to be laminar. Also, assuming MR fluids are incompressible, the continuity equation of the fluid can be written as: where u is the vector field of the mixture velocity. The momentum of the fluid is governed by the Navier-Stokes equation: where t is the flow time; ν e f f and ρ e f f are the respective effective kinematic viscosity and density of the MR fluid; ∇p is the flow pressure gradient, which contributes to the internal source term; and S is the external source term due to the magnetic response of the MR fluid. The effective kinematic viscosity and density of the MR fluid is considered to be the weighted average between the carrier fluid (subscript f ) and the particle suspensions (subscript p). The effective density of the MR fluid is given as: where α is the volume fraction of the magnetic-responsive particles. The effective dynamic viscosity µ e f f of a two-phase mixture containing rigid spherical particles was first predicted by Einstein [43] as: where µ f is the dynamic viscosity of the continuous phase (carrier fluid), and k 1 is the Einstein shape factor. For rigid spheres, the Einstein shape factor is given by k 1 = 2.5. However, this approximation is only ideal for low particle volume fractions (α < 0.05). For higher particle volume fractions, Vand [44] proposed the following expression, which accounts for both mutual hydrodynamic interactions and particle collisions: where k 1 is the Einstein shape factor, k 2 is the shape factor for collision doublets, θ is the collision time constant, and B H is the hydrodynamic interaction constant. For spherical particles [44], k 1 = 2.5, k 2 = 3.175, θ = 4, and B H = 0.609. Once the effective dynamic viscosity µ e f f is determined, the effective kinematic viscosity ν e f f can be calculated via:

Magneto-Hydrodynamic Model
The magneto-hydrodynamic (MHD) model simulates the interaction between an electrically conducting fluid flow and the applied electromagnetic field. There are two main aspects of the coupling between the fluid flow field and the magnetic field: electromagnetic induction and Lorentz force. Electromagnetic induction is the effect of electric current induction as a result of the movement of conductive material in a magnetic field. Lorentz force is generated as a result between the interaction of the electric current and the magnetic field.
In MR fluid dampers, the macroscopic physics can be understood as occurring due to a Lorentz force. As illustrated in Figure 3, when the electromagnetic coil is in the on-state, a magnetic field B is generated in the horizontal direction. When the damper is actuated, the MR fluid experiences a changing magnetic field. According to Faraday's law of induction, an electromotive force (EMF) is induced and a current will flow around the annular orifice with current density j. The interaction between the current density and magnetic field creates a force that opposes the motion of the fluid [45]. This force is called a Lorentz Force and is equal to the cross product of the current density and magnetic field j × B. It is this interaction of the fluid and magnetic field that creates a damping force and will result in a higher pressure differential across the annular orifice. where is the dynamic viscosity of the continuous phase (carrier fluid), and k1 is the Einstein shape factor. For rigid spheres, the Einstein shape factor is given by k1 = 2.5. However, this approximation is only ideal for low particle volume fractions (α < 0.05).
For higher particle volume fractions, Vand [44] proposed the following expression, which accounts for both mutual hydrodynamic interactions and particle collisions: where k1 is the Einstein shape factor, is the shape factor for collision doublets, θ is the collision time constant, and BH is the hydrodynamic interaction constant. For spherical particles [44], k1 = 2.5, k2 = 3.175, θ = 4, and BH = 0.609.
Once the effective dynamic viscosity is determined, the effective kinematic viscosity can be calculated via:

Magneto-Hydrodynamic Model
The magneto-hydrodynamic (MHD) model simulates the interaction between an electrically conducting fluid flow and the applied electromagnetic field. There are two main aspects of the coupling between the fluid flow field and the magnetic field: electromagnetic induction and Lorentz force. Electromagnetic induction is the effect of electric current induction as a result of the movement of conductive material in a magnetic field. Lorentz force is generated as a result between the interaction of the electric current and the magnetic field.
In MR fluid dampers, the macroscopic physics can be understood as occurring due to a Lorentz force. As illustrated in Figure 3, when the electromagnetic coil is in the on-state, a magnetic field B is generated in the horizontal direction. When the damper is actuated, the MR fluid experiences a changing magnetic field. According to Faraday's law of induction, an electromotive force (EMF) is induced and a current will flow around the annular orifice with current density j. The interaction between the current density and magnetic field creates a force that opposes the motion of the fluid [45]. This force is called a Lorentz Force and is equal to the cross product of the current density and magnetic field . It is this interaction of the fluid and magnetic field that creates a damping force and will result in a higher pressure differential across the annular orifice.  In the MHD model, Maxwell's equations in the following form are used to describe the electromagnetic fields: where B is the magnetic field in T, and E is the electric field in V/m. In addition: where H and D are the induction fields for the magnetic and electric fields, respectively. q is the electric charge density in C/m 3 , and j is the electric current density vector in A/m 2 . For media with sufficiently high conductivity, the electric charge density q and the displacement current ∂D ∂t are considered negligible.
The induction fields H and D can be determined by: where µ is the magnetic permeability and ε is the electric permittivity. The electric current density vector j may be determined using two different approaches, namely the magnetic induction method and the electric potential method. The magnetic induction method is derived from Ohm's law and the Maxwell's equation, whereas the electric potential method solves the electric potential equation to determine the current density using Ohm's law. In this study, the electric potential method is used.
The current density is defined by Ohm's law using: where σ is the electrical conductivity of the medium. In a flow under an external magnetic field B, Ohm's law can be expressed as: where U is the velocity of the fluid flow. In the electric potential approach, the electric field E is given as: where ϕ is the scalar potential and A is the vector potential. When the magnetic field is static and much less than B 0 , the Ohm's law shown in Equation (14) can be expressed as: Due to the principle of conservation of electric charges in media with sufficiently high conductivity: Therefore, the electric potential equation can be written as: Finally, the following equation is derived based on the Ohm's law and Maxwell's equations: The coupling between the MHD model and the Navier-Stokes equation is achieved via the additional source term S in Equation (2) due to the Lorentz force, given by:

Material Properties
The MR fluid used in this computational model is assumed to be a mixture of water at 20 • C (ρ f = 999 kg/m 3 , µ f = 1.003 cP) and spherical iron particles (ρ f = 7870 kg/m 3 ) with a particle volume fraction α = 0.20. The effective density ρ e f f and viscosity µ e f f were determined using Equations (3) and (5), respectively. The magnetic permeability of the MR fluid was taken to be 4 H/m, and the electrical conductivity of the mixture was σ = 15,000 S/m.

Computational Domain
In order to simplify the underlying problem and reduce the computational power required, only MR fluid around the annular orifice region (see Figure 2) was modeled in this study. A two-dimensional axisymmetric configuration was implemented to reduce the computational effort. The computational domain is shown in Figure 4. Numerical results were extracted at six equally spaced locations throughout the domain marked as 1 -6 The shaded region between locations 3 and 4 represents the annular orifice where the MR fluid is affected by the external magnetic field.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 17 Finally, the following equation is derived based on the Ohm's law and Maxwell's equations: (19) The coupling between the MHD model and the Navier-Stokes equation is achieved via the additional source term in Equation (2) due to the Lorentz force, given by:

Material Properties
The MR fluid used in this computational model is assumed to be a mixture of water at 20 °C ( = 999 kg/m 3 , = 1.003 cP) and spherical iron particles ( = 7870 kg/m 3 ) with a particle volume fraction = 0.20. The effective density and viscosity were determined using Equations (3) and (5), respectively. The magnetic permeability of the MR fluid was taken to be 4 H/m, and the electrical conductivity of the mixture was = 15,000 S/m.

Computational Domain
In order to simplify the underlying problem and reduce the computational power required, only MR fluid around the annular orifice region (see Figure 2) was modeled in this study. A twodimensional axisymmetric configuration was implemented to reduce the computational effort. The computational domain is shown in Figure 4. Numerical results were extracted at six equally spaced locations throughout the domain marked as ①-⑥The shaded region between locations ③ and ④ represents the annular orifice where the MR fluid is affected by the external magnetic field.

Mesh
To numerically obtain the solutions for the aforementioned governing equations, a structured mesh consisting of rectangular elements was created for the two-dimensional computational domain since such elements are considered to be the most suitable due to its high conformity to the geometry and low aspect ratio.
A mesh sensitivity analysis was performed in order to arrive at a mesh-independent solution. To this end, based on the unmagnetized case, a mesh convergence study was conducted by gradually decreasing the element size and observing the changes in the maximum velocity in the fluid flow, as shown in Figure 5. It can be seen that the maximum velocity reached its asymptotic value (to four significant figures) once the mesh reached a 0.25 mm element size. For this element size, there were a total number of 48,000 mesh elements, and the maximum velocity within the domain was 0.3737 m/s.

Mesh
To numerically obtain the solutions for the aforementioned governing equations, a structured mesh consisting of rectangular elements was created for the two-dimensional computational domain since such elements are considered to be the most suitable due to its high conformity to the geometry and low aspect ratio.
A mesh sensitivity analysis was performed in order to arrive at a mesh-independent solution. To this end, based on the unmagnetized case, a mesh convergence study was conducted by gradually decreasing the element size and observing the changes in the maximum velocity in the fluid flow, as shown in Figure 5. It can be seen that the maximum velocity reached its asymptotic value (to four significant figures) once the mesh reached a 0.25 mm element size. For this element size, there were a total number of 48,000 mesh elements, and the maximum velocity within the domain was 0.3737 m/s.

Numerical Scheme
The finite volume method was used to discretize the fluid governing equations. Semi-Implicit Method for Pressure Linked Equations (SIMPLE) and Pressure-Implicit with Splitting of Operators (PISO) schemes were tested for the coupling between pressure and velocity, whilst first-order and second-order upwind methods were tested for spatial discretization. The underrelaxation factor for the pressure correction using the SIMPLE scheme was set to a value of 0.3. For the PISO scheme, the underrelaxation factor was set to unity, taking advantage of the solutions of the pressure correction in two stages. Owing to the simplicity in the geometry setup and mesh grids, the differences in numerical results between these two schemes were observed to be minimal (≈0.1%). For simplicity, the SIMPLE and second-order upwind schemes were adopted. Transient solutions were attained in time by marching the computations with a fixed time step of 10 −4 s. At each time step, the governing equations were iteratively solved and deemed to be converged until the Root-Mean-Square (RMS) residuals fell below 10 −6 .

Boundary Conditions
The boundary conditions applied in this study are shown in Figure 6. The top and bottom surfaces of the computational domain were defined as non-slip stationary walls to realistically represent the surfaces of the piston and housing of the MR fluid damper. The inlet and outlet of the domain were defined by specifying the pressure at these locations. The outlet pressure was set to be the gauge pressure for all the simulations conducted in this study, while three different cases of inlet pressure were investigated, as shown in Figure 7. The first case was a constant pressure of = 500 Pa. The second case involved the application of an oscillating sinusoidal pressure with respect to time: where 2 and f = 10 Hz.

Numerical Scheme
The finite volume method was used to discretize the fluid governing equations. Semi-Implicit Method for Pressure Linked Equations (SIMPLE) and Pressure-Implicit with Splitting of Operators (PISO) schemes were tested for the coupling between pressure and velocity, whilst first-order and second-order upwind methods were tested for spatial discretization. The underrelaxation factor for the pressure correction using the SIMPLE scheme was set to a value of 0.3. For the PISO scheme, the underrelaxation factor was set to unity, taking advantage of the solutions of the pressure correction in two stages. Owing to the simplicity in the geometry setup and mesh grids, the differences in numerical results between these two schemes were observed to be minimal (≈0.1%). For simplicity, the SIMPLE and second-order upwind schemes were adopted. Transient solutions were attained in time by marching the computations with a fixed time step of 10 −4 s. At each time step, the governing equations were iteratively solved and deemed to be converged until the Root-Mean-Square (RMS) residuals fell below 10 −6 .

Boundary Conditions
The boundary conditions applied in this study are shown in Figure 6. The top and bottom surfaces of the computational domain were defined as non-slip stationary walls to realistically represent the surfaces of the piston and housing of the MR fluid damper.

Numerical Scheme
The finite volume method was used to discretize the fluid governing equations. Semi-Implicit Method for Pressure Linked Equations (SIMPLE) and Pressure-Implicit with Splitting of Operators (PISO) schemes were tested for the coupling between pressure and velocity, whilst first-order and second-order upwind methods were tested for spatial discretization. The underrelaxation factor for the pressure correction using the SIMPLE scheme was set to a value of 0.3. For the PISO scheme, the underrelaxation factor was set to unity, taking advantage of the solutions of the pressure correction in two stages. Owing to the simplicity in the geometry setup and mesh grids, the differences in numerical results between these two schemes were observed to be minimal (≈0.1%). For simplicity, the SIMPLE and second-order upwind schemes were adopted. Transient solutions were attained in time by marching the computations with a fixed time step of 10 −4 s. At each time step, the governing equations were iteratively solved and deemed to be converged until the Root-Mean-Square (RMS) residuals fell below 10 −6 .

Boundary Conditions
The boundary conditions applied in this study are shown in Figure 6. The top and bottom surfaces of the computational domain were defined as non-slip stationary walls to realistically represent the surfaces of the piston and housing of the MR fluid damper. The inlet and outlet of the domain were defined by specifying the pressure at these locations. The outlet pressure was set to be the gauge pressure for all the simulations conducted in this study, while three different cases of inlet pressure were investigated, as shown in Figure 7. The first case was a constant pressure of = 500 Pa. The second case involved the application of an oscillating sinusoidal pressure with respect to time: where 2 and f = 10 Hz.  The inlet and outlet of the domain were defined by specifying the pressure at these locations. The outlet pressure was set to be the gauge pressure for all the simulations conducted in this study, while three different cases of inlet pressure were investigated, as shown in Figure 7. The first case was a constant pressure of P inlet = 500 Pa. The second case involved the application of an oscillating sinusoidal pressure with respect to time: P inlet (t) = 500 + 250 sin(ωt) (21) where ω = 2π f and f = 10 Hz.
The third case is based on the Northridge earthquake in 1994 [46]. According to the displacement measured during the earthquake, the inlet pressure was scaled and approximated using a piecewise function, as shown in Equation (22).

Results and Discussion
Simulation results for the constant, sinusoidal, and seismic response cases are presented. In the constant pressure case, the results for different magnetic field strengths are presented, while in the sinusoidal pressure and seismic response cases, the results for a single magnetic field strength of B = 1.0 T are presented. Figure 8 shows the static pressure measured across the annular orifice when a constant inlet pressure of 500 Pa and magnetic field strengths from 0 T to 1.5 T were applied. In the absence of a magnetic field, the resistance to shear flow experienced by the damper was due to the viscosity and the yield strength of the composite fluid. In this case, a linear pressure drop was observed between the inlet and the outlet of the computational domain, and the pressure drop across the annular orifice,

Constant Pressure
The third case is based on the Northridge earthquake in 1994 [46]. According to the displacement measured during the earthquake, the inlet pressure was scaled and approximated using a piecewise function, as shown in Equation (22).

Results and Discussion
Simulation results for the constant, sinusoidal, and seismic response cases are presented. In the constant pressure case, the results for different magnetic field strengths are presented, while in the sinusoidal pressure and seismic response cases, the results for a single magnetic field strength of B = 1.0 T are presented. Figure 8 shows the static pressure measured across the annular orifice when a constant inlet pressure of 500 Pa and magnetic field strengths from 0 T to 1.5 T were applied. In the absence of a magnetic field, the resistance to shear flow experienced by the damper was due to the viscosity and the yield strength of the composite fluid. In this case, a linear pressure drop was observed between the inlet and the outlet of the computational domain, and the pressure drop across the annular orifice, between locations 3 and 4 , was 53.82 Pa.

Results and Discussion
Simulation results for the constant, sinusoidal, and seismic response cases are presented. In the constant pressure case, the results for different magnetic field strengths are presented, while in the sinusoidal pressure and seismic response cases, the results for a single magnetic field strength of B = 1.0 T are presented. Figure 8 shows the static pressure measured across the annular orifice when a constant inlet pressure of 500 Pa and magnetic field strengths from 0 T to 1.5 T were applied. In the absence of a magnetic field, the resistance to shear flow experienced by the damper was due to the viscosity and the yield strength of the composite fluid. In this case, a linear pressure drop was observed between the inlet and the outlet of the computational domain, and the pressure drop across the annular orifice, between locations ③ and ④, was 53.82 Pa. In the cases where the annular orifice was magnetized, as the piston moved, a current was induced in the MR fluid. While in microscopic models, it has been shown that chain-like aggregates form and resist shear flow [2], in this macroscopic model, it was the resulting Lorentz force acting in opposition to the fluid flow that resulted in damping. This was demonstrated by the greater pressure differential between locations 3 and 4 . The pressure drop increased to 177.55 Pa, 333.63 Pa, and 409.7 Pa for the B = 0.5 T, 1.0 T, and 1.5 T cases, respectively.

Constant Pressure
When the magnetic field was increased from B = 0 T to B = 0.5 T, and again from B = 0.5 T to B = 1 T, there was a linear increase in the pressure differential across locations 3 and 4 . However, as the magnetic field was increased to B = 1.5 T, the pressure difference asymptotically increased. This resulted from the MR fluid approaching its magnetic saturation point. Increases in magnetic field strength above this point will only serve to marginally increase the damping coefficient of the damper.
In Figure 9, the unmagnetized case had a fluid velocity of 0.374 m/s at location 6 downstream from the inlet condition. In the magnetized cases, the fluid velocity was much lower downstream from the inlet with a deacceleration upon entering the annular orifice after location 3 and an acceleration upon exiting after location 4 . This deceleration and acceleration of the fluid in the magnetized cases were noticeably absent from the unmagnetized case. The deceleration of the fluid was caused by the Lorentz force generated after location 3 and the re-acceleration was caused by the fluid entering a low-pressure zone after location 4   In the cases where the annular orifice was magnetized, as the piston moved, a current was induced in the MR fluid. While in microscopic models, it has been shown that chain-like aggregates form and resist shear flow [2], in this macroscopic model, it was the resulting Lorentz force acting in opposition to the fluid flow that resulted in damping. This was demonstrated by the greater pressure differential between locations ③ and ④. The pressure drop increased to 177.55 Pa, 333.63 Pa, and 409.7 Pa for the B = 0.5 T, 1.0 T, and 1.5 T cases, respectively.
When the magnetic field was increased from B = 0 T to B = 0.5 T, and again from B = 0.5 T to B = 1 T, there was a linear increase in the pressure differential across locations ③ and ④. However, as the magnetic field was increased to B = 1.5 T, the pressure difference asymptotically increased. This resulted from the MR fluid approaching its magnetic saturation point. Increases in magnetic field strength above this point will only serve to marginally increase the damping coefficient of the damper.
In Figure 9, the unmagnetized case had a fluid velocity of 0.374 m/s at location ⑥ downstream from the inlet condition. In the magnetized cases, the fluid velocity was much lower downstream from the inlet with a deacceleration upon entering the annular orifice after location ③ and an acceleration upon exiting after location ④. This deceleration and acceleration of the fluid in the magnetized cases were noticeably absent from the unmagnetized case. The deceleration of the fluid was caused by the Lorentz force generated after location ③ and the re-acceleration was caused by the fluid entering a low-pressure zone after location ④. The downstream fluid velocities at location ⑥ were 0.284 m/s, 0.163 m/s, and 0.0946 m/s for the B = 0.5 T, 1.0 T, and 1.5 T cases, respectively. The lower velocity in the magnetized cases was characteristic of MR fluids and was caused by the increases in apparent viscosity and viscoelasticity of the fluid in the presence of magnetic fields [5].

Sinusoidal Pressure Profile
The model was shown to simulate the expected characteristics of MR fluid dampers in the constant pressure case. The response of the system to a variable sinusoidal pressure inlet condition with a peak-to-peak amplitude of 500 Pa was tested with magnetic field strengths of B = 0 T, 0.5 T, 1.0 T, and 1.5 T applied to the annular orifice. This was to simulate the response of the MR fluid damper when used in civil structures to dampen the vibration induced from strong wind conditions.
In Figure 10a-d, between locations ① to ③, there was a large pressure variation, and between locations ④ to ⑥, there was a much smaller pressure variation across all cases, indicating that the

Sinusoidal Pressure Profile
The model was shown to simulate the expected characteristics of MR fluid dampers in the constant pressure case. The response of the system to a variable sinusoidal pressure inlet condition with a peak-to-peak amplitude of 500 Pa was tested with magnetic field strengths of B = 0 T, 0.5 T, 1.0 T, and 1.5 T applied to the annular orifice. This was to simulate the response of the MR fluid damper when used in civil structures to dampen the vibration induced from strong wind conditions.
In Figure 10a-d, between locations 1 to 3 , there was a large pressure variation, and between locations 4 to 6 , there was a much smaller pressure variation across all cases, indicating that the fluid velocity was damped after the annular orifice. The largest pressure drops once the flow has developed across locations 3 and 4 at t = 1.32 s were 88.46 Pa, 131.99 Pa, 228.45 Pa, and 320.01 Pa for the B = 0 T, 0.5 T, 1.0 T, and 1.5 T cases, respectively. As expected, the magnetized cases had a much higher pressure drop across the annular orifice than the unmagnetized case. When the magnetic field strength was increased across the cases, the pressure drop across the annular orifice could be seen to almost linearly increase with the peak-to-peak pressure differences across locations 3 and 4 , increasing as the field strength was increased, as illustrated in Figure 11.
Meanwhile, in the constant pressure case, the pressure drop tapered off as the field strength was increased toward the saturation point of the particles; this could not be seen in the sinusoidal case. In the sinusoidal case, there was a constant increase as the magnetic field strength was increased from B = 1.0 T to 1.5 T. While this seems to suggest that the saturation field strength was not being approached, at B = 1.5 T, there should be asymptotic behavior. However, this did not occur in the sinusoidal case as the velocity of the fluid through the annular orifice was lower than in the constant pressure case. This was evident in the significantly lower pressure drops for the sinusoidal case shown in Figure 11, despite the theoretical mean inlet pressure of the two conditions being the same, and in the vertical shift of the inlet pressures profile. Although the inlet pressure profile was set to 500 + 250 sin(ωt), which should create a mean pressure identical to the constant pressure case, because of the resistance due to the magnetic field, the actual mean inlet pressure achieved was 434.54 Pa for the sinusoidal case. This lower mean pressure meant that the fluid velocity in the annular orifice did not approach the requisite velocity required for the current channeled through the fluid to result in saturation. This indicates that if the frequency of the wave was decreased, the pressure drop of the damper under sinusoidal loading would approach the values in the constant pressure case. Actual wind loading conditions would have a much lower frequency than the inlet condition used in the model and this accounted for the poorer performance of the damper when loaded with a sinusoidal profile. In general, it appears that the damper was more effective at damping lower frequency profiles as these allowed the flow in the annular orifice to fully develop. Regardless, the damper still achieved significant damping and the model was shown to capture the relevant characteristics of MR fluid dampers.
fluid velocity was damped after the annular orifice. The largest pressure drops once the flow has developed across locations ③ and ④ at t = 1.32 s were 88.46 Pa, 131.99 Pa, 228.45 Pa, and 320.01 Pa for the B = 0 T, 0.5 T, 1.0 T, and 1.5 T cases, respectively. As expected, the magnetized cases had a much higher pressure drop across the annular orifice than the unmagnetized case. When the magnetic field strength was increased across the cases, the pressure drop across the annular orifice could be seen to almost linearly increase with the peak-to-peak pressure differences across locations ③ and ④, increasing as the field strength was increased, as illustrated in Figure 11.
Meanwhile, in the constant pressure case, the pressure drop tapered off as the field strength was increased toward the saturation point of the particles; this could not be seen in the sinusoidal case. In the sinusoidal case, there was a constant increase as the magnetic field strength was increased from B = 1.0 T to 1.5 T. While this seems to suggest that the saturation field strength was not being approached, at B = 1.5 T, there should be asymptotic behavior. However, this did not occur in the sinusoidal case as the velocity of the fluid through the annular orifice was lower than in the constant pressure case. This was evident in the significantly lower pressure drops for the sinusoidal case shown in Figure 11, despite the theoretical mean inlet pressure of the two conditions being the same, and in the vertical shift of the inlet pressures profile. Although the inlet pressure profile was set to 500 + 250sin ( ), which should create a mean pressure identical to the constant pressure case, because of the resistance due to the magnetic field, the actual mean inlet pressure achieved was 434.54 Pa for the sinusoidal case. This lower mean pressure meant that the fluid velocity in the annular orifice did not approach the requisite velocity required for the current channeled through the fluid to result in saturation. This indicates that if the frequency of the wave was decreased, the pressure drop of the damper under sinusoidal loading would approach the values in the constant pressure case. Actual wind loading conditions would have a much lower frequency than the inlet condition used in the model and this accounted for the poorer performance of the damper when loaded with a sinusoidal profile. In general, it appears that the damper was more effective at damping lower frequency profiles as these allowed the flow in the annular orifice to fully develop. Regardless, the damper still achieved significant damping and the model was shown to capture the relevant characteristics of MR fluid dampers.

Seismic Response
Another possible application of MR fluid dampers is in the control of the vibrational characteristics of civil structures under seismic conditions. The suitability of the model in simulating the response of MR fluid dampers to real earthquake conditions was examined using a piecewise pressure inlet condition that approximated the 1994 Northridge earthquake, and the results are shown in Figure 12. As expected, the system experienced the largest pressure drop between locations ③ and ④ among all the equally-spaced locations, with a maximal difference of 419.84 Pa occurring at t = 4.89 s. It was evident that the system was effectively damped during all stages of the earthquake and exhibited the expected characteristics during both the large and small inlet pressure variations. In a similar manner to the sinusoidal pressure profile, the models also exhibited similar behavior with magnetic field strengths of B = 0.5 T and 1.5 T applied. Those results are omitted for conciseness.

Seismic Response
Another possible application of MR fluid dampers is in the control of the vibrational characteristics of civil structures under seismic conditions. The suitability of the model in simulating the response of MR fluid dampers to real earthquake conditions was examined using a piecewise pressure inlet condition that approximated the 1994 Northridge earthquake, and the results are shown in Figure 12. As expected, the system experienced the largest pressure drop between locations 3 and 4 among all the equally-spaced locations, with a maximal difference of 419.84 Pa occurring at t = 4.89 s. It was evident that the system was effectively damped during all stages of the earthquake and exhibited the expected characteristics during both the large and small inlet pressure variations. In a similar manner to the sinusoidal pressure profile, the models also exhibited similar behavior with magnetic field strengths of B = 0.5 T and 1.5 T applied. Those results are omitted for conciseness.

Conclusions
In this study, a numerical model was developed to simulate the behavior of magnetorheological (MR) fluid dampers by utilizing a computational fluid dynamics approach in conjunction with the magnetohydrodynamics model. Ohm's law and Maxwell's equations were coupled with the Navier-Stokes transport equations through the inclusion of the Lorentz force as a source term.
The simulations in this study were conducted based on a two-dimensional axisymmetric MR

Conclusions
In this study, a numerical model was developed to simulate the behavior of magnetorheological (MR) fluid dampers by utilizing a computational fluid dynamics approach in conjunction with the magnetohydrodynamics model. Ohm's law and Maxwell's equations were coupled with the Navier-Stokes transport equations through the inclusion of the Lorentz force as a source term.
The simulations in this study were conducted based on a two-dimensional axisymmetric MR fluid damper geometry with a magnetized area representing the annular orifice. Three test cases were conducted to investigate the capability of the present model, including a constant pressure case, a sinusoidal pressure case, and a seismic pressure case. Various plots of pressure and velocity were generated against time and displacement, and the model was found to accurately capture the macroscopic characteristics of MR fluids. In the constant pressure case, the MR fluid was found to have a lower velocity when the damper was in the on-state. In this state, it was also observed that the MR fluid decelerated after entering the annular orifice and then accelerated after exiting the annular orifice. This was due to the increased apparent viscosity and viscoelasticity of the fluid created by the resisting Lorentz force in the annular orifice.
It was found that MR fluids were most effective at damping constant pressure profiles with the highest pressure drop measured across the annular orifice. However, this was due to the high frequency of the sinusoidal profile and it is expected that the performance of the MR fluid between the two cases would be similar in real world conditions where the frequency of the load would be reduced.
Despite this, in both the sinusoidal and seismic pressure cases, there were large pressure variations measured before the magnetized annular orifice, followed by significantly smaller ones after the annular orifice, signifying that the MR fluid dampers were effective at damping these pressure profiles.
In conclusion, the presented numerical methodology and simulations in MR fluid modeling is beneficial in reducing the time and cost required in optimizing the design and performance of such devices. The present model is advantageous over the microscopic computational models, such as the discrete phase model, as it is significantly less computationally expensive to effectively model the behavior of MR fluid dampers in their full scale. Some limitations of the present study include the two-dimensional simplification of the MR fluid damper geometry. However, the underlying physical principles are identical between two-dimensional and three-dimensional simulations, thus the present study provides a solid foundation for larger-scale, three-dimensional simulations. The focus of the present model was the fluid behavior within a MR fluid damper when subjected to pressure and external magnetic fields. It may be combined with structure-focused computational models, such as the finite element method [47,48], to form a fluid structure interaction (FSI) study for further examinations of the effectiveness of MR fluid dampers in protecting various structures against vibration. Lastly, in order to further validate the computational model presented, physical experiments are to be conducted in the future for direct quantitative comparison. Funding: This research was supported by the Australian Research Council (ARC Project ID DP150101065 and DP160100021). The support of NVIDIA Corporation with the donation of the Titan Xp and Quadro P6000 GPUs used for this research is also gratefully acknowledged.