#### Abstract

So far, most previous studies on the nonlinear hysteresis analysis of ER/MR dampers have been limited to one-dimensional modeling techniques. A two-dimensional (2D) axisymmetric CFD model of MR dampers is developed in this work. The main advantage of the proposed 2D model of MR dampers lies in that it can be used to predict dynamic behavior of MR devices of arbitrary geometries. The compressibility of MR fluids is the main factor responsible for the hysteresis behavior of MR dampers, and it has been rarely investigated in previous multidimensional modeling of MR damper. In our model, the compressibility of MR fluids is taken into account by the two-dimensional constitutive model which is extended from a previous one-dimensional physical model. The model is solved by using the finite element method, and the movement of the piston is described by the moving mesh technique. The MR damper in a reference is simulated, and the model predictions show good agreement with the experimental data in the reference.

#### 1. Introduction

Magnetorheological (MR) fluid damper shows great promises for semiactive vibration controls of civil structures [1–3], aerospace automobiles [4, 5], and aircraft [6–8] due to their advantages of mechanical simplicity, high dynamic range, low power requirements, large force capacity, and robustness.

The nonlinear hysteresis analysis of ER/MR dampers has been concentrated on one-dimensional modeling techniques because of their simplicity and easy physical interpretation, including quite maturing phenomenological modeling [9–12] and more recent rapid developing physical modeling [13–19]. Fruitful achievements have been made, which greatly deepens our understanding and broadens applications of ER/MR dampers. Notably, the hysteretic relationship between the piston velocity and damping force is found to be related to the compressibility and inertia of MR fluids [16, 17, 20, 21], while the hysteresis between the excitation current and damping force results from the hysteretic response of ferromagnetic materials [22].

Multidimensional analysis is receiving increasing attention because it is more accurate for general modeling of MR devices with arbitrary geometries. Using a generalized viscosity expressed in terms of a low degree polynomial function of the applied electric fields, Ursescu [23] simulated the ER flow in a channel with a prescribed inlet flow velocity and the free outlet. A two-dimensional (2D) computational fluid mechanics (CFD) model was constructed by Sahin et al. [24] for a MR valve having complex flow region. It showed apparent advantages of better agreements with the experimental data than the traditional 1D analytical model because the pressure drops in the elbows, expansions, entrance, and exit sections, which are neglected in the latter, are all considered in the CFD model. For MR valves with short orifice lengths, errors of analytical models become significantly large, especially at higher velocities, so a CFD model will be indispensable. Zhou and Bai [25] conducted a three-dimensional (3D) numerical FEM analysis for MRF seal technology applied on a circular cooler. By defining the apparent viscosity as the ratio of the fluid’s yield stress over the local shear rate, Gołdasz and Sapiński [26] studied a squeeze mode MR damper with a CFD model, and the well-known fact was confirmed that the compressive loads increase with the decreasing gap height. More recently, using the finite volume method on a two-dimensional moving grid, Syrakos et al. [27] successfully captured the hysteresis of a damper caused by the inertia of fluid under high-frequency loadings.

A multidimensional finite element model also helps us to better understand multiphysical behavior of MR dampers by including other physical aspects such as electromagnetic response. For example, Caterino et al. [28, 29] showed that the promptness of MR devices is almost exclusively dependent on the effectiveness of electromagnetic performances which directly determines the response time of MR dampers. Recently, Zheng et al. [22] established a more sophisticated multiphysics model which considered the magnetic, temperature, and flow fields, respectively, by the inverse Jiles–Atherton hysteresis model, conjugate heat transfer equations, and Navier–Stokes equations. With the piston movements described by a deformed mesh, Case et al. [30] developed a multiphysics finite element model for a small scale MR damper, and the coupling of the magnetic field with the flow field was achieved through a modified Bingham plastic material model which relates the fluid’s dynamic viscosity to the intensity of the induced magnetic field. A similar work was conducted by Sternberg et al. [31] which couples a three-dimensional (3D) magnetostatic analysis with the flow analysis by using a bilinear Newtonian fluid with a field-dependent dynamic viscosity.

