Two-Dimensional CFD Modeling of Hysteresis Behavior of MR Dampers

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. ,e 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. ,e 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. ,e model is solved by using the finite element method, and the movement of the piston is described by the moving mesh technique. ,e MR damper in a reference is simulated, and the model predictions show good agreement with the experimental data in the reference.

e 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][10][11][12] and more recent rapid developing physical modeling [13][14][15][16][17][18][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 threedimensional (3D) magnetostatic analysis with the flow analysis by using a bilinear Newtonian fluid with a fielddependent 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.

Problem Definition and Governing
Equations.e 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.
e working gap has a width of g 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.e damper has a stroke of s 0 , and a two-dimensional cylindrical coordinate system is adopted as shown in Figure 1.

Conservation of Momentum.
e flow of MR fluids in MR dampers is governed by the well-known Navier-Stokes equations: where i and j are tensor indexes, x 1 � r and x 3 � z in cylindrical coordinate systems, v i is the flow velocity component along x i -coordinate, b i is the body force along coordinate x i and will be neglected in this study, σ ij is the total stress tensor, (σ) ij,j is the divergence operator in the cylindrical coordinate system, and ρ is the density of MR fluids.

Conservation of Mass.
e 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.en, the conservation of mass is expressed simply by the vanishing divergence of the flow velocity: (2)

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 (z/zt).When the displacement gradients of the flow are large, the time derivatives, z/zt, 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 preand 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 by where _ ε is the total strain rate, _ ε k and _ ε η are the elastic and viscous parts of the total strain rate, k is the stiffness of the spring representing the compressibility of MR fluid, 2 Shock and Vibration 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 D ij is the strain rate tensor defined by where α is a small positive real number to avoid error of division by zero and it is set to 10 −8 ; the shear rate magnitude e 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.
en, the total stress tensor in equation ( 1) is where η N is the postyield viscosity and δ ij is the identity tensor.

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 zdirection, the interpolation is performed only in z-coordinates, with r-coordinates held constant.

Movement of the Piston. e displacement of the piston (z p ) is
where z 0 is the amplitude of sinusoidal piston displacement, f is the excitation frequency, and t is time.

Mesh of the Gap Domain.
e mesh of the working gap moves as a rigid translation in z-direction: where z gap is the Euler coordinates of the material points (fluid particles) in the gap domain.e mesh of the chambers is stretched or compressed, which can be described by the following linear interpolating functions: where z is the Euler coordinates of the material points occupying the chambers, Z top1 and Z bot1 are the Lagrangian coordinates at which the node displacements equal to the piston displacement, and Z top0 and Z bot0 are the Lagrangian coordinates at which the node displacement are zeros.For the coordinate system in Figure 1, these key coordinates are Z top1 � L/2, Z bot1 � −L/2, Z top0 � (s 0 + L/2), and Z bot0 � −(s 0 + L/2).

Boundary Conditions.
e 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 as where n is the outward normal vector of the boundary.

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

Evaluation of Damping Force.
e damping force, F, can be calculated by integrating the pressure along the piston boundaries as where A p is the effective area of the piston,

Implementation of the Proposed Model in COMSOL Multiphysics
e model equations are implemented in COMSOL Multiphysics software.
e constitutive equation, the mass conservation equation, and the momentum conservation equation are implemented, respectively, by using three builtin PDE modules which have the following general form PDE: where u is the vector of variables to be solved, Γ is the conservative flux vector, f is the source term vector, e a is the mass coefficient matrix, and d a 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 w), four shear stress components (tau11, tau13, tau22, and tau33), and a pressure field (p).ese 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.

PDE Module 1 in COMSOL: Constitutive Equation (Variables: tau11, tau13, tau22, and tau33).
In COMSOL, tensors are implemented as matrixes or vectors.e constitutive equation was implemented in the source term f as a fourth-order vector, i.e., where the domain-dependent auxiliary variable, eta, is the COMSOL implementation of equation ( 6). e damping coefficient matrix, d a , is the fourth-order identity matrix.
e conservative flux Γ is a zero vector and the mass coefficient, e a , is a zero square matrix.e whole input of the constitutive equation in COMSOL is shown in Figure 4.

PDE Module 2 in COMSOL: Mass Conservation Equation (Variable: p).
e mass conservation equation is implemented in the source term of the general form PDE as f � vdiv, where the auxiliary variable, vdiv, is defined in Table 1.All the coefficient matrices are zero matrices.e complete input of the mass conservation equation in COMSOL is shown in Figure 5.

PDE Module 3 in COMSOL: Momentum Conservation Equation (Variables: u and w).
e Cauchy stress is implemented in the conservative flux vector as Apparently, the damping coefficient matrix is simply a second-order unit matrix.
e complete input of the momentum conservation equation in COMSOL is shown in Figure 6.

Moving Mesh (ALE) Module: Motion Equation of Mesh.
e mesh motion can be readily implemented in COM-SOL by just prescribing the mesh motion equation in the built-in moving mesh (ALE) module, as shown in Figure 7.

Solver Configuration.
e built-in implicit backward differentiation formula (BDF) solver in COMSOL is used to solve the equation system for the model.e absolute solver tolerance is set to 1 × 10 −5 , and the maximum of iterations is raised up to 1 × 10 3 .All other solver configurations use default settings.Shock and Vibration 5 effect of mesh sensitivity of the model results.

Validation of the Proposed Model
ree 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). e   Shock and Vibration differently refined meshes are shown in Figure 8, with their close-up views presented in Figure 9.
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.

Experimental Validation.
e experimental data of the MR damper in the study by Yang [34] under different   Shock and Vibration operating conditions (exciting current: 0.25 A, 2 A; amplitude/frequency: 5.08 cm/0.05Hz, 2.54 cm/0.1 Hz, 1.27 cm/0.2Hz, 0.254 cm/1 Hz) will be used to validate the proposed model.
e 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.en, 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.Shock and Vibration 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). is 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.en, it accumulates constantly until it becomes large enough to overcome the yield strength and makes the fluid in the gap start to flow.e change in the damping force lags behind the change in the piston velocity, due to the compressibility of MR fluid.us, 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), excitation amplitude decreases with increasing excitation frequency.Accordingly, the compression or the pressure of the fluid in the pressurized chamber is small under high-frequency excitations, the MR damper works mainly in preyield zone.For this reason, it 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.
e 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.
e pressure in each chamber is nearly constant, with the large pressure gradient observed across    Shock and Vibration  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).

Conclusions
is 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.e solution of the model shows a slight mesh sensitivity.Except the time when the piston changes its direction, the highspeed 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.e hysteresis of the damping force with respect to the piston velocity is frequency-dependent, and it becomes apparently wider with increasing excitation frequency.
and K is the material parameter related to the compressibility of MR fluids.e viscosity of the modified Bingham model, η, is defined to be dependent on the yield stress of MR fluid, i.e., η � τ y D ij           1 − exp −m D ij              , i, j � 1, 3, MR fluids in the working gap, α, MR fluids in the chambers,

Figure 6 :
Figure 6: Input of momentum conservation equation in COMSOL.

Figure 5 :
Figure 5: Input of mass conservation equation in COMSOL.

Figure 7 :
Figure 7: Input of mesh motion equation in COMSOL.(a) Top chamber.(b) Working gap.(c) Bottom

Figure 10 :
Figure 10: Effect of mesh refinement on FEM results.

Table 1 :
Auxiliary variables defined in COMSOL for implementation convenience.According to the expression conventions of COMSOL, the spatial derivative of a variable, for example, zu/zr, is represented by u/r.e postyield viscosity implemented in COMSOL, etan, is 1.5 Pa•s. Note.