Predicting the dynamic behaviour of hydrobushings

A novel and promising approach for the prediction of the dynamic stiffness of hydrobushings is presented, combining Finite Element and CFD methods. The rubber structure of the mount is modelled in ABAQUS and the flow of fluid through the inertia track is calculated in FLUENT. The obtained results from the latter simulation are incorporated in the finite element code for the final stiffness prediction. The calculation is very sensitive to both rubber and fluid properties. The dynamic behaviour of rubber material has accurately been characterised with a new simple shear specimen in a forced non-resonant test. Satisfactory results are obtained when comparing numerical simulations to experimental tests in a practical application. Discrepancies between simulations and tests are mainly due to the simplifications assumed when creating the model. Nevertheless, stiffness of the mount is well predicted and so is the damping, although the frequency at which its maximum value is achieved is underestimated by 4–6 Hz, result that could be improved if non-stationary boundary conditions were considered when solving the fluid flow and incorporating it to the finite element code.


Introduction
Rubber is, without doubt, the most widely used and cheap solution when the connection between two parts of a structure, a vehicle or a machine, requires the combination of stiffness and vibration isolation.In automobiles the main vibration sources are the engine and road inputs through the wheels.Transfer path analysis quantifies the origins and parts contributing to structure borne noise in the passenger compartment.Rubber elements play here an important role as they are designed to reduce the transmissibility through a particular path.
Nevertheless, the designer is forced to find a relative optimum solution to ensure not only the desired vibration isolation from the different sources but also the adequate dynamic behaviour and manoeuvrability of the car.It is not only a question of solving the dilemma soft/stiff mounts, but having the correct damping factor.The use of isolators combining rubber and fluid cavities -hydromounts -is giving additional options to improve both vehicle confort and riding.
Elastomeric or rubber mounts have been used to isolate vehicle structure from undesired vibrations since the 1930s [41].Elastomeric mounts can be designed so that desired dynamic stiffness characteristics are achieved in all directions for proper vibration isolation.However, the dynamic behaviour of these mounts is such that both stiffness and damping increase slightly with frequency in the operation frequency range (usually not over 200 Hz in automotive applications), not providing extra damping or stiffness when necessary.
The problem has been partly solved via passive fluid mounts.Since in 1962 Richard Rasmussen patented a hydraulic mount for increasing damping [41], these type of devices have received greater attention, especially due to their highly tuneable characteristics.A passive fluid mount consists of two rubber cavities, full of fluid, that are connected through a channel (usually referred to as the inertia track), so that flow of fluid from one cavity to the other is permitted.The inertia track causes the dynamic stiffness of the mount to be more frequency dependent than in the case of simple rubber mounts, allowing the device to be tuned to have high damping at a desired frequency.Sometimes, the inertia track is not present and both chambers are connected through an orifice, resulting in a simpler system.Mainly in engine supporting applications, a decoupler may be incorporated to fluid mounts to provide amplitude dependent stiffness.
The present paper is focused on the mounts located in the suspension system of a car, where quite high amplitude vibrations are present, and thus, only variation with frequency of the dynamic stiffness of the mount is of interest.No decoupler is incorporated to the bushing.Such mounts are known as hydrobushings and are differentiated from the hydromounts, which are elements focused on isolating the engine from the structure, in those hydrobushings the amplitude dependence is neglected.The latter ones need to incorporate a displacement control via a decoupler, converting them into more complicated devices.
The behaviour of a hydrobushing is quite simple.At low frequencies fluid passes freely from one chamber to the other through the channel.The stiffness of the hydrobushing is approximately equal to that of the empty rubber mount (without fluid inside).As the frequency increases, the fluid travelling through the inertia track creates additional damping.Dynamic stiffness of the mount increases until a maximum value for the damping is reached, called the resonant frequency.Above the resonance frequency, the channel essentially closes off.At high frequencies, the stiffness of the mount is approximately equal to the sum of the elastomeric stiffness and the volumetric stiffness, higher than that of the empty elastomeric mount.Esentially, the dynamic stiffness of the hydrobushing is increasing from a base value (where it is similar to that of the empty mount) to a top value (augmented by the volumetric stiffness when no fluid is circulating through the inertia track) while the damping is showing its maximum close to a value called the resonant frequency.
The most widely used flexible chamber fluid mounts with inertia tracks have been thoroughly investigated by several researchers.Yu et al. [41], for example, provide a summary of different engine mounting systems, where fluid mounts are included, with references to papers where mathematical models are described.In the same way, Brach et al. [8] hold a thorough discussion on engine vibrations and the desired engine mount characteristics.They also present a complete literature survey of models prior to the work by Singh et al. [32], which gives a good description of hydraulic mount operation and characteristics.Min Lu [25] has conducted a thorough work on the characteristics of hydromounts, leading to interesting conclusions in regards to the influence of different parameters in the predicted dynamic stiffness of the bushing.
As far as development of models are concerned, the work by Flower [11] provides an excellent physical understanding of modelling techniques.Until recently, most attempts to model fluid mounts assumed linearity and were restricted to operating conditions (test conditions).Singh et al. [32], i.e., employed linear system theory to estimate dynamic stiffness of hydromounts.They carried out a linear mathematical analysis of hydromounts with long orifice and a decoupler using continuity and momentum equations.
Nevertheless, prediction of dynamic behaviour of such elements involves a variety of non-linear characteristics and is subjected to different excitations.Taylor [34] presented dynamic stiffness curves for a long orifice mount with and without a decoupler, experimentally obtained under different amplitudes and excitations.He showed that the analysis of the fluid flow through a long orifice in the hydromount, subjected to oscillatory excitation, is quite complex.After that, many authors, i.e. [9,14,16,19,20,37], have investigated and developed lumped non-linear models that predict the dynamic behaviour of hydromounts, most of them including the decoupler switching characteristics.Even more, i.e., Lee et al. [21] or Vahdati [38], have modelled the hydromounts using the bondgraph method, which is shown to be an efficient way to develop the multi-disciplinary model equations.
In this study, a different approach to the prediction of dynamic stiffness of hydrobushings is presented, combining Finite Element and CFD methods.The methodology proposed herein,compared to lumped models already mentioned, allows discussing about how the geometry of both, the mount and the channel, is affecting the dynamic stiffness of the system.Furthermore, it provides the designer with a tool that helps predicting, at the design stage, which is the effect that changes in the geometry of the bushing or any parameter of the inertia track (section, length, etc) as well as type of rubber or fluid will have in the final stiffness.The costly and time consuming trial and error process that usually takes place when physically having to produce a prototype for a specific application is so skipped.Bushings with inertia track and without decoupler are considered exclusively.The rubber structure of the mount is modelled in ABAQUS, the coupling between deformation and pressure in the cavities is considered through especial elements, and the flow of fluid is calculated in FLUENT.Results from CFD are incorporated into the FE model to obtain final dynamic stiffness of the hydrobushing with quite satisfactory results.Discrepancies with experimental measurements are mainly due to the simplifications assumed when creating the model.This approach is a novel and promising one in the prediction of the dynamic stiffness of hydrobushings.

Building the model of a hydrobushing. Methodology
A primary difficulty in modelling applications, in which cavities filled of fluid are present, is the coupling between the deformation of the structure and the pressure exerted by the fluid on it.The response of the structure depends not only on the external loads applied on it, but also on the pressure of the fluid over the cavity boundaries, which, in turn, is affected by the deformation of the structure.Furthermore, if different connected cavities exist within the application, the flow of fluid from one cavity to the other must be considered, as it also influences the overall response of the system.
The coupling between the deformation of the structure (due to external loads) and the pressure inside the cavities is provided in ABAQUS [1] through special elements called hydrostatic fluid elements.These elements completely cover the inner boundary of the cavities, sharing nodes at the cavity boundaries with structure elements.This way a membrane that defines the cavity is created, allowing the calculation of its volume and the change in the fluid pressure.
Figure 1 schematically shows the main idea of how these elements work.The applied exterior load deforms the structure, which makes the cavity filled of fluid also to deform.Hydrostatic fluid elements that enclose the cavity equal the normal force exerted on them varying the inner pressure P. A new degree of freedom is, thus, introduced in the problem: the pressure inside the cavity.This new degree of freedom is measured in a common node that is shared by all hydrostatic fluid elements associated to the cavity: the cavity reference node.This node is also used to calculate the volume of the cavity.
One of the most interesting capabilities of ABAQUS, as long as modelling fluid filled cavities is concerned, is the possibility of modelling the flow of fluid from one cavity to another in applications where connected cavities are present.Each cavity is modelled in the way described in the previous paragraphs, where a cavity reference node, measuring the pressure inside, is considered for each of them.However, ABAQUS does not model the connecting channel itself, neither the fluid flowing from one to the other: a special element is used -fluid link element -which accounts the mass flow rate of the fluid as a function of the pressure difference between cavities (see Fig. 2).That is, the only data needed to run the model are the mass flow rates (q) through the channel for different pressure values (∆P 12 = P 1 − P 2 ).The way this data is obtained is up to the user and it is further explained in Section 5. Once loads are applied and deformations calculated, ABAQUS estimates the pressure inside each cavity.The mass flow rate from one cavity to the other is known through the data introduced in the FE code.At each increment, the pressure in the cavities is actualised, in an iterative process till equilibrium is reached.
Relationship between mass flow rate (q) and pressure difference (∆P 12 ) may be specified either in a tabular form (mass flow rate, pressure difference) or introducing the viscous resistance coefficient (C V ) and the hydrodynamic resistance coefficient (C H ) according to Eq. ( 1): Nevertheless, in dynamic calculations (as the ones we are interested in when predicting the dynamic behaviour of a hydrobushing), it must be taken into account that ABAQUS assumes that the amplitude of the dynamic excitation is sufficiently small so that the problem may be treated as a lineal perturbation about a previous pressurized state.This fact supposes a linealisation of Eq. ( 1) in the point of the previous fluid state: where index "0" in Eq. ( 2) corresponds to that previous fluid state.If the system is not pressurized before the dynamic load is applied, data introduced for the fluid link element will be linealised around zero value of pressure, considering the value of the slope at the origin for C V and zero for C H .This fact may lead to significant errors in the calculus, being therefore much more effective the user to introduce the value of C V , which is considered to be most adequate, given the real range of pressures at which the system is working, as will be discussed in Section 6 when all this technique is applied to a real case.In summary, presented methodology allows splitting the problem into two calculations: estimation of the relationship between mass flow rate (q) and pressure difference (∆P 12 ) -which is suggested to be solved by the CFD code FLUENT, as further commented in Section 5 -and final prediction of the dynamic stiffness of the hydrobushing in ABAQUS, once ∆P 12 = f (q) data from the CFD calculation is provided to ABAQUS (according to Eqs (1) and ( 2)).
All this modelling technique is subjected to various simplifications that need to be taken into account.First of all, even though dynamic calculations may be carried out, the inertia of the fluid in the cavities is not considered in such calculations.It is up to the user to include this effect, if it is considered to be of interest, by distributing mass elements around the cavity boundaries.Secondly, there is also a limitation in the model when steady-state dynamic calculations are carried out, concerning the flow of fluid from one cavity to the other.The relation between the mass flow rate and the pressure difference Eq. ( 1) is always considered as linear in such calculations, even though that is not always the case.Finally, the inertia of the fluid in the channel is usually ignored when solving the inertia track in the CFD code to obtain the ∆P 12 = f (q) relationship, as steady boundary conditions are considered.This topic and its possible effect in the final dynamic stiffness results are discussed at the end of Section 6.
Furthermore, this methodology is very sensitive to a correct definition of the properties of both, the fluid and the deformable structure (natural rubber in the case we are interested in), as well as to the correct data set for the mass flow rate that takes place in the channels connecting cavities.Next sections are dedicated to these matters.

Fluid properties
The properties of the fluid must be specified to run the coupled fluid/structure problem.Not all properties of fluid need to be known; only density, viscosity and compressibility of the fluid are of interest concerning the FE code.If required, the density of the fluid may be considered to be dependent on temperature and pressure.
As for compressibility, the bulk or compressibility modulus of the fluid is required exclusively.It is assumed that this modulus is independent of fluid density changes, although it may also be specified as a function of temperature if needed.
Furthermore, ABAQUS allows making a distinction between hydraulic fluids (which are considered to be incompressible or nearly incompressible) and pneumatic fluids (which are modelled as an ideal gas).The first ones are the most common ones in applications such as hydrobushings, when usually the fluid is considered to be incompressible (liquid oils).Nevertheless, the user is responsible for confirming this point and the errors that may arise if it is not the case.
When different cavities filled of fluid are connected, it is supposed that fluid properties are the same for both fluids in each cavity.Since the connecting channel is not modelled (data corresponding to the mass flow rate as function of pressure difference is provided in a tabular form) properties for the fluid in each cavity are separately introduced in the model, being the user's responsibility to ensure that meaningful fluid transfers are taking place.

Modelling the dynamic behaviour of rubber
A correct definition of the properties of the structure is also critical if accurate predictions of dynamic stiffness of hydrobushings want to be obtained.Despite rubber being so common in this kind of applications, knowledge of material properties is quite often poor.Probably, because of its complex behaviour, as its dynamic mechanical properties are strongly dependent on temperature, frequency and strain amplitude.
Because rubbers are not perfectly elastic, the strain during cyclic deformation always lags slightly behind the stress.So, the application of a sinusoidal deformation will result in a sinusoidal force of the same frequency, but displaced in phase by an amount known as loss angle.It is very convenient to consider both, the elastic (in phase) response and the out of phase response, in terms of two moduli.The overall response can then be expressed as a complex modulus: In Eq. ( 3) G 1 is the in phase, storage modulus, and G 2 is the out of phase, loss modulus.The phase or loss angle δ is given by tgδ = G 2 /G 1 .Both are very dependent on frequency, temperature and dynamic strain amplitude.Therefore, to accurately define the dynamic behaviour of rubber, the complex shear modulus and its dependence on frequency, temperature and dynamic strain amplitude has to be investigated and measured.
A new test sample, based on the quadruple simple shear sample proposed by the standard ISO 1827 [18], is presented here to easily and accurately obtain the dynamic properties of natural rubber, properties that will be afterwards introduced in the FE code ABAQUS.

Measurements on simple shear samples to characterise rubber dynamic properties
Many techniques have been developed and verified for measuring the complex Young's modulus (or shear modulus) of viscoelastic and rubber-like materials.Payne [28], Ferry [10] and Mark [23], for example, propose and reference a large number of different type of tests to obtain the complex dynamic modulus.The set-up of some of these tests is considerably complex and some post-processing calculations have to be made to obtain the values really needed for the FE code.Some of them also have limitations in the frequency, temperature or amplitude range in which results can be obtained, not being adequate for the application of interest.The American Society for Testing and Materials (ASTM) has recently edited a guide [6] focused on describing several of the procedures for determining the dynamic properties of vulcanised rubber materials.As a guide, it is intended to provide descriptions of options available rather than to specify the use of any method in particular.Forced non-resonant vibration methods are the most widely used.Authors like Gent [15], Gade et al. [13], Tomlinson et al. [36] or Austrell [7] have tried them in different manners and with different test pieces with satisfactory results.Forced resonant vibration tests, compared to the non-resonant ones, have the disadvantage that there is not usually a continuously variable frequency range, although the range that can be covered is fairly wide.The only standard test that has been so far proposed to measure vibration damping properties of materials [5] falls into the forced resonant vibration tests.Nevertheless, it is not a useful test if the dynamic strain amplitude dependence needs to be investigated, as it is not controlled during the test procedure.Sim et al. [31], Pritz [30] or Tomlinson [35] also propose forced resonant methods to determine the complex Young's modulus of materials using either transmissibility functions or transfer function approach.Free vibration methods (like the one covered in ASTM D945-92 [4]) are not widely used, despite being the equipment needed comparatively simple, as often non-linear behaviour is observed and interpretation of results is then difficult.Moreover, these kind of tests are not able to cover wide range of frequencies unless test pieces of different sizes are used and the dynamic strain amplitude cannot be investigated.The possibility of using various tests depending on the frequency or strain amplitude range [24] is also of interest, though the need of different test set-ups and the cost of it could make it less attractive compared to using only a simple and unique test.
The authors have determined the complex shear modulus of natural rubber with great accuracy from simple shear dynamic tests at different frequencies and dynamic strain amplitudes through a forced non-resonant simple shear vibration test, using the sample shown in Fig. 3.Note that this test sample used for the dynamic measurements is half the one used in quasi-static tests, as proposed in the standard ISO 1827 [18].This new configuration overcomes the dynamic problems of the standardised quadruple simple shear sample due to inertia effect of the two metal parts that connect the four rubber blocks.Its main advantage, compared to other configurations, is that the test set up is exactly the same as for the simple shear quasi-static test.Thus, the same testing machine and clamping devices could be used for both tests.
This test sample has been used with success to investigate the dependence of rubber properties on frequency, temperature and strain amplitude.Some of the test results obtained are illustrated in Figs 4 to 8, comparing samples of different hardness (Shore A 40, 50, 60 and 70).In all cases, the material is natural rubber filled with two types of carbon black (SRF-N-772 and NEGROMEX N-330).In order to avoid Mullins' effect [26] samples have been mechanically conditioned and some initial cycles have been applied to them before starting the measurement.
It can be seen in Fig. 4 that G * module increases with hardness and slightly with frequency in the region of engineering interest.Anyway, till 1000 Hz [22] there seems to be moderate influence of frequency on the dynamic properties of natural rubber, a limit that exceeds the range of frequency of interest in the current application.It is also important to point out that a "jump" has been observed to exist between the quasi-static and the dynamic (even at low frequencies) of G (see Table 1).For unfilled rubbers the "jump" may not be very important (less than 25%), but as the amount of filler increases the "jump" becomes bigger.For heavily filled rubbers the ratio of dynamic complex modulus to the quasi-static one may be more than double.This "jump" is, of course, more noticeable  when amplitude is small.The reason for it is not yet understood, although it has also been observed by some other researchers (Sj öberg [33] for example).The damping (loss angle) is rather more sensitive to frequency than the modulus, but these variations follow a similar pattern, increasing damping as frequency rises.
Modulus of G * increases as temperature decreases, as shown in Fig. 5.This characteristic is common for rubbers of different hardness, even if additional results are not presented here.Natural rubbers containing a high proportion of reinforcing fillers show the largest variation of the dynamic modulus with temperature.Below 10 • C-0 • C there is a more rapid increase, as the material leaves the rubbery region and enters the transition region (Fig. 6).At about −60 • C, which is close to its glass transition temperature; natural rubber becomes glass-like and brittle.Similarly, loss angle of G* is almost constant in the rubbery region, increasing rapidly in the transition region up to a maximum value (which has not been achieved in the temperature range considered within these tests) (Fig. 6).
The dynamic modulus of rubber also decreases as the amplitude of the applied strain increases, as illustrated in Fig. 7.This effect is associated with the breakdown of interactions within the filler and between the filler and the rubber matrix and it is also known as Payne effect [29].Therefore, filled rubbers show a stronger dependence on this parameter than unfilled ones, which could be considered not dependent on amplitude.
However, the breaking of filler structure, described usually as friction behaviour, causes an increment in damping (loss angle of G * ) (Fig. 8).Nevertheless, the increase in loss factor or damping is not maintained in all range of amplitudes: the curve of modulus versus amplitude smoothes at large amplitudes.It seems to be related with remaining polymer chains and hydrodynamic effects [29].Finally, it must be taken into account that the properties of the material are very dependent on the rubber compound (especially type and content of reinforcing filler -carbon black -).It is important that samples and bushings modelled in the FE code are compounded from the same rubber material, so that useful material characteristics are obtained from tests in samples.

Introduction of dynamic properties of rubber in the FE code ABAQUS
It is not easy to define the complex modulus of rubber material in ABAQUS, being the use of the viscoelasticity the only way to do it.However, difficulties arise as ABAQUS uses a linear theory to define the time-dependent or viscoelastic properties of rubber [2].It is assumed that material properties are only dependent on time (or frequency),  but not on the dynamic strain amplitude, which is not completely true, especially for filled rubbers, as shown by experimental tests performed on simple shear samples (Section 4.1).Therefore, the dependence of dynamic stiffness on the dynamic strain amplitude cannot be directly introduced in the code.This fact may become critical in cases where the amplitude effect is relevant.
At present work is being done in that direction and some solutions on how to implement the amplitude dependence on a commercial FE code are being investigated, though they are out of the scope of the present paper.The methodology followed herein consists in introducing the frequency dependent properties of the rubber material at a specific dynamic strain amplitude: the working dynamic strain amplitude of the hydrobushing, that has been estimated.Simulations of the behaviour of the rubber mount at very low and high frequencies have validated the usefulness of the data considered, as will be illustrated in Section 6 when the practical case is described.
Frequency dependent properties, for one dynamic strain value, are introduced in ABAQUS in a tabular form: ωRe(g), ω Im(g), ω Re(k), ω Im(k), freq (Hz) Where: ω Re(g G ∞ and K ∞ are the long-term shear and bulk moduli, obtained from quasi-static simple tests [18] at low speed.If rubber is considered to be fully incompressible (which is almost always the case, being its bulk modulus about one thousand times larger than the shear modulus), terms involving the bulk modulus are ignored (zero value is introduced).

Characterisation of the channel and fluid flow
To model the fluid-structure interaction in the hydrobushing, it is necessary to complete the mathematical model provided by ABAQUS with the relationship between the pressure difference or pressure drop in the cavities (∆P 12 ) and the flow-rate in the channel (q) which links them.It is possible to obtain this relationship experimentally or by solving a mathematical model of fluid flow within the channel.In this paper the latter option has been considered since experimental methods are laborious to set up and economically expensive.
Three channels have been considered in the study (Fig. 9, channels F, G and H), all of them having a constant rectangular cross-section.Figure 9 schematically shows the different configurations selected for the current application: channel F and G (left) only differ in the size of the section, the configuration being the same for both, while the configuration of channel H (right) is different from the other two, as well as the size of its cross-section.
In order to solve the pressure drop in the cavities as a function of the mass flow rate three different approaches (i.e. three different mathematical models) have been used.These are shown below in order of complexity.
The first and easiest approach is to express the relationship between ∆P 12 and q by considering the channel as a straight conduit and the flow inside it as incompressible, steady, fully-developed and in the laminar regime, when Re < 2300.It is calculated using the hydraulic diameter of the section.This relationship can be taken from the bibliography [40]: Where ∆P 12 , µ, L AB , a and b and η are the pressure drop, the fluid viscosity, the length of channel, the half-lengths of the sides of the channel and a known function of b/a, respectively.Using Eq. ( 4), the simulated values of dynamic stiffness of hydrobushing do not fit well the experimental data.This is because it has been assumed that the channel is a straight conduit, which is not a valid supposition since drops in pressure, due to changes in the flow direction and the curvature of the channel, have not been taken into account.
Therefore, in the second approach, the effects on drops in pressure due to changes in the direction of the flow inside the channel have been considered.this reason the channel has been considered as being made up of a set of special pieces (elbows and doubled elbows) linked by means of curved conduits (see Fig. 10 for a schematic representation of channel F and the different parts it can be divided into).The latter approximation is a common practice in the analysis of hydraulic networks, where the pressure drop between two points of the network is equal to the sum of energy losses suffered by the fluid as it passes through the special pieces and pipes which connect the two points.In this way the pressure drop between extremes of the channel is expressed as the sum of the pressure drops (i.e.: energy losses) of fluid as it passes through the different parts in which the channel has been divided.In a type F channel case (Fig. 10), the total pressure drop can be written as: Where ∆P 1A , ∆P BC and ∆P D2 are the pressure drops in these special pieces in which the channel has been divided, and where the flow changes its direction suddenly; and ∆P AB and ∆P CD are those produced in the curved conduits which connect the special pieces.
Although dynamic stiffness of the hydrobushing obtained by simulation using the second approach is better than the result obtained with the first one, important discrepancies still exist with respect to the results obtained experimentally.This second approach provides better results in hydraulic networks because the losses produced by special pieces are less important than those due to friction in the pipes.Moreover, in hydraulic networks special pieces are usually separated by long pipes, along which it is possible to consider the flow as fully-developed, allowing the energy loss in each special piece to be handled independently and the total energy loss to be obtained as a sum of the individual losses.However, since the channels considered in this paper are very short and the flow does not fully develop in any significant part of their length, it is incorrect when evaluating the pressure drop to consider the pieces of the channel independently.In addition, another problem associated with this approach is the fact that the most widely available bibliographical data on energy losses in special pieces is related to turbulent rather than laminar regime [17].
Finally, Computational Fluid Dynamics (CFD) techniques have been applied to obtain the pressure drop vs. flow rate relationship of inertia tracks of hydrobushings.CFD includes the modelling and numerical techniques to simulate phenomena in which, fluid flow, heat transfer and chemical reactions are present.Nowadays, as a consequence of increasing computers speed and storage capabilities, CFD is becoming increasingly sophisticated, being applied more and more in research and in industry [39].
The mathematical model of fluid flow consists of governing equations and mathematical definition of the boundary and initial conditions.The fluid flow into the channel is governed by well-known Navier-Stokes equations: where u i (x, t), p(x, t), ρ and µ are, respectively, the velocity component of the fluid in the spatial i-direction, the pressure, the density and the molecular viscosity of the fluid.These equations express the law of mass conservation and Newton's second law.The flow into the channel is considered incompressible and therefore ρ and µ are considered constant.Initial and boundary conditions must be imposed to complete the mathematical model of flow model.Non-slip boundary conditions (u i = 0) have been imposed at the walls of the channel.The respective values of pressure have been imposed at the ends of the channel.Since channel flow is steady, it is not necessary to establish initial conditions.
In order to solve the continuous mathematical model of the flow numerically, it is necessary to turn this model from a continuous to a discrete model that consists of a system of algebraic equations.Differential equations, flow domain and boundary conditions, all of which form the continuous mathematical model, are merged and transformed into a discrete mathematical model by applying a discretization technique.The unknowns of this system of equations correspond to the values of flow variables (components of velocity and pressure) at a discrete group of positions within the spacial domain and at particular instants in time.
The discretization technique is initially applied to the spatial and time domains, dividing both into discrete subdivisions.If the flow is steady, it is not necessary to discretize the time domain.The spatial flow domain is divided into a discrete number of small and non-overlapping volumes, named elements or cells, with established shape (hexahedral, tetrahedral, pyramid or wedge).This process is known as meshing, and is a key process in all CDF analysis.
Finally, the discrete mathematical model is obtained from the mesh and the discretized time domain, the partial differential equations, and the boundary and initial conditions by applying a further discretization technique.Several techniques exist to obtain this discrete mathematical model; the most well known are the Finite Differences Method, the Finite Elements Method and the Finite Volume Method [27].
Fluent v5.5 [12], the commercial CFD code used in this work to solve the mathematical model of the flow, implements the Finite Volume Method.The latter discretization method begins by integrating Eq. ( 6) about the volume of each cell of the mesh.These integral equations are discretized leading to a system of equations that in a steady case can be written as: where φ P and φ nb are the unknown values of a flow variable (velocity component or pressure) at the centre of cell P and at the centre of cell P's neighbouring cells, respectively.S φ is the independent term.Second order schemes have been used to evaluate diffusive and convective terms in order to discretize integral equations in the Finite Volume Method.The SIMPLE [27] method has been used to couple the calculation of pressure and velocity since the flow is considered incompressible.
Once the discrete mathematical model of flow has been solved numerically for several pressure drops in the channel, the relationship between pressure drop and flow-rate can be obtained (Fig. 11).The simulated results of hydrobushing dynamic stiffness using the aforementioned relationship equate well with experimental data (see Section 6).
This CFD approach provides the best agreement between simulated and experimental data since it is capable of exploiting a more complicate mathematical model of the flow that takes into account all of its main features.However, this approach is computationally more expensive than those specified earlier and when the hydrobushing has different geometric characteristics it is necessary to solve the mathematical model of the flow again to obtain a new pressure drop/flow rate relationship.

Application to a real case
All this modelling technique has been applied to a practical case, intending to predict the dynamic stiffness of an axial hydrobushing.The volume of cavities and the pressure inside them are modified when it is axially loaded, the oil flowing from one cavity to the other in such circumstances.Three different channels have been considered connecting both cavities (Fig. 9).
The mesh generated in ABAQUS may be seen in Fig. 12.Only the rubber part is modelled, as the rest of the elements (inner and outer cylinders, channel and closing elements) are considered to be rigid compared to rubber.The  boundaries of both cavities are covered with hydrostratic fluid elements (to ensure coupling between the deformation of the structure and the fluid) and the cavity reference nodes are connected through the fluid link element (FLINK) where the information on the fluid flow is introduced (relation between mass flow rate and pressure drop).This relation has been obtained from modelling the fluid in FLUENT, in the way described in Section 5, and with the results summarised in Fig. 11.
As for the boundary conditions, the outer cylinder, channel and closing elements are fixed and the inner one is subjected to a harmonic excitation of 0,4 mm amplitude.The material is natural rubber of hardness Shore A 60.As it is a filled rubber, the amplitude dependence may not be neglected and should be considered within the calculations.Nevertheless, it cannot be directly introduced in the code.
It is expected to obtain accurate results from the simulations if properties of rubber at the correct amplitude value are introduced.The deformation of the hydrobushing is close to simple shear.Therefore, the amplitude of the dynamic strain has been calculated as if it was simple shear and measurements on samples have been performed for amplitude values close to the result obtained.Afterwards, the real properties have been estimated by fitting the numerical results to experimental ones in two special cases: 1) when the cavities of the hydrobushing are empty (instead of full of fluid); and 2) when the cavities are full but they are not connected through a channel.In both cases the dynamic stiffness of the hydrobushing is only dependent on the properties of the rubber and the fluid (in the latter one), but not on the flow of fluid through the channel.The relationship between pressure drop and flow-rate, either in tabular form or through the coefficients C V and C H (see Section 2 and Eq. ( 1)), is introduced in ABAQUS through the definition of FLINK element.Nevertheless, the user must not forget that when dynamic calculations are involved, ABAQUS assumes that the amplitude of the dynamic excitation is sufficiently small so that the problem may be treated as a lineal perturbation about a previous pressurized state.The immediate effect of this assumption is that the curves are linealised, according to Eq. ( 2), in the point of the previous fluid state, assigning zero value to C H and the value of the slope of the curve in the previous pressure state to C V .In order to control which value is assigned to C V it is recommended to fit the curves to straight lines, minimising the error, and directly introduce the value of C V in the FE code.This has been the procedure adopted here, with satisfactory results.
The results from the numerical calculations compared to the measurements to validate the models can be seen in Figs 13 and 14, in the case that channel F or G are considered.Results for the remaining case (channel H) are not presented, though similar conclusions may be drawn from all of them.All cases show a good correlation between the measured and predicted module of the stiffness (N/mm) with errors less than 10% (which is considered acceptable among industrial manufacturers).The stiffness of the hydrobushing falls somewhere between two asymptotic values: at low frequencies, the stiffness value is close to that of the hydrobushing working without fluid; on the other hand, at high frequencies, the value of the stiffness should be close to what happens when both cavities are not connected at all, as if there was no fluid flow through the channel.The model shows this effect with accuracy, especially in the case of channel G (see Fig. 14).
Nevertheless, damping is not so well predicted.The maximum value of the loss angle falls into the acceptable error (± 2.5 • ), although it is not achieved at the same frequency.The maximum of the predicted phase is displaced 4-6 Hz from the measured one.It is thought that this fact could be due to considering steady boundary conditions, instead of harmonic ones, when calculating the channels in FLUENT, as the inertia of the flowing fluid is not taken into account.Several authors have discussed the effect of considering the oscillation of the fluid in the results.Min Lu [25] demonstrates that the maximum loss factor occurs near the inertia track fluid resonant frequency.Ahmed et al. [3] have performed different analysis, with and without the consideration of fluid oscillation, being discrepancies attributed to that effect.
In fact, channel F has been solved in FLUENT for harmonic boundary conditions (fluid oscillation) and it has been seen that the pressure drop (∆P 12 ) and the mass flow rate (q) are not exactly in phase as it is assumed when stationary boundary conditions are applied.This effect is much more important as frequency increases (Figs 15 and  16).Nevertheless, this out of phase effect cannot be so far included in the FE code (it is not considered so in Eq. ( 2)).

Concluding remarks
A novel approach to the prediction of dynamic stiffness of hydrobushings has been presented by the authors, combining Finite Element and CFD methods.The dynamic properties of the rubber have been determined using a test sample derived from the standardised static sample, which is not suitable for dynamic tests.Test sample suggested in this study is half the one used in quasi-static tests, which makes it very easy to test in the same testing machine and with the same clamping devices.Thus, a simple shear non-resonant test is suggested.
Concerning the flow of fluid in the inertia track, the use of CFD provides better results than any other procedure, as it is capable of solving a more complicated model of the flow, taking into account all of its main features.Two are the main assumptions of the model.First of all, and even if CFD results show it is not so, the relationship between pressure drop and flow rate must be considered linear when introduced in the FE code.This problem has been overcome by linealising the ∆P 12 = f (q) curves in the operation point.Secondly, the fluid flow is considered steady when solving it in FLUENT, so that the pressure drop and the flow rate are in phase.That is, the inertia of the fluid in the inertia track is not considered, damping being provided only by friction in the channel.As it has been stated by different researches and confirmed by the authors, the inertia of the fluid is affecting the resonant frequency at which the maximum damping appears.
The dynamic stiffness of a real hydrobushing has been predicted making use of the methodology proposed with satisfactory results, despite of the simplifications assumed in the model.The stiffness is well predicted and so is the damping, although the frequency at which its maximum value is achieved is underestimated by 4-6 Hz.

Fig. 1 .
Fig. 1.Schematic representation of the coupling between the deformation of the structure and the pressure exerted by the fluid on the boundaries of the cavity.

Fig. 2 .
Fig. 2. Representation of the fluid flow between two cavities connected by a channel.

Fig. 3 .
Fig. 3. Simple shear test sample proposed for characterising dynamic behaviour of rubber.

Fig. 9 .
Fig. 9. Schematic representation of the three channels (F, G and H) considered in the study and their sections (in mm 2 ).

Fig. 10 .
Fig. 10.Schematic representation of channel F and the parts it can be divided into.

Fig. 11 .
Fig. 11.Pressure drop as function of mass flow rate for the three channels (F, G and H) considered within the practical case.

Fig. 13 .
Fig. 13.Numerical results vs. experimental tests for the hydrobushing with channel F.

Fig. 14 .
Fig. 14.Numerical results vs. experimental tests for the hydrobushing with channel G.

Fig. 15 .
Fig. 15.Sinusoidal pressure applied to boundaries of Channel F in FLUENT and sinusoidal mass flow rate obtained.Frequency: 2 Hz.

Fig. 16 .
Fig. 16.Sinusoidal pressure applied to boundaries of Channel F in FLUENT and sinusoidal mass flow rate obtained.Frequency: 50 Hz.

Table 1
Comparison between the quasi-static and dynamic G * modulus (MPa) at room temperature and 0,05 mm amplitude Shore A |G Quasi-static values have been obtained by standardised quasistatic tests and considering linear relationship between stress and strain, which is not completely true due to the non-linearity of material.