Investigating the characteristics of a magnetorheological fluid damper through CFD modeling

Computational Fluid Dynamics (CFD) analysis is conducted on mono-tube vehicle MRF damper investigated experimentally in a previous study. In this study, the fluid of the type MRF-132DG was inserted inside a damper of a car rear suspension system. The CFD analysis describes the fluid flow through the internal orifices between the compression and the rebound chambers. Averaged Navier–Stokes equations were solved by the SIMPLE method, and the RNG k-ε was used to model the turbulence at the fluid crossing through the orifices. All the CFD model boundary conditions’ values were set to the same values reported in the previous experimental study, except for the viscosity values. When varying the applied magnetic field density, the changes of MRF’s viscosity values were assessed by using a viscosity meter. Results showed a viscosity increase of 70% when the magnetic field excitation current was elevated from 0 A to 5 A. The damping forced and damping values were calculated using the rebound and compression static pressures obtained from the contour plots. It was also observed that the damping values exponentially increase with the increase in viscosity. The results of the CFD simulation were compared against those from the experiments, and good matching was observed.


Introduction
Noise and vibration are socio-environmental problems that call for solutions by engineers. Controlling noise and vibration requires a good knowledge of the dynamic systems and their main components, i.e., mass, stiffness, and damping. The primary application of a damper is to prevent discomfort from vibrations, outright structural failure, or even noise emissions (Li and Darby 2006). Vibration isolation countermeasures are categorized as passive, active, or semi-active. Passive control usually takes place within the design phase. It allows the damping device to absorb the vibration generated in the absence of a power supply and feedback sensors. Active control dampers require a feedback loop (i.e., sensor, actuators, a power supply, etc). In semi-active control, the damper's characteristics are modified by an active source of power that influences the system's dynamics (Sassi et al 2018). Viscous dampers are widely utilized as a countermeasure in different engineering applications such as passive vibration reduction, isolation, and control. The viscous damper consists of a stainless-steel piston and a steel cylinder separated into two chambers by the piston head. The damper's cylinder contains compressible hydraulic oil. As a result of absorbing the vibrating body's oscillating movement, the damper rod moves accordingly, forcing the oil inside the cylinder to move from one chamber to another (compression and rebound) through a small area (orifices). This movement leads to head loss because of the pressure difference as the fluid flow between the two chambers (Agrawal and Amjadian 2015). Viscous dampers are used in a wide range of applications. Figure 1 shows a schematic representation of a viscous damper along with it is internal parts.
Many studies investigate the use of viscous dampers in civil and automobile applications (Narkhede and Sinha 2014, Maikulal et al 2016, Galluzzi et al 2018, Lu et al 2018. Viscous dampers are considered as a high Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. energy dissipating density device (dissipating high energy compared to damper's small size) (Narkhede and Sinha 2014). In civil applications, viscous dampers enhance the structure performance for both shock load and long duration load (earthquake and wind load) (Narkhede and Sinha 2014). In automobiles, the vehicle suspension traditionally consists mainly of a spring and a damper mounted in parallel. A car suspension function maximizes the friction between the tires and the road surface, provides steering stability with proper handling, and ensures the passengers' comfort. Smart materials are composite materials whose properties dramatically change under external stimuli such as temperature, electric field, or magnetic field (Syam and Muthalif 2020). Magneto-Rheological (MRF) represents a particular type of smart materials whose rheological behavior is rapidly varied or controlled by applying a magnetic field (Jolly et al 1999). Without that excitation, the MRFs have low viscosity and behave as Newtonian fluids. However, under a magnetic field, the MRF micrometallic particles form a column-structure that resists shear deformation or flow. This change in the material induces a rapid increase in fluid viscosity or the development of a semisolid state. MRFs are usually modeled using the Bingham Plastic Model (Jolly et al 1999). An informative review about the MR fluid was provided by (Olabi and Grunwald 2007), discussing the MR fluid main components; base oil, metallic particles, and special additives. They also investigated the design and the use of MRFs in different applications through three modes; shear mode (used in clutches and brakes), valve mode (used in dampers), and squeeze mode (a combination mode) (Olabi and Grunwald 2007).
To quickly transform a standard viscous damper into an MRF one, the conventional hydraulic oil should be replaced by an MR fluid. Since there is no significant change in the damper's mechanical system, almost the same numerical and experimental investigation methods of the viscous damper are applicable for testing the modified MRF damper. In most experimental conditions, MRF-based dampers' dynamic testing is usually performed using a device that supplies a repetitive axial motion to move the piston rod back and forth. Axial fatigue machine, damper testing machine, shaker, or any mechanical slider attached to a motor will be a suitable actuator for such an experiment (Badri et al 2020). Guan et al (2019) developed a test bench for a twin-tube shock absorber consisting of a motor, a load rack, a static frame, and a 25 mm-radius crank. The crank transforms the rotation movement into reciprocating motion. The loading rack is connected to the crank through a linkage and then attached to the damper. A load cell is located at the bottom of the bench to record the damping force using a data acquisition (DAQ) system (Guan et al 2019). Hemantha and Arun (2018) developed a test setup for MRF dampers by connecting an exciting mechanical system and a linear variable differential transformer (LVDT) sensor to the damper's mounted end. A load cell was attached to the moving rod to record the variation of resistive force. The coils were excited by a DC power supply to provide the excitation current needed to create the magnetic field. All the sensors were connected to a DAQ device to record both displacement and force. For specific automobile applications, (Iglesias et al 2014) investigated, theoretically and experimentally, the use of MRF damper under external excitation. The experimental setup consisted of an oscillating testing machine with a maximum displacement of 22.5 cm. The frequency was varied using the 'sweep' feature to evaluate the damper response for a wide range of frequencies. The experimental results were used to validate the nonlinear hysteresis model proposed for predicting the damping force.
Numerous studies were conducted to explore the computational fluid dynamic (CFD) modeling approach in characterizing the viscous and MRF dampers. In particular, (Jiao et al 2018) developed an analytical model and an ADINA-based CFD model to investigate the nonlinear characteristics of micro-vibration fluid viscous damper, commonly used in aerospace applications. The fluid and the structure (orifice) models were created separately, while a harmonic motion was assigned as input to the structure model. A developed analytical model was used to validate the CFD by comparing the pressure gradient force for the two models. The CFD model was established to investigate the orifice openings numbers' effect in the damping coefficient value (Maikulal et al 2016). Three different scenarios were considered, starting with two orifice openings, then six, followed by having ten openings. The model results, for different velocity values, showed an essential decrease in the damping force for the ten-opening-orifice damper piston design. This particular investigation (Maikulal et al 2016) will be considered as a reference for the present study. However, several modifications will be needed to make the model suitable for MRF application. Goldasz, (2015) presented a theoretical model for a semi-active MRF damper. This model included different characteristics that contribute to the damping force, such as fluid chamber compressibility, fluid inertia, cylinder elasticity, and friction. Gedik et al (2012) studied the 2D steady laminar flow of an MR fluid between parallel plates using a CFD approach. The MR fluid flows through the plate's gap with velocities that change according to the magnetic flux variation. It was concluded that the fluid flow speed decreases with the increase of the external excitation magnetic values. Hemantha and Arun (2018) characterized a twin-tube MR damper experimentally and using coupled CFD-FE analysis. To have the ability to conduct a dynamic FE analysis along with a CFD analysis. CFD analysis was based on Navier-Stokes equations, and the FE used Maxwell equations. The field-dependent yield stress t b was obtain using the fitting correlation curve for each magnetic flux density value (B) (Han, Kim and Choi 2009). The magnetic-dependent yield stress t b value was then used in the Bingham plastic model to find the corresponding viscosity. The cylinder piston was assigned a sinusoidal motion using a moving mesh, to simulate the fluid cells movement along with time. The literature showed a wide range of CFD studies, but none of them rely on experimental results of viscosity values regarding magnetic field. Moreover, they all suffer from a common pitfall, which is the simplified geometries of the orifice shapes used in their CFD models. Sassi et al (2018) investigated the damping characteristics of MRF damper experimentally. They loaded the damper with dead weight and recorded the piston rod's travel time in compression and rebound cycles. Based on a linear force-velocity model, their damping values were obtained by using the travel distance of the rod, the travel time, and the driving force. In the present study a CFD model is used to simulate the behavior of the same MRF damper designed and used in (Sassi et al 2018). The same boundary conditions and parametric study provided in (Sassi et al 2018) are applied in the CFD model.
Hence, the present investigation introduces a CFD model for a twin-tube MRF damper experimentally investigated in previous work. This work's novelty lies in finding damping coefficients numerically depending on the experimental relationship between the viscosity and the applied magnetic field. This paper is, therefore, organized into six sections. In section 1, a brief introduction is presented. In section 2, the CFD fluid flow modeling methodology in finding the rebound and compression damping coefficient is developed. The experimental procedures used in finding the corresponded viscosity values to each magnetic field density value are illustrated in section 3. The computational analysis using ANSYS-17 Fluent flow numerical software is established in section 4. The orifice contour pressure gradient is shown and discussed in section 5. Section 6 summarizes the conclusions and the possible directions for future work.

Methodology
In this work, the approach considered to develop the numerical model is in line with the procedures previously used in (Maikulal et al 2016) and (Hemantha and Arun 2018). The CFD model reported in (Maikulal et al 2016) obtained the damping values for a viscous damper by introducing a fluid flow with a given velocity through an orifice then calculating its pressure drop. The fluid viscosity used in (Maikulal et al 2016) was considered constant since the damper's viscous fluid was deemed Newtonian. However, the MRF is a non-Newtonian fluid, and therefore its viscosity value may change with the variation of the applied magnetic field. A typical change of an MRF viscosity is depicted in the results generated by an FE-CFD coupled model (Hemantha and Arun 2018). The FE magnetostatic model was developed to calculate the magnetic flux density inside the MR damper flow gap. Different viscosities were found and then implemented in the CFD model based on the magnetic flux density values for each excitation current (Hemantha and Arun 2018). In this section, the MRF damper will be modeled using the model previously considered in (Maikulal et al 2016), which calculates the damping coefficient by calculating the forces acting on the compression and the rebound sides of the piston for constant velocities as shown in figure 2. Viscosity values illustrated in section 3 will then be used and identified in the material properities section in the Fluent flow model.
The relationships between force and velocity in compression and rebound (extension) strokes of the viscous damper are usually nonlinear. For moderate values of forces and velocities, the relationship between them could be reasonably considered as linear (Badri et al 2020): where F d is the damping force, C is the damping coefficient, and V is the piston-damper linear velocity. According to (Maikulal et al 2016), the damping force is defined as the difference between the rebound and the compression forces. When the damping force considers the fluid's friction effect, Guan et al (2019) proposed the following equation for a twin-tube shock absorber: where F d is the damping force, P c is the pressure in the compression chamber, P r is the pressure in the rebound chamber, A c is the cross sectional area of the orifice plate in the compression chamber, A r is the cross sectional area of the orifice plate in the rebound chamber, and F f is the frictional force. Similar to model [13], no friction effect will be counted in the present work. Although the friction will is negligible for low magnetic field values, the micro metallic chain structure will develop considerable friction at high magnetic flux density values. Both rebound and compression forces depend on the static pressure value acting on the orifice surface during the piston's oscillating motions. During the piston's reverse motion, the fluid is forced to flow from the rebound to the compression chamber, across a tiny orifice inside the piston. Therefore, a pressure build-up, known as rebound pressure, takes place around that orifice. The value of that rebound force F r is depicted in the following equation (3): where p r is the pressure in the rebound chamber, and A r is the cross-sectional area of the orifice plate in the rebound chamber.
On the other hand, during the damper rod's extension stroke, the fluid tends to move back from the compression to the rebound chamber. This motion generates a high-pressure value at the compression side. The compression force F c is given by: where p c is the pressure in the compression chamber, and A c is the cross-sectional area of the orifice plate in the compression chamber. The resultant damping force is due to the subtraction of the rebound force from the compressing force. In the case of constant velocity V, the damping coefficient for the rebound and the compression are given by: . . 7 2 The Fluent Flow ANSYS-17 provides different models to simulate turbulence. In this work, and for the sake of simplicity, the standard k-ε model was used. This choice was considered because this model does not include complex structures or flows involving (Jiao et al 2018). The standard k-ε is a two-equation model that includes two transported variables: k the turbulent kinetic energy, and ε the dissipating rate of the kinetic energy.
Although the realizable k-ε represents an improvement in the standard k-ε model (Harlow 2004), a negligible pressure difference was found when comparing the two models' results. For the material type setting, the fluid was considered viscous with a density of 3090 Kg m −3 . Six viscosity values were assigned to that fluid, as is going to be shown in section 3. The velocity was set to 0.08 m s −1 as mentioned in the experimental reference (Badri et al 2020). For coupling pressure and viscosity, a simple method was used along with a second-order upwind spatial solver to increase accuracy.

Experimental work
The relationship between the magnetic field excitation current and the fluid viscosity was assessed experimentally in two stages. First, we investigated the effect of the current on the magnetic field density. Then a viscosity meter, equipped with a magnetic excitation device, was used to obtain the corresponding viscosity value.

Magnetic field measurements
The magnetic field excitation system used in this study consists of a coil holder with twelve cores placed inside special openings precisely created for this purpose. The coil holder was 3D printed from Acrylonitrile-Butadiene-Styrene (ABS) material as shown in figure 3(a). Figure 3(b) shows the detailed dimensions of the magnetic excitation system. The excitation system surrounds the damper cylinder along a distance of 10 cm. The cores were arranged in four equally separated columns of three cores each (Sassi et al 2018). Copper wires were  wound around each coil, with 250 turns per set. Different pole arrangements were investigated, and the best performance was found to be obtained when two neighboring poles had opposite charge signs.
Once the coils were finished and correctly inserted into their holding support, the change in the magnetic flux density could be measured when the excitation current is varied. The magnetic flux density is recognized as: where m is the magnetic permeability of the core, N number of turns in the coil, and I is the excitation current. The magnetic flux density was measured using the previously described excitation system and a Gauss meter as in figure 4. The twelve cores were connected to a power supply in a particular arrangement that matches the requirements specified in (Sassi et al 2018). The meter needle was placed at a reference point located precisely in the center of the core holder. As the center is the most distant point from the cores (circumference), the reading taken at that point may have the minimum value of flux density because of the low air permeability.
The excitation current was periodically increased, from 0 A to 5 A, with a 1 A increment. The magnetic flux density was measured at each current value, as shown in figures 4 and 5.

Viscosity measurements
The dynamic viscosity value depends on the shear stress field-dependent t , b which also depends on the magnetic field density value (B). For the oil considered in the present study (MRF-132DG), its characteristic was  investigated previously by (Hidayatullah et al 2019). In the present study a rotational viscosity meter is used to obtain the viscosity values regarding the applied current. Rotational viscosity meters are devices that provide dynamic viscosity values by measuring the required torque to rotate a spindle with a known RPM inside a fluid container (Barnes 2001). Usually, the viscosity meter device consists of a motor, spindle, motor controller, rotational motion sensor, and force sensor. The spindle rotates inside the fluid container forcing the fluid to move accordingly, and an essential torque is generated. By measuring the force at the oil reservoir outer cylinder, the spindle RPM speed, and knowing the distance from the rotation axis, the dynamic viscosity is accordingly obtained. This study's advanced viscosity meter provides the viscosity value based on the torque and the RPM velocity of the spindle. The main feature for the actual device is that the torque and RPM are being measured using a servo system. The servo system determines the required current for moving the spindle inside the fluid with a set speed. When setting a velocity to the spindle, the fluid started to decrease its velocity. That's when the servo system determined the required extra current needed to contain the same rotation speed. For testing the viscosity of the MRF-132DG, a particular setup has been developed using the Viscosity meter (Viscometer 2300 RV). Such equipment comes with different spindles types, depending on the range of viscosity to be measured and its desired resolution. The spindle 'L1' was used to examine the MRF-132DG viscosity under different magnetic field values. The chosen spindle has a minimum recording value of (1 mPa.S), and a range of viscosity values (30-20000) mPa.s. Using a glass container to accommodate the fluid during viscosity testing may decrease the magnetic field effect and bring a reading error in viscosity value. Hence, the procedure was modified to measure the viscosity in situ by including the damper steel cylinder itself as a fluid container. The cores were placed around the damper, itself filled with MRF, as shown in figure 6. Once the spindle (L1) was placed inside the cylinder and the coils were excited, then the viscosity meter could be started to measure the viscosity values. The readings were done with an accuracy of 1%, and a repeatability error of 0.2%. The reading values are displayed in figure 7. The relationship between viscosity and the current variation is expressed by equation (9).

Orifice features
The previous studies' orifices were commonly considered small thickness cylinders with many holes (Guan et al 2019, Narkhede and Sinha 2014). The fundamental orifices' detailed small design features were frequently neglected. These small design features (like non-uniform edges) may substantially contribute to the pressure drop and the damping effect. In this study, the damper under investigation was previously used in a former experimental research work (Sassi et al 2018). The damper was disassembled, and its internal parts were checked to get their accurate dimensions and shapes. The graphical details are displayed in figures 8 and 9.

Computational analysis
The CFD numerical investigation requires three significant steps: pre-processing, solving, and post-processing. The pre-processing analysis starts by identifying all geometry dimensions to create the fluid domain. It is followed by an appropriate meshing of the elements, with a careful selection of element types and sizes. In the pre-processing analysis, initial phenomenon values of temperature and velocities may be calculated. The main steps in the solving step include assigning the appropriate model for the study while applying all the boundary conditions needed; material, viscosities, velocity, and pressure to the fluid domains. The post-processing phase will show the contour plot of the static pressure related to the compressing and rebound side of the viscous damper piston.

Meshing
In the CFD model, the meshing is concerned with the flow area, which started from the inlet through the throttling orifices, all the way to the outlet. Unstructured automatic fine meshing with 3D tetrahedral cells was used. The meshing smoothing was chosen to be high while having a normal curvature angle of 18°. The mesh responded to high gradient and curvature areas by providing more cells to those areas. Increasing the number of cells at the throttling area while having fewer elements at areas near the inlet and the outlet ensures accurate results with less computational time. The grid is then generated having 58169 elements with 12151 nodes, as shown in figure 10. The mesh quality is identified through different procedures. A mesh dependence study is the most common method in finding the best size of the element that should be used in the analysis. Although this method is giving a good result, it was found to be a time-consuming method. Three various parameters located in the meshing statistics of Fluent ANSYS provided an indication measure for the mesh quality. They are the skewness, the orthogonal quality, and the aspect ratio. Skewness is defined as the difference between the shape of the cell and an equilateral cell of equivalent volume. Its value ranges between 0 and 1. For a tetrahedral element, any value below 0.95 is acceptable (ANSYS Inc. 2008). The skewness of this study was found to be 0.72, which lead to a faster solution convergence. Another meshing quality parameter is the orthogonal quality, which refers to how close are the angles between adjacent faces of cells. The minimum value of the orthogonal quality was found to around 0.2. When such a value is above the minimum limit of 0.01, it is considered (ANSYS Inc. 2008). The aspect ratio for a tetrahedral element is the ratio between the sphere radius inscribed in the element and the sphere radius circumscribed in the same element. The maximum aspect ratio for the given mesh is equal to 14.77. Such a value is below the allowable max limit of 500, usually obtained when the mesh is of good quality (ANSYS Inc. 2008).

Simulation and modeling
After achieving the grid of the fluid area, the standard k-e model was used, as discussed in the Methodology section. The viscosity values, previously displayed in figure 6, were assigned to the liquid-water material characteristics to investigate the magnetic field density variation effect.
The rebound and compression pressures should be both considered for finding the damping coefficients. Both scenarios were considered in the CFD model to find the pressure at both sides of the orifices. Firstly, the fluid flow was supposed to be from the rebound chamber to the compression chamber with an initial velocity at the inlet boundary layer. This scenario will have a maximum pressure value at the orifice surface directly facing the inlet boundary. Reversing the flow direction by changing the inlet and the outlet boundary while maintaining the same velocity will result in the compression pressure. Figure 11 is showing the inlet, outlet, and three orifice holes of the fluid flow. Rebound and compression areas are name-selected for later pressure contour representation. A constant pressure condition was assigned to the outlet. In the meantime, as a boundary condition, the cylindrical partitions were modeled as walls. Only the inlet boundary is giving a constant value of 0.08 m s −1 representing the damper rod velocity that delivered the first damping value in the reference experimental study (Sassi et al 2018).

Results and discussion
After computing the results, the contours post-processing tool in CFD enabled showing the characteristic gradients, namely velocity, total pressure, dynamic pressure, and shear strain. In this study, only static pressure is required to identify the forces on both orifice sides (rebound, compression). Figure 12 displays an example of the static pressure contour plot for the rebound and compressing sides of the piston for an applied current of 1 A.
The rebound-compression forces are calculated for the maximum pressure value observed at the contour plot range. These forces are also related to the piston side surface areas. The sides areas are (´-3.71 10 4 m 2 ) at the compression side, and (´-3.40 10 4 m 2 ) at the rebound side, respectively. Consequently, the rebound force is more significant than the compression one because of the difference in the area values. As already explained in equations (5) and (6), the pressure drop generates different forces on both sides of the piston which difference,   divided by the inlet constant velocity, results in damping effect. All the results involved in the damping values are shown in table 1. Figure 13 depicts the variation of the rebound and compression damping coefficient obtained by CFD simulation, and those obtained experimentally (Sassi et al 2018), for different values of excitation current. At low currents, not exceeding 2.5A, one can see that the CFD model seems to provide overestimated numerical results. The difference between both sets of results is neighboring 10%. Such trend remains valid until the neighboring of 2.5 A, where the experimental results become higher. The gap between simulation and experimental results keeps increasing untill reaching 17%, for an excitation current of 5A.
In the present analysis, the CFD model takes into consideration the fluid chamber compressibility and the variation of fluid viscosity depending on the yield stress. However, the contribution of the internal friction in the damping force was omitted. The friction between the internal walls and the fluid, and between the internal particles of the fluid, are probably negligible, in the absence of external excitation. However, for higher values of magnetic fields, generated by a current more than 3 A, the friction effect could not be ignored anymore. In other words, at high values of current (3 A-5 A), the omission of internal friction, makes the CFD model fail to accurately capture the true behavior inside the damper.
Moreover, another reason behind the difference between simulation and experimental results is the way used in (Sassi et al 2018) to estimate the damping effect. Usually, the investigation of damping devices requires the use of appropriate testing machines equipped with speed controllers and force sensors to measure the damping coefficient. However, in the absence of such a machine, authors decided to load the damper with different dead weights and record the time in which the rod travels back and forth as in figure 14. Using the travel distance, the time was converted into speed, and consequently, the relationship between the force and the speed was easily obtained for both cycles of rebound and compression as hown in figures 14(a) and (b), respectively.
In their approach to measure the motion speed of the rod (Sassi et al 2018), authors considered the motion as uniform, which was not the case. Indeed, for a falling object, subjected to a resisitive force which is proportional to the speed, the equation of motion is of the type: Therefore, in the absence of external excitation, the damping constant is moderate, its effect on the rod's speed remains limited, and the gap between both sets of results is around 12% as in table 2. However, for a large value of current excitation neighboring 5 A, the damping coefficient is high enough to create a large slowdown of the rod's motion. For this case, the gap between true and estimated values of speed is around 25%.

Conclusion
This paper developed a CFD model which was used to characterize the change of the damping coefficient for an MRF damper, using ANSYS 17 Fluent. The change of the viscosity values was obtained experimentally, for five magnetic field values generated by the external excitation system, and then used to update the MR fluid features of the fluid flow model. The comparison of the CFD model and experimental damping coefficients results reveals the following: The MRF damper coefficient increases exponentially as the viscosity value of the rheological fluid increase.
For low values of excitation current, particularly not exceeding 2 A, the internal friction forces induced to the MRF are insignificant, and therefore the CFD model results are fairly close to the experimental ones.
For high excitation current and hence high magnetic field, the discrepancy between the model and the experimental results are significant. The way that experimental results were computed is one of the main reasons for the discripency between results.
In the absence of friction and cavitation force components, the CFD simulation generates accurate values that fairly match with the experimental results. In its actual state, the proposed model is only applicable for magnetic field values that give a viscosity range of 0.1 Pa.s −1.4 Pa.s.