As mentioned above, the inertia and compressibility of fluids are the two main factors affecting the hysteresis behavior of MR dampers. However, the flow analyses in all the above 2D/3D models are still limited to the inertia effect of fluid, with the main differences lying in the different modifications or regulations of the 2D/3D Bingham model. So far, the 2D modeling of the hysteresis of ER/MR dampers arisen from the compressibility of fluids has not been reported yet, and this motivates the present work.

#### 2. 2D Moving Mesh Finite Element Modeling of Hysteresis of MR Dampers

##### 2.1. Problem Definition and Governing Equations

The layout of a MR damper is shown in Figure 1. A piston of radius *R*_{p}, with a shaft of radius *R*_{s}, at its center, divides the damper into the top and bottom chambers. Movement of the piston causes MR fluids to flow through the annular gap between the house cylinder and the piston. The working gap has a width of and effective length of *L*. Because this work focuses on the flow analysis, the whole space filled with MR fluid between the piston and the house cylinder is the physical domain of the finite element model to be constructed later, with the magnetic field related geometries such as coils omitted. The damper has a stroke of s_{0}, and a two-dimensional cylindrical coordinate system is adopted as shown in Figure 1.

##### 2.2. Conservation of Momentum

The flow of MR fluids in MR dampers is governed by the well-known Navier–Stokes equations:where and are tensor indexes, and in cylindrical coordinate systems, is the flow velocity component along -coordinate, is the body force along coordinate and will be neglected in this study, is the total stress tensor, is the divergence operator in the cylindrical coordinate system, and is the density of MR fluids.

##### 2.3. Conservation of Mass

The effect of the compressibility of MR fluid, which will be taken into account in the constitutive model, is assumed to be negligible on the fluid density. Then, the conservation of mass is expressed simply by the vanishing divergence of the flow velocity:

##### 2.4. A 2D Constitutive Model Inspired from 1D Physical Model

Problems of fluid viscoelasticity is notorious for being extremely difficult and challenging because of numerical difficulties such as the advective nature of the constitutive equations and the interaction of multiple discrete unknown fields (viscoelastic stress, velocity, and pressure) [32].

In the case where the displacement gradients of the flow are infinitesimal, the time derivatives on the stress and strain tensors are simply the common partial derivative with respect to the time (). When the displacement gradients of the flow are large, the time derivatives, , must be replaced with the convected derivatives of the stress and strain tensors, also known as contravariant convected time derivative [33]. Because the preyield MR fluid behaves like a rigid solid, the displacement gradients are infinitesimal. And for simplicity, we assume that the stress variation with time for the whole MR fluid flow, including both the pre- and postyield zones, can still be approximated by a common partial time derivative.

1D phenomenal or physical models [15, 16, 20] have been successfully used to predict the hysteresis of MR dampers, so it is naturally to expect a 2D extension from them. Guo et al. [15] developed a simple physical model as shown in Figure 2, which can be mathematically described bywhere is the total strain rate, and are the elastic and viscous parts of the total strain rate, is the stiffness of the spring representing the compressibility of MR fluid, and is the viscosity characterizing the quasistatic model of MR dampers.

Equation (3) can be rewritten as

Inspired from the above 1D constitutive for the physical modeling of hysteresis of MR dampers, the following 2D extension constitutive model in terms of the symmetric viscous stress tensor is proposed:where is the strain rate tensor defined by , is the symmetric viscous tensor () and is the material parameter related to the compressibility of MR fluids. The viscosity of the modified Bingham model, , is defined to be dependent on the yield stress of MR fluid, i.e.,where is a small positive real number to avoid error of division by zero and it is set to ; the shear rate magnitude is defined by

The simplicity of the constitutive model (equation (5)) comes at the cost of use of a nonobjective derivative, and the validation of this constitutive model will be demonstrated later by its comparison with experimental data.

Then, the total stress tensor in equation (1) iswhere is the postyield viscosity and is the identity tensor.

##### 2.5. Moving Mesh

Moving mesh technique is required to include the arbitrary motion of the piston. Because the ratio of width to length of the working gap size is extremely large, mesh movement should be carefully handled. Linear interpolating moving equations are adopted in this work to avoid severe element distortions, thus improving the converging ability of the finite element model. Because the movement of the piston involves only the movement in *z*-direction, the interpolation is performed only in *z*-coordinates, with *r*-coordinates held constant.

###### 2.5.1. Movement of the Piston

The displacement of the piston () iswhere is the amplitude of sinusoidal piston displacement, is the excitation frequency, and is time.

###### 2.5.2. Mesh of the Gap Domain

The mesh of the working gap moves as a rigid translation in *z*-direction:where is the Euler coordinates of the material points (fluid particles) in the gap domain.

###### 2.5.3. Mesh of Chamber

The mesh of the chambers is stretched or compressed, which can be described by the following linear interpolating functions:where is the Euler coordinates of the material points occupying the chambers, and are the Lagrangian coordinates at which the node displacements equal to the piston displacement, and and are the Lagrangian coordinates at which the node displacement are zeros. For the coordinate system in Figure 1, these key coordinates are .

##### 2.6. Boundary Conditions

The Dirichlet boundary conditions (BCs) for the flow velocity are shown in Figure 3. And natural (i.e., zero flux) BCs are implemented for the viscous tensor and pressure aswhere is the outward normal vector of the boundary.

##### 2.7. Initial Conditions

The initial values of viscous stress, the pressure, and the flow velocity are all assumed to be zeros.

##### 2.8. Evaluation of Damping Force

The damping force, *F*, can be calculated by integrating the pressure along the piston boundaries aswhere is the effective area of the piston, .

#### 3. Implementation of the Proposed Model in COMSOL Multiphysics

The model equations are implemented in COMSOL Multiphysics software. The constitutive equation, the mass conservation equation, and the momentum conservation equation are implemented, respectively, by using three built-in PDE modules which have the following general form PDE:where is the vector of variables to be solved, **Γ** is the conservative flux vector, is the source term vector, is the mass coefficient matrix, and is the damping coefficient matrix. For the model of MR dampers in this work, a total of 7 variables need to be defined in COMSOL, i.e., two velocity components (*u* and ), four shear stress components (tau11, tau13, tau22, and tau33), and a pressure field (p). These three PDE modules share the same fluid domain. For the convenience of implementation, we defined several auxiliary variables as listed in Table 1.

Detailed implementation of the proposed model is given in the following sections.

##### 3.1. PDE Module 1 in COMSOL: Constitutive Equation (Variables: tau11, tau13, tau22, and tau33)

In COMSOL, tensors are implemented as matrixes or vectors. The constitutive equation was implemented in the source term as a fourth-order vector, i.e.,where the domain-dependent auxiliary variable, eta, is the COMSOL implementation of equation (6). The damping coefficient matrix, , is the fourth-order identity matrix. The conservative flux **Γ** is a zero vector and the mass coefficient, , is a zero square matrix. The whole input of the constitutive equation in COMSOL is shown in Figure 4.

**(a)**

**(b)**

##### 3.2. PDE Module 2 in COMSOL: Mass Conservation Equation (Variable: )

The mass conservation equation is implemented in the source term of the general form PDE as , where the auxiliary variable, , is defined in Table 1. All the coefficient matrices are zero matrices. The complete input of the mass conservation equation in COMSOL is shown in Figure 5.

##### 3.3. PDE Module 3 in COMSOL: Momentum Conservation Equation (Variables: and )

The Cauchy stress is implemented in the conservative flux vector as

The convective part of the time derivative of velocity is included in the source term as

Apparently, the damping coefficient matrix is simply a second-order unit matrix. The complete input of the momentum conservation equation in COMSOL is shown in Figure 6.

##### 3.4. Moving Mesh (ALE) Module: Motion Equation of Mesh

The mesh motion can be readily implemented in COMSOL by just prescribing the mesh motion equation in the built-in moving mesh (ALE) module, as shown in Figure 7.

**(a)**

**(b)**

**(c)**

##### 3.5. Solver Configuration

The built-in implicit backward differentiation formula (BDF) solver in COMSOL is used to solve the equation system for the model. The absolute solver tolerance is set to , and the maximum of iterations is raised up to . All other solver configurations use default settings.

#### 4. Validation of the Proposed Model

##### 4.1. Mesh Sensitivity Analysis

Before validating the FEM model by using experimental data, it is necessary to study the effect of mesh sensitivity of the model results. Three refinement levels, i.e., coarse, normal, and fine meshes with 13656, 50698, and 127018 degrees of freedom (DOFs), respectively (Figure 8), were used to simulate the performance of the MR damper under a medium-frequency sinusoidal displacement excitation (2.54 cm, 0.1 Hz, 2 A). The differently refined meshes are shown in Figure 8, with their close-up views presented in Figure 9.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

As shown in Figure 10, the difference between the model results is negligible, and we will choose the normal mesh in the following validation of the model. It took about 8 hours to solve the model with the fine mesh on a personal computer (i5 CPU, 4G RAM), while the models with coarse or normal meshes only took less than 1 hour.

**(a)**

**(b)**

##### 4.2. Experimental Validation

The experimental data of the MR damper in the study by Yang [34] under different operating conditions (exciting current: 0.25 A, 2 A; amplitude/frequency: 5.08 cm/0.05 Hz, 2.54 cm/0.1 Hz, 1.27 cm/0.2 Hz, 0.254 cm/1 Hz) will be used to validate the proposed model. The main structural parameters of the damper are listed in Table 2. By fitting the experimental data in the case of 5.08 cm, 0.05 Hz, 0.25 A, the material constants, *K* and *m*, are identified to be 0.2 × 10^{4} and 0.25, respectively. Then, the performances of the damper in the other seven cases are predicted with this set of values. As shown in Figures 11–14, the model predictions are in good agreement with the experimental data, so the proposed model is reliable.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

Due to an applied magnetic field, the MR fiuid in the working gap consists of chain structures, leading to a very large preyield viscosity (i.e., the shear yield strength). This rheology character of MR fluids gives rise to the stick-slip friction-like response of the damping force with respect to the piston velocity. When a MR damper works in the preyield zone, for example, when the piston slowly starts moving from the rest, the pressure in the pressurized chamber gradually grows due to the compressibility of MR fluid in chamber. Then, it accumulates constantly until it becomes large enough to overcome the yield strength and makes the fluid in the gap start to flow. The change in the damping force lags behind the change in the piston velocity, due to the compressibility of MR fluid. Thus, a hysteresis between the damping force and the piston velocity will be observed when a MR damper is working in the preyield zone, even if the pressure in the pressurized chamber is relatively low.

While keeping the maximum piston velocity unchanged (about 0.015 m/s as shown in Figures 11–14), the excitation amplitude decreases with increasing excitation frequency. Accordingly, the compression or the pressure accumulation of the fluid in the pressurized chamber is small under high-frequency excitations, and the MR damper works mainly at in preyield zone. For this reason, it is indeed observed in Figures 11–14 that the hysteresis of the damping force with respect to the piston velocity becomes apparently wider with increasing excitation frequency.

The typical field results of the flow velocity and pressure are presented in Figures 15–17. As shown in Figure 15, except the time when the piston changes its direction, the high-speed flow mainly occurs in the working gap due to the small ration of the gap width to the piston radius. The pressure in each chamber is nearly constant, with the large pressure gradient observed across the working gap (Figure 16). Moreover, the flow pattern can be more clearly seen from the steam lines in Figure 17, and large vortexes are observed when the piston reverses its direction (*t* = *T*/4).

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

#### 5. Conclusions

This paper develops a two-dimensional axisymmetric finite element model for simulating the dynamic behavior of MR dampers. By describing the MR fluid as a viscoelastic material, the model successfully captures the hysteresis in the damping force with respect to the piston velocity. The solution of the model shows a slight mesh sensitivity. Except the time when the piston changes its direction, the high-speed flow occurs in the working gap due to the small ratio of the gap width to the piston radius. Large vortexes are generated in the chambers of MR dampers during the piston reversing its direction. The hysteresis of the damping force with respect to the piston velocity is frequency-dependent, and it becomes apparently wider with increasing excitation frequency.

#### Data Availability

All data generated or analyzed during this study are included in this paper.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (no. 51308450), Foundation of Key Laboratory of Structures Dynamic Behavior and Control (Ministry of Education) in Harbin Institute of Technology (HITCE201708), and Research Foundation of Xi’an University of Architecture and Technology (no. RC1368).