# Two-Dimensional CFD Modeling of Hysteresis Behavior of MR Dampers

Authors: Pengfei Guo; Jing Xie
Journal: Shock and Vibration
 (2019)

---

## 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.


---

## Body

## 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 Figure1. A piston of radius Rp, with a shaft of radius Rs, 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 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. The damper has a stroke of s0, and a two-dimensional cylindrical coordinate system is adopted as shown in Figure 1.Figure 1
Typical configuration of a MR damper.


### 2.2. Conservation of Momentum

The flow of MR fluids in MR dampers is governed by the well-known Navier–Stokes equations:(1)ρ∂vi∂t+vj∂vi∂xj=σij,j+ρbi,i=1,3,where i and j are tensor indexes, x1=r and x3=z in cylindrical coordinate systems, vi is the flow velocity component along xi-coordinate, bi is the body force along coordinate xi 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.


### 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)vi,i=0,i=1,3.


### 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 (∂/∂t). When the displacement gradients of the flow are large, the time derivatives, ∂/∂t, 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 by(3)ε˙=ε˙k+ε˙η=s˙k+sη,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, and η is the viscosity characterizing the quasistatic model of MR dampers.Figure 2
A physical model of MR dampers.Equation (3) can be rewritten as(4)s˙=kε˙−kηs.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:(5)∂τij∂t=KDij−Kητij,i,j=1,3,where Dij is the strain rate tensor defined by Dij=vi,j+vj,i, τij is the symmetric viscous tensor (τij=τji) and K 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.,(6)η=τyDij1−exp−mDij,i,j=1,3,MR fluids in the working gap,α,MR fluids in the chambers,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 Dij is defined by(7)Dij=12DijDij,i,j=1,3.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) is(8)σij=2ηNDij−pδij+τij,i,j=1,3,where ηN is the postyield viscosity and δij 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 inz-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 (zp) is(9)zp=z0sin2πft,where z0 is the amplitude of sinusoidal piston displacement, f is the excitation frequency, and t is time.


#### 2.5.2. Mesh of the Gap Domain

The mesh of the working gap moves as a rigid translation inz-direction:(10)zgap=zp,zgap∈gap domain,where zgap 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:(11)z=ztop,z∈top chamber domain,zbot,z∈bottom chamber domain,ztop=zpZ−Ztop0Ztop1−Ztop0,zbot=zpZ−Zbot0Zbot1−Zbot0,where z is the Euler coordinates of the material points occupying the chambers, Ztop1 and Zbot1 are the Lagrangian coordinates at which the node displacements equal to the piston displacement, and Ztop0 and Zbot0 are the Lagrangian coordinates at which the node displacement are zeros. For the coordinate system in Figure 1, these key coordinates are Ztop1=L/2,Zbot1=−L/2,Ztop0=s0+L/2,andZbot0=−s0+L/2.


### 2.6. Boundary Conditions

The Dirichlet boundary conditions (BCs) for the flow velocity are shown in Figure3. And natural (i.e., zero flux) BCs are implemented for the viscous tensor and pressure as(12)∂τij∂n=0,∂p∂n=0,where n is the outward normal vector of the boundary.Figure 3
Boundary conditions for flow velocity.


### 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 as(13)F=∫z=L/2p2πrdr−∫z=−L/2p2πrdrAp,where Ap is the effective area of the piston, Ap=πRp2−Rs2.


## 2.1. Problem Definition and Governing Equations

The layout of a MR damper is shown in Figure1. A piston of radius Rp, with a shaft of radius Rs, 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 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. The damper has a stroke of s0, and a two-dimensional cylindrical coordinate system is adopted as shown in Figure 1.Figure 1
Typical configuration of a MR damper.


## 2.2. Conservation of Momentum

The flow of MR fluids in MR dampers is governed by the well-known Navier–Stokes equations:(1)ρ∂vi∂t+vj∂vi∂xj=σij,j+ρbi,i=1,3,where i and j are tensor indexes, x1=r and x3=z in cylindrical coordinate systems, vi is the flow velocity component along xi-coordinate, bi is the body force along coordinate xi 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.


## 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)vi,i=0,i=1,3.


## 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 (∂/∂t). When the displacement gradients of the flow are large, the time derivatives, ∂/∂t, 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 by(3)ε˙=ε˙k+ε˙η=s˙k+sη,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, and η is the viscosity characterizing the quasistatic model of MR dampers.Figure 2
A physical model of MR dampers.Equation (3) can be rewritten as(4)s˙=kε˙−kηs.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:(5)∂τij∂t=KDij−Kητij,i,j=1,3,where Dij is the strain rate tensor defined by Dij=vi,j+vj,i, τij is the symmetric viscous tensor (τij=τji) and K 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.,(6)η=τyDij1−exp−mDij,i,j=1,3,MR fluids in the working gap,α,MR fluids in the chambers,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 Dij is defined by(7)Dij=12DijDij,i,j=1,3.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) is(8)σij=2ηNDij−pδij+τij,i,j=1,3,where ηN is the postyield viscosity and δij 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 inz-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 (zp) is(9)zp=z0sin2πft,where z0 is the amplitude of sinusoidal piston displacement, f is the excitation frequency, and t is time.


### 2.5.2. Mesh of the Gap Domain

The mesh of the working gap moves as a rigid translation inz-direction:(10)zgap=zp,zgap∈gap domain,where zgap 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:(11)z=ztop,z∈top chamber domain,zbot,z∈bottom chamber domain,ztop=zpZ−Ztop0Ztop1−Ztop0,zbot=zpZ−Zbot0Zbot1−Zbot0,where z is the Euler coordinates of the material points occupying the chambers, Ztop1 and Zbot1 are the Lagrangian coordinates at which the node displacements equal to the piston displacement, and Ztop0 and Zbot0 are the Lagrangian coordinates at which the node displacement are zeros. For the coordinate system in Figure 1, these key coordinates are Ztop1=L/2,Zbot1=−L/2,Ztop0=s0+L/2,andZbot0=−s0+L/2.


## 2.5.1. Movement of the Piston

The displacement of the piston (zp) is(9)zp=z0sin2πft,where z0 is the amplitude of sinusoidal piston displacement, f is the excitation frequency, and t is time.


## 2.5.2. Mesh of the Gap Domain

The mesh of the working gap moves as a rigid translation inz-direction:(10)zgap=zp,zgap∈gap domain,where zgap 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:(11)z=ztop,z∈top chamber domain,zbot,z∈bottom chamber domain,ztop=zpZ−Ztop0Ztop1−Ztop0,zbot=zpZ−Zbot0Zbot1−Zbot0,where z is the Euler coordinates of the material points occupying the chambers, Ztop1 and Zbot1 are the Lagrangian coordinates at which the node displacements equal to the piston displacement, and Ztop0 and Zbot0 are the Lagrangian coordinates at which the node displacement are zeros. For the coordinate system in Figure 1, these key coordinates are Ztop1=L/2,Zbot1=−L/2,Ztop0=s0+L/2,andZbot0=−s0+L/2.


## 2.6. Boundary Conditions

The Dirichlet boundary conditions (BCs) for the flow velocity are shown in Figure3. And natural (i.e., zero flux) BCs are implemented for the viscous tensor and pressure as(12)∂τij∂n=0,∂p∂n=0,where n is the outward normal vector of the boundary.Figure 3
Boundary conditions for flow velocity.


## 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 as(13)F=∫z=L/2p2πrdr−∫z=−L/2p2πrdrAp,where Ap is the effective area of the piston, Ap=πRp2−Rs2.


## 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:(14)ea∂2u∂t2+da∂u∂t+∇⋅Γ=f,where u is the vector of variables to be solved, Γ is the conservative flux vector, f is the source term vector, ea is the mass coefficient matrix, and da 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). These three PDE modules share the same fluid domain. For the convenience of implementation, we defined several auxiliary variables as listed in Table 1.Table 1
Auxiliary variables defined in COMSOL for implementation convenience.
Physical meaning
Name
Expression in COMSOL
Shear rate
D11
2 ∗ur
D13
uz + wr
D22
2 ∗u/r
D33
2 ∗ wz
Magnitude of shear rate
sr
Sqrt (0.5 ∗ 4 ∗ (u/r)2 + 2∗ (uz + wr)2 + 4 ∗ (wz)2 + eps)
Total stress
Sigma11
etan ∗ D11 + tau11-p
Sigma13
etan ∗ D13 + tau13
Sigma22
etan ∗ D22 + tau22-p
Sigma33
etan ∗ D33 + tau33-p
Divergence of velocity
vdiv
u/r + ur + wz
Note. According to the expression conventions of COMSOL, the spatial derivative of a variable, for example, ∂u/∂r, is represented by u/r. The postyield viscosity implemented in COMSOL, etan, is 1.5 Pa·s.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 termf as a fourth-order vector, i.e.,(15)f=K∗D11−K/eta∗tau11K∗D13−K/eta∗tau13K∗D22−K/eta∗tau22K∗D33−K/eta∗tau33,where the domain-dependent auxiliary variable, eta, is the COMSOL implementation of equation (6). The damping coefficient matrix, da, is the fourth-order identity matrix. The conservative flux Γ is a zero vector and the mass coefficient, ea, is a zero square matrix. The whole input of the constitutive equation in COMSOL is shown in Figure 4.Figure 4
Input of constitutive equation in COMSOL. (a) Conservative flux term. (b) Source term, damping coefficient matrix, and mass coefficient matrix.
(a)
(b)


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

The mass conservation equation is implemented in the source term of the general form PDE asf=vdiv, where the auxiliary variable, vdiv, 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.Figure 5
Input of mass conservation equation in COMSOL.


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

The Cauchy stress is implemented in the conservative flux vector as(16)Γ=−sigma11−sigma13−sigma13−sigma33.The convective part of the time derivative of velocity is included in the source term as(17)f=−rho∗u∗ur+w∗uz+sigma11−sigma22/r−rho∗u∗wr+w∗wz+sigma13/r.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 Figure6.Figure 6
Input of momentum conservation equation in COMSOL.


### 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 Figure7.Figure 7
Input of mesh motion equation in COMSOL. (a) Top chamber. (b) Working gap. (c) Bottom chamber.
(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 to1×10−5, and the maximum of iterations is raised up to 1×103. All other solver configurations use default settings.


## 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 termf as a fourth-order vector, i.e.,(15)f=K∗D11−K/eta∗tau11K∗D13−K/eta∗tau13K∗D22−K/eta∗tau22K∗D33−K/eta∗tau33,where the domain-dependent auxiliary variable, eta, is the COMSOL implementation of equation (6). The damping coefficient matrix, da, is the fourth-order identity matrix. The conservative flux Γ is a zero vector and the mass coefficient, ea, is a zero square matrix. The whole input of the constitutive equation in COMSOL is shown in Figure 4.Figure 4
Input of constitutive equation in COMSOL. (a) Conservative flux term. (b) Source term, damping coefficient matrix, and mass coefficient matrix.
(a)
(b)


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

The mass conservation equation is implemented in the source term of the general form PDE asf=vdiv, where the auxiliary variable, vdiv, 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.Figure 5
Input of mass conservation equation in COMSOL.


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

The Cauchy stress is implemented in the conservative flux vector as(16)Γ=−sigma11−sigma13−sigma13−sigma33.The convective part of the time derivative of velocity is included in the source term as(17)f=−rho∗u∗ur+w∗uz+sigma11−sigma22/r−rho∗u∗wr+w∗wz+sigma13/r.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 Figure6.Figure 6
Input of momentum conservation equation in COMSOL.


## 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 Figure7.Figure 7
Input of mesh motion equation in COMSOL. (a) Top chamber. (b) Working gap. (c) Bottom chamber.
(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 to1×10−5, and the maximum of iterations is raised up to 1×103. 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 (Figure8), 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.Figure 8
FEM meshes with different refinement levels. (a) Coarse mesh (13656 DOFs). (b) Normal mesh (50698 DOFs). (c) Fine mesh (127018 DOFs).
(a)
(b)
(c)Figure 9
Close view of FEM meshes. (a) Coarse mesh. (b) Normal mesh. (c) Fine mesh.
(a)
(b)
(c)As shown in Figure10, 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.Figure 10
Effect of mesh refinement on FEM results.
(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 × 104 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.Table 2
The configuration of the damper in Yang’ study [34].
Variable
Physical meaning
Value
G
Gap width
2.057 mm
S
0
Stroke
80 mm
R
p
Piston radius
101.63 mm
R
s
Piston shaft radius
40 mm
L
Effective axial pole length
84.428 mm
η
N
Postyield viscosity of MR fluid
1.5 Pa·sFigure 11
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 5.08 cm, 0.05 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 12
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 2.54 cm, 0.1 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 13
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 14
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 0.254 cm, 1 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(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 Figures11–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 Figures15–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).Figure 15
Magnitude of flow velocity (m/s) in the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(a)
(b)
(c)Figure 16
Pressure field (Pa) in the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(a)
(b)
(c)Figure 17
Stream lines within the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(a)
(b)
(c)


## 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 (Figure8), 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.Figure 8
FEM meshes with different refinement levels. (a) Coarse mesh (13656 DOFs). (b) Normal mesh (50698 DOFs). (c) Fine mesh (127018 DOFs).
(a)
(b)
(c)Figure 9
Close view of FEM meshes. (a) Coarse mesh. (b) Normal mesh. (c) Fine mesh.
(a)
(b)
(c)As shown in Figure10, 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.Figure 10
Effect of mesh refinement on FEM results.
(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 × 104 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.Table 2
The configuration of the damper in Yang’ study [34].
Variable
Physical meaning
Value
G
Gap width
2.057 mm
S
0
Stroke
80 mm
R
p
Piston radius
101.63 mm
R
s
Piston shaft radius
40 mm
L
Effective axial pole length
84.428 mm
η
N
Postyield viscosity of MR fluid
1.5 Pa·sFigure 11
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 5.08 cm, 0.05 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 12
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 2.54 cm, 0.1 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 13
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(a)
(b)Figure 14
Model predictions (lines) and experimental data (Yang [34]) of hysteretic behavior of MR dampers under sinusoidal displacement excitations of 0.254 cm, 1 Hz. (a) Displacement-force hysteresis. (b) Velocity-force hysteresis.
(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 Figures11–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 Figures15–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).Figure 15
Magnitude of flow velocity (m/s) in the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(a)
(b)
(c)Figure 16
Pressure field (Pa) in the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(a)
(b)
(c)Figure 17
Stream lines within the MR damper under sinusoidal displacement excitations of 1.27 cm, 0.2 Hz, 2 A: (a)t = T/12; (b) t = 2T/12; (c) t = 3T/12.
(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.



---

*Source: 9383047-2019-04-10.xml*