Tuned inerter dampers with linear hysteretic damping

This paper explores the influence of linear hysteretic damping on the performance of passive tuned‐inerter devices. An inerter is a device that produces a force proportional to the relative acceleration across its two terminals; devices incorporating inerters have received widespread attention in the earthquake engineering community, because they offer the ability to improve the seismic response of structures. However, the majority of this research has assumed that the damping components within the tuned‐inerter device exhibit viscous, rather than hysteretic, damping. This restriction imposes an essential question on how the hysteretic damping model will change the performance of the device compared with the viscous damping model. It is shown that the response of viscous and hysteretic inerter systems have significant differences in displacement amplitude due to the frequency dependency of the damping. Therefore, a new formulation for obtaining the optimum loss factor of the hysteretic damping in the inerter system is proposed. Next, the challenges associated with accurately predicting the time‐response of a hysteretically damped system are discussed. A numerical time‐integration method is extended to address these challenges, using a new formulation that has the benefit of being broadly applicable to multidegree‐of‐freedom hysteretic linear systems and nonstationary random signals. The results show that the earthquake responses from the hysteretic damping model can differ significantly from the ones obtained via the viscous model.


INTRODUCTION
In 1997, Okumura Atsushi 1 patented the so-called gyro-mass device, which operated in parallel with a spring and a damper. This gyro-mass, consisting of two disks acting as inertia elements, rotated on a geared rod in such a way that the force generated by the gyro-mass, F, is equivalent to the relative acceleration between its two terminals, x A and x B . This relation can be expressed as Here, J is the moment inertia of the disks, and r d is the distance from the centre of the disks to the point where the rod is attached (otherwise called the radius of gyration). The parameter b d is called the inertance, which is measured in kilogram, and can in fact be achieved by several mechanisms, such as rotational motion of a flywheel 2 or fluid flow. 3 Following this initial concept, several other types of gyro-mass or inerter devices have been proposed. Smith 4 developed an inerter, consisting of flywheel, rack, pinions and gears. The force generated by these types of inerter can be expressed in the same way as Equation (1). Recently, many studies on the realisation of this inerter concept have been conducted based on three main mechanisms: mechanical geared, ball screw and fluid-based inerters. [5][6][7][8][9] The application of the inerter for vibration suppression in civil engineering structures has been studied by proposing several configurations of three basic elements: inerter, viscous damper and spring (or equivalent stiffness element). One of the first inerter-damper configurations to be introduced was the tuned-viscous-mass-damper (TVMD), described in detail in Ikago et al. 2 The damper consists of a parallel-connected ball-screw inerter-viscous damper in series with a spring element. Of particular note is that a device of this type has been put into service in a real structure. 10 The generalisation of this device layout is called a parallel-connected-viscous-inerter-damper (PVID). For example, Pan et al. 11 derived the optimum solution for the PVID in a damped single-degree-of-freedom (SDOF) structure with a numerical optimisation technique, using allowable parametric bounds. More recently, the feasibility of using the PVID concept for vibration suppression of civil structures using a fluid inerter has been studied. 12,13 Lazar et al. 14 proposed a tuned-inerter-damper (TID) for vibration suppression of civil engineering structures. The TID consists of a parallel connected spring and viscous damper in series with an inerter. This arrangement is similar to that of a tuned-mass-damper (TMD) with the mass element being replaced by an inerter. The TID has been proven to have similar performance to the TMD for the targeted resonance, but it also offers additional benefits compared with the TMD. As is well known, the TMD typically requires a large amount of mass to provide a large mass ratio and to achieve optimum performance. For the same mass ratio, the actual mass of the TID is far less than that of the TMD because the inerter element can produce inertance several times higher than its physical mass. 4 Additionally, it has been also shown that the TID achieves its best performance when placed at the base of the main structure, which is contrary to the TMD which must be placed at the top of the structure, where the maximum displacements occur. To obtain its optimum damping parameter, Lazar et al. 14 proposed a tuning rule that is based on the fixed-point theory by Den Hartog 15 for the TMD. In a similar manner, Hu et al. 16 studied five configurations employing an inerter including the TID layout and gave a tuning example for each of them when attached to an undamped SDOF structure. More recently, the TID has been studied for vibration suppression of cables 17 and for enhancing the performance of base-isolated structures subjected to earthquake base-excitation. 18, 19 Marian and Giaralis 20 proposed an alternative device called a tuned-mass-damper-inerter (TMDI). The device consists of a traditional TMD with an inerter attached to the mass element. The advantage of the device is its capability of reducing the total mass of the TMD. The concept of changing the inertance by using a control system to provide larger mass ratio has also been investigated. 21 The frequency of the TMD can be easily tuned to fit the frequency of the excitations. Optimal design of the TMDI for dynamic vibration suppression by assuming white noise process as base input was studied by Pietrosanti et al. 22 The TMDI has also been studied for application to base-isolated structures. 23,24 Following the initial analytical studies on the application of inerter-based dampers for vibration suppression of civil structures, a selection of experimental studies have been conducted in order to validate those well established concepts. Brzeski et al. 25 conducted an experimental study of the TMDI on an SDOF system. An experimental study of a PVID concept using a mechanical geared inerter can be found in Hessabi and Mercan. 26 More recently, Hessabi 27 conducted experimental work testing a mechanical geared inerter-based damper for vibration suppression of a civil structure using the real-time hybrid simulation method.
Despite the successful studies of the inerter-based damper systems both analytically and experimentally, many research questions remain. In particular, models such as Equation (1) are highly idealised, and as the development of these concepts progress, it is important to have models that can capture more physically realistic cases. For example, some studies have suggested the importance of considering the nonlinearities of both the damping and the inerter devices. Buelga et al., 28 amongst others, 29-31 studied the effect of the ball-screw inerter nonlinearity on the performance of the structure to which the inerter is attached. This paper proposes the use of the hysteretic damping as a new alternative to viscous damping in the idealised inerter-based damper model, and in particular, we will focus on the tuned-inerter-hysteretic-damper (TIhD) and the tuned-mass-hysteretic-damper-inerter (TMhDI). Hysteretic damping has the benefit of smaller response amplitude amplification at higher vibration frequencies due to its energy dissipation being frequency independent. Systems of this type have been studied extensively. [32][33][34][35] The hysteretic damping model also better describes the behaviour of certain materials. 36 Mathematically, one of the most frequently used ways of modelling hysteretic damping is as part of a complex stiffness function, which has the problem of noncausality-an issue that shall be discussed in detail below.
Equivalent viscous damping discussed in Chopra 37 has been widely used to model the loss factor or damping in the complex stiffness case. However, it is well known that for both single-degree-of-freedom (SDOF) and multi-degree-of-freedom (MDOF) systems, that is, civil structures, the use of equivalent viscous damping may only be valid for a particular frequency range near the frequency of interest, which is usually the resonance frequency. 38 In addition, concerning a nonharmonic excitation case, such as earthquake inputs, it is essential to evaluate the structural response in a time-domain basis rather than a frequency-domain basis. Therefore, in the present work, we investigate the hysteretic damping by developing a new method in a time domain.
The noncausal nature of the complex stiffness function leads to the occurrence of unstable poles in the time domain analysis of the system. As a result, implementing this type of model is not as straightforward as the viscous model. The idea of a causal hysteretic element has been developed as an alternative, 39 although this too has significant difficulties with computation. 40 To deal with the unstable poles issue, several methods have been proposed in the case of SDOF systems.
The use of analytic signals, along with the Hilbert transform and a time reversal technique, is one of the first methods introduced to solve such systems in the time domain; see Inaudi and Makris. 41 In this paper, the method of Inaudi and Makris 41 is modified in order to solve the time domain response of the structures with TIhD and TMhDI. This modification takes advantage of some key simplifying assumptions that are justified for the types of system under consideration. As a result, the proposed method makes it possible to solve the problem using standard numerical integration algorithms such as the Runga-Kutta family of solvers that are available within MATLAB, rather than the more specialist methods required previously. A discussion about the method is presented in Section 3.
A recent investigation in Málaga-Chuquitaype et al. 42 presented material modelling of an inerter damping system by means of finite elements. Both SDOF and MDOF systems were investigated, and responses of harmonic excitations and earthquake cases were evaluated. However, for a preliminary design stage, numerical simulation via finite element analysis could be computational costly. Hence, the present numerical method via Hilbert transform and time-reversal technique is proposed and implemented to support analysis on TIhD and TMhDI designs. Both the TIhD and TMhDI use a hysteretic damping element instead of a viscous damper. The motivation for using this damping model is to try and more accurately represent the case when a material damping (e.g., rubberised component) is present in the inerter-based device. When the mass of the inerter is also included in the model, then the TMhDI model becomes more appropriate. It should be also noted that in this work, both TIhD/TMhDI and the host structure are assumed linear. This is an important assumption of the models used in this paper.
The TMhDI has a small secondary mass between the inerter and the complex damping elements, analogous to the TMDI concept. In this paper, the secondary mass is considered to be approximately 5% of the inertance, to represent the mass of the inerter device only. The idea is that this TMhDI will capture both the effects of the inerter mass and the hysteretic damping and therefore offer a potentially more physically realistic model than the idealised concepts of the TID and the TMDI.

ANALYTICAL MODELLING AND OPTIMISATION
A generalised model of a n-DOF lumped mass system is shown in Figure 1. The structure is separated into three parts: bottom storey, i = 1; middle storeys, i th , where i ∈ [2 ∶ n − 1]; and top storey, i = n. In this study, the performance comparison of four different inertial damper systems when placed at the bottom story level was examined. As a well established concept, the TID was used as a benchmark to the performance of the other devices. The TMDI was studied to see the effect of the small additional mass element in the TID. Then the TIhD and TMhDI were studied to gain insight into the effect of hysteretic damping.
The equations of motion of the n-DOF structure as shown in Figure 1 in absolute coordinates can be written as where m i and k i−1,i , i ∈ [1 ∶ n] represent the mass and stiffness between storeys i − 1 and i; X i represents the Laplace transform of the ith storey displacement, when i = 0, then X 0 = R, which represents the Laplace transform of the base displacement; s represents the Laplace transform variable and F 1,0 represents the force transferred to the structure by the inerter damper in the Laplace domain given in Table 1. The inerter damper systems in Table 1 have the additional parameters b d , m d , k d , c d and which are the device inertance, mass, stiffness, viscous damping and hysteretic damping loss factor, respectively, and = √ −1.

Optimisation in SDOF structures
In this section, the tuning procedures to obtain the optimum parameters of the hysteretic inerter-based devices in an SDOF structure will be derived based initially on the traditional concept using viscous damping, which is in this case being replaced by complex stiffness term k d (1 + j ) representing the material damping. Here, k d is the stiffness of the inerter dampers and is the loss factor of the hysteretic damping element. It should be noted that this model is a noncausal model, which means it is physically not realisable. However, this model has been widely accepted in analysis 43 to accurately represent a class of nonlinear damping, 44 as well as the phenomena of energy dissipation in a variety of materials such as rubber and viscoelastic polymers. [45][46][47][48] As a result, there is a limitation that the model proposed in this work is currently only applicable to linear structures. Although the equivalent viscous damping is often directly used to obtain the approximate equal response between viscous and hysteretic damping, in this study, we found that the result is different after the optimisation process. This limitation has also been reported in Wong 49 for the hysteretic dynamic vibration absorber (HDVA). Therefore, it is important to derive a formulation for obtaining the optimum parameters of an inerter device with hysteretic damping.

Optimisation of the TIhD
Consider now a TIhD attached to a 1-DOF structure as shown in Figure 1 for n = 1. The transfer function X∕R can be derived analytically and written in the frequency domain as where is the loss factor of the hysteretic damping, = b d ∕m 1 is the inertance-to-mass ratio, q = ∕ n is the frequency ratio, = k 0,1 ∕k d is the stiffness ratio, is the forcing frequency (assuming a sine wave input) and n = √ k 0,1 ∕m 1 . The approximate optimum parameters of the TIhD can be analytically obtained based on fixed-point theory via an algebraic solution. This technique has been derived by Hu et al. 16 for the TID optimisation. Following the same procedure, and by replacing the viscous damping ratio with the loss factor , one obtains the optimum for the TIhD as P and Q are given by where The frequency ratio at the first two fixed-points P and Q are the solution of the following equation: The optimum stiffness ratio ( opt ) can be obtained by using either where q R is the frequency ratio at the third fixed-point R given by Figure 2A,B shows the difference between the optimised TID and TIhD using fixed-point theory. It can be seen that for viscous damping, Figure 2A, the underlying peaks of P and Q are quite low in amplitude, and therefore, the peaks of the combined optimal curve are very close to the fixed points P and Q. However, in the hysteretic damping case, the underlying peaks of P and Q are now quite high in amplitude, and therefore, the peaks of the combined optimal curve are much further away from the fixed points P and Q. It is therefore suggested that the optimised TIhD will require a further fine tuning process to make both peaks closer to equal amplitude. The effect of increasing is also obvious. It can be seen when is increased from 0.3 to 0.9, the TIhD shows a potential benefit in the higher frequency region where it has smaller amplification effect. Table 2 presents the obtained optimum values of and after the optimisation process. As previously mentioned, equivalent viscous damping is often used to obtain an approximation for the response of a hysteretically damped system. However, this kind of approximation will not persist after the optimisation process. In this case, the equivalent viscous damping relationship = 2 given in Chopra 37 cannot be used for the TIhD. The approximate relationship can instead be expressed as = 2 . It should be noted that this expression does not give the exact optimum solution-a point that was also noted by Wong 49 for the HDVA system. Therefore, it is important to derive the optimum of the TIhD separately from the derivation of of the TID.

Optimisation of the TMhDI
TMhDI is basically a TIhD with a small secondary mass element m d between the inerter b d and the complex-stiffness element as illustrated in Figure 1C. The TMhDI model can potentially be used to represent an even more realistic device than the TIhD device. In this analysis, the inertance is designed to be dominant, where m d is assumed to be just 5% of the inertance to represent the mass of the inertial device. Therefore, the tuning of the TMhDI follows the TIhD tuning rule based on fixed point analysis presented in the previous subsection. Figure 3 shows the performance comparison in the frequency domain between the four considered inertial damper systems in an SDOF structure subjected to harmonic base displacement. Although the response at low frequencies is close for all four devices, it is clear from Figure 3A that the devices with hysteretic damping give significantly lower amplitude responses at higher frequencies. In this higher frequency range, the hysteretic damping also changes the phase angle as can be seen in Figure 3B.

Optimisation in MDOF structures
First, we comment on the optimal location of the TMDI in MDOF structures. It is well known that the optimal location for a TMD is at the top of a MDOF structure and that the TID is optimum at the base of the structure. 14 The passive control High was chosen to clearly see the difference between the structural response with hysteretic and viscous damping above the resonance frequency forces in the Laplace domain can be written as follows when a TMDI is mounted at the ith storey level.
where X i represents the Laplace transform of the displacement of mass m i and b d is the inertance. The equation of motion for the y-DOF TMDI system can be written as Substituting Equation (12) into Equation (11a) gives where We now discuss the optimum placement of the TMDI in a n-DOF structure as shown in Figure 1. In matrix form, the equation of motion of the system can be written as where M and K are the mass and stiffness matrices and Z = X i − R represents the vector of relative storey displacement.
The above equation can also be written in form of modal matrices where is the eigenvector matrix and  = Φ T MΦ and  = Φ T KΦ are the modal mass and stiffness matrices, respectively. Now, we consider the cases where the TMDI is mounted at either the bottom or top storey level. Assuming that only the first vibration mode is significant, the following transfer functions can be obtained:

TMDI at bottom storey level
TMDI at top storey level It is obvious that when the mass element of the TMDI is small or close to zero, the above equations will be identical to the TID transfer functions given in Lazar et al. 14 Figure 4A shows how the optimum location of the TMDI changes between top and bottom story level by increasing m d and decreasing b d at the same time using Equations (17) and (18). It can be seen that increasing m d will make the TMDI at the bottom story less effective. Conversely, decreasing m d makes its performance better when placed on the bottom story. This aligns with our assumption of making the m d only 5% of the device inertance, to represent the mass of the inerter device.
Furthermore, the optimum parameters of the TID in a MDOF structure were obtained by using a method proposed by Lazar et al. 14 Figures 2A,B and 4A have already demonstrated that the tuning procedure of the TID can be used for estimating the optimum parameters of the TIhD, TMDI and TMhDI. However, as previously noted, additional fine tuning is required to get the actual optimum values for these three systems, and this is also the case for MDOF systems.

Concluding Remarks
From the discussion above, it can be concluded that when the damping in the tuned-inerter devices (TID and TMDI) exhibits linear hysteretic (TIhD and TMhDI) rather than viscous, the behaviour of the devices will be different, especially at the higher frequency of excitations. This is obvious from the frequency response analysis of the systems. However, obtaining the structural response in the time domain will be challenging due to the inherent occurrence of unstable poles. Hence, a robust time domain technique is required to capture the hysteretic behaviour rather than using an equivalent viscous damping method.
The following section will describe this new method for obtaining the time domain response of structures with linear hysteretic damping. Furthermore, some validations of the method are presented.

TIME DOMAIN MODELLING
In this section, for the first time in the time domain analysis, the hysteretic-based inerter systems are evaluated by means of analytic signals and time reversal technique. The method used here is based on that proposed by Inaudi and Makris 41 in an extended and modified format that allows for the evaluation of the TIhD and TMhDI systems. The hysteretic damping is modelled in the form of a complex stiffness element of the inerter device. Hilbert functions are used to represent the analytic signals of the structural responses.
It should be noted that this time domain method does not change the noncausality of the linear hysteretic damping. It has been a long standing problem that complex stiffness model is a noncausal model and therefore very difficult to deal with. However, this model has been widely used 43 and has been proven to be very accurate and practical for modelling several engineering applications, for example, linearisation of nonlinear damping. 41 Reference to the causal linear hysteretic damping model can be found in Makris. 39 Furthermore, the causality of the inerter-based damper devices has been extensively discussed in Makris 50 and Makris (2018). 51 Specifically, in Makris, 50 the concept of inertoviscoelasticity is introduced. This is a form of model built from combinations of springs, inerters and viscous damping elements. These type of mechanical networks were discussed in detail in that paper and have been proven to be causal. In Makris, 51 the time-response function of three-element (spring, dashpot, inerter) networks with different arrangements were considered. The method was adopted from Makris 50 to ensure the causality of the mechanical systems in the time domain.
In this paper, the focus is not to resolve the noncausality of the linear hysteretic damping (complex stiffness) model, but to give a better representation of that model in the time domain in addition to the frequency domain. The time reversal technique is adopted and extended from Inaudi and Makris 41 to resolve the unstable poles from the differential equations. A new efficient numerical method is proposed, which is also applicable for MDOF systems and nonstationary input signals.
However, there are some important nuances to consider. Specifically, Bae et al. 52 showed that the Hilbert transform of the external impulse signal using the discrete Fourier transform (FT) and inverse Fourier transform (IFT) produces a totally different result compared with the one obtained from its exact formula. Not only this, Bae et al. 53 also showed that the Hilbert transform of a signal may produce a nonzero value at the beginning of the signal although its original real signal has zero value at time zero. To solve these problems, Bae et al. 52,53 introduced a technique called discrete convolution Hilbert transform (DCHT) and extended-time-duration using a Newton-Raphson iteration method. In the DCHT technique, the external amplitude of the force signal is divided into a finite number of rectangular signals. The Hilbert transform of the force signal can then be obtained from the Hilbert transform of each of the rectangular signals. The application of this method to a free and transient response of a hysteretic damping system was presented in Bae et al. 54,55 A system with hysteretic damping was studied in the case of triangular external force. A more complex system consisting of a mixed viscous-hysteretic damping was studied in Bae et al. 56 using this method.
However, unlike these previous methods that require some special integration techniques based on zero or first order hold methods, the modified method proposed in this paper makes it possible to solve the system in state space form using standard numerical integration algorithms, such as the Runga-Kutta family of solvers that are readily available within MATLAB. A detailed discussion about the method, in particular the time reversal technique, is presented next.

Time reversal technique
In general, a continuous structural system can be observed by applying a discretisation process. Thus, the system can be simplified into a set of elements with a certain number of degrees of freedom (DOFs). The dynamic equation of motion for an MDOF system can be expressed in matrix form as The sizes of the matrices in Equation (19) are defined by the number of DOFs, n. The first and second derivatives with respect to the time are denoted by the dot,
The mass, stiffness and damping matrices are associated with the properties of the structure, that is, material weight, elasticity and damping coefficient. In particular, there are various ways to define the damping of the system, that is, viscous damping and hysteretic damping. The viscous damping force is represented in Equation (19) as the term C .
x. However, in a linear hysteretic damping case, the damping force is a displacement-dependent term. It is usually modelled proportional to the complex value of the stiffness. Hence, Equation (19) becomes where G is the matrix of damping coefficient that represents the loss factor of each element of the system. It is a common practice to avoid using complex coefficients for transient response or time-domain analysis, because in this case the damping model does not have a direct physical realisation. For this reason, the linear hysteretic damping is often represented by an equivalent viscous damping. Assuming constant amplitude oscillatory motion with a frequency , the two damping forces are identical if Hence, an equivalent viscous damping can be selected for a particular frequency of interest as The selection of is typically the dominant frequency at which the damping is active, that is, the first natural frequency. However, in reality, the individual element damping occurs at different frequencies. Thus, accurate solutions can only be obtained if the excitation frequency is at or near the selected .
In the present work, it is essential to evaluate hysteretic damping in its original complex form to obtain an optimum design of TIhD and TMhDI. Moreover, the investigation of earthquake excitation force is one of the interests in the current work. Therefore, to address this issue, the time reversal technique of Inaudi and Makris 41 is adopted and modified.
Inaudi and Makris suggested that the displacements and the forces can be evaluated in analytic forms, that is, complex signals constructed by real and imaginary components. Hence,
It is important to note that the vectors and matrices denoted by x a (2n × 1), A (2n × 2n), B (2n × 1) and L F (n × 1) are all in complex form. The force vector in analytic form, F a , is defined as the product of the complex time function, r a , with a coefficient vector of elements' load factor, L F , whereas 0 (n × 1) and I (n × 1) are the zeros and identity matrices, respectively.
Despite seeming that Equation (25) can be solved directly involving a standard integration procedure, that is, convolution integral, the eigenvalues of A are pairs of stable and unstable roots. 41 Therefore, a special treatment is required to obtain the stable solutions of Equation (25). Inaudi and Makris 41 introduced a backward integration technique that reverses the independent variable (time) in the modal coordinate of the unstable parts. The analytic vector x a can be rewritten by defining analytic modal coordinates, q k a and q l a , for the stable and unstable parts, respectively.
where q k a = q k + H[q k ] and q l a = q l + H[q l ], whereas (2n × 2n) is the eigenvector matrix of A. By introducing a reversed time, z = −t, and a functionq l (z) = q l (t), hence, where [ ′ ] defines the derivative in the reversed time. Therefore, Equation (25) can be rewritten as where A * = −1 A and B * = −1 B. To be noted here, A * is a diagonal matrix; thus, Equation (29) becomes uncoupled. Hence, components of q k a and q l a can be evaluated separately by means of forward time integration and backward time integration, respectively. Inaudi and Makris 41 numerically solved both a SDOF problem and a 2-DOF problem by applying discrete-time sampling to the analytical integral solutions of Equation (29). However, this approach may be difficult to implement in a more complex and higher DOFs system as it requires the analytical integral solution for each different system. Therefore, in the present work, a more robust computational approach is proposed. Herein, a variable separation procedure is applied to separate the real and imaginary components in Equation (29); hence, which can be evaluated separately for forward integration and backward integration In Equations (31) and (32), the matrices A * * (2n × 2n), Ã * * (2n × 2n), B * * (2n × 1) andB * * (2n × 1) consisted of components that provide coupling between the real and imaginary parts of each modal coordinate.
The key aspect in the present work is that Equations (31) and (32) are sets of ordinary differential equations (ODEs); thus, they can be straightforwardly solved via a standard numerical procedure, that is, by using the Runge-Kutta method. Therefore, it provides a more robust approach to evaluate any structural system as long as the mass, stiffness, damping matrices and the force vector are known. Further discussion on the computational algorithm of the present approach is given in Appendix B1.

Validation
In this subsection, validations of the new computational approach are presented. An example of a SDOF system excited by a harmonic sine-wave input is now considered. The results of the present computational approach are compared with those obtained via the zero-order hold method of Inaudi and Makris. 41 Figure 5A,B shows the time domain response of a simple SDOF system with one mass and one coupled spring-damping element represented by a complex stiffness subjected to sinewave ground displacement. It is obvious that the solution using the present method are in a good agreement with the zero-order hold method by Inaudi and Makris 41 around the resonance ( = n ) and away from resonance ( = 30 n ). Here, is the frequency of excitation and n is the natural frequency of the system. In addition, a comparison with a common method using equivalent viscous damping is also presented when the equivalency is made at the resonance. As can be seen in Figure 5B, the equivalent viscous damping method becomes inaccurate at frequencies away from the natural frequency.
The following section will discuss further case studies where the time reversal method is applied for analysis of structures with TIhD and TMhDI.

Harmonic base excitation
An SDOF structure as shown in Figure 1A for n = 1, mass m and stiffness k 0,1 set to 1kNs 2 ∕m and 0.0987kN∕m, respectively was used in this example. Figure 6A shows the performance comparison between the four inerter-based damper systems applied to an SDOF structure. Their optimised parameters are given in Table 3. It can be seen from this figure that the performance of the TID and the TMDI is similar. This can be understood because the effect of m d , which is only 5% of the inertance, is small across the frequencies. However, the TMDI has slightly higher response around the resonant frequency after the optimisation. Its time domain response is given in Figure 7A,B obtained by using the method presented in the previous section. A similar result also can be seen for the TIhD compared with the TMhDI, which indicates that the effect of m d is very small.
Although the response amplitude around the resonant frequency is slightly higher with the hysteretic damping, it has a potential benefit in the higher frequency region as previously discussed. The attenuation difference between TID/TMDI and TIhD/TMhDI in the high frequency region is 20 dB/dec. In civil engineering applications, especially for low frequency structures such as base-isolated structures, this is a valuable benefit because the response around the resonance frequency needs to be reduced to deal with long-period earthquakes, while at the same time maintaining lower response at the higher frequency region.
To explore MDOF structures, a 3-DOF structure as shown in Figure 1A for n = 3 was selected for a case study. The parameters of the structure were designed to be the same with the 3-DOF structure presented in Lazar et al., 14 where m 1 = m 2 = m 3 = 1kNs 2 ∕m and k 0,1 = k 1,2 = k 2,3 = 1500kN∕m. All of the inerter-damper devices shown in Figure 1B -D were located at the bottom storey as this is their optimum location. Figure 6B shows the transfer function X 3 ∕R of a 3-DOF structure in the frequency domain. Its time domain response is given in Figure 7C,D. The optimised parameters of the inerter-damper devices are given in Table 3 derived based on a fixed value of = 0.18. The structural response amplitude around the first resonance frequency is similar for all cases. As mentioned previously for SDOF structures, the hysteretic damping gives a slightly higher response after the optimisation. However, around the second and third resonance frequencies, the response of the structure with TIhD and TMhDI is significantly higher when compared with the TID and TMDI with no hysteretic damping. This is because the force transferred to the structure by the TIhD and the TMhDI is larger than that of TID and TMDI. This also indicates that modelling a system that has hysteretic damping with a viscous damping model may significantly overestimate the level of damping that can be achieved in higher modes of vibration. As a result, using the model proposed in this paper may have considerable benefits for modelling physical systems of this type.    From Figure 9, it can be observed that overall the optimisation by SADE algorithm works better than the FPT, particularly for low-frequency earthquakes (Earthquakes [11][12][13][14][15][16][17][18][19][20]. For Earthquakes 9, 29 and 30, the FPT is actually better than SADE. It is not clear what the reason for this is, but one possible explanation is that the objective function given in Equation (33) is designed to seek for the minimum overall average, not the minimum for each of the earthquakes. Overall, the average is −24.66% for the bars below FPT and 5.18% for the bars above the FPT. This means that the optimisation using the SADE algorithm gives a much better response compared with the FPT across all 30 earthquake cases. In terms of applications, computations of this type can be used during the design stage, in order to give estimates of the level of displacement responses in a structure where the damping in the inerter device is hysteretic. Figure 10 shows a comparison of the performance of the four devices for each earthquake, for the three different cases considered. Each of the device was optimised by using the SADE algorithm. The comparison was made relative to the TID performance, which was the best overall. So in this figure, each bar represents the difference between each device and TID (which is taken as 0%) for each earthquake.
An example of the time domain response is given in Figure 11 for both SADE and FPT optimisations. The trends observed in both Figures 9 and 10 can also be observed in these time domain plots. Specifically, the overall amplitudes  Figure 11B for SADE are significantly less that the amplitudes in Figure 11A for FPT. It is also possible to see that the hysteretic damping devices have higher amplitude responses than the TID or TMDI in both plots.
From Figures 10, we can infer that the effect of the secondary mass m d is small due to the fact that the TMDI performance is only up to 6% different to the TID performance. The effect of the hysteretic damping was much larger, because the TIhD and TMhDI were up to 33% different to the TID results. This implies that when the physical damping behaviour in an inerter-based device is hysteretic, then the viscous damping approximation could lead to results that overestimate the effectiveness of the device.
The methodology presented in this paper offers an alternative modelling approach that can capture the effect of the hysteretic damping in the inerter-based device for arbitrary input signals such as earthquakes. As a result, it should be possible to obtain more accurate estimates of the device performance during a design procedure.

CONCLUSIONS
This paper has introduced two novel passive inerter-based dampers, namely, the TIhD and TMhDI. Both devices were based on the well-established idealised concepts of the TID and TMDI with the viscous damping element being replaced by hysteretic damping element. The primary motivation was in order to create models that were able to capture the physical behaviour of structural systems in a more realistic way.
Among all the devices considered, the TMhDI is believed to be potentially the most realistic because it has both hysteretic damping and a small secondary mass, which captures the effect of the mass of the inerter device itself. It has been shown how the optimum parameters of both devices can be obtained via fixed-point theory although one limitation is that additional fine tuning is required for the devices with hysteretic damping. In addition, it was shown that there are regions of the frequency domain where the hysteretic damping leads to different transmissibility values compared to viscous damping. This was observed for both the SDOF and MDOF examples considered. For example, for the SDOF structure with harmonic excitation, both TIhD and TMhDI have the potential benefit of reducing displacement response in the higher frequency region, and the attenuation difference is of the order of 20 dB/dec. However, in the MDOF example, both hysteretic devices exhibit significantly higher amplitudes of displacement response for the higher modes of vibration compared with the viscous damping model. This is something to be aware of in systems where the physical behaviour is closer to hysteretic than viscous.
For the earthquake signals considered here, the SADE algorithm was used to optimise the devices' parameters, and it was shown that the effect of additional small m d , in this case is 5%, is not significant for the example considered. The TIhD and the TMhDI performance were shown to have a larger response than the TID and TMDI in the earthquake example. Moreover, the obtained optimum stiffness and damping parameters of the TIhD and TMhDI for the example considered do not satisfy the standard equivalent viscous damping relationship that would enable a conversion from the TIhD to TID or from TMhDI to TMDI. This implies that the TIhD and the TMhDI are better treated without using and equivalent viscous damping assumption when performing seismic design analysis. This study has been made possible by first developing a time domain method via time reversal technique for the analysis of structural systems with linear hysteretic damping. The TIhD and TMhDI models subsequently developed offer an approach that can capture the systems' performance in a more realistic way when the physical behaviour of the system is closer to hysteretic than viscous. The method has also been validated by comparing the results with the zero-order hold method of Inaudi and Makris and with equivalent viscous damping. Future work will include verification of the method using experimental test results. 56. Bae S, Jeong W, Cho J, Lee S. Transient response of vibration systems with viscous-hysteretic mixed damping using Hilbert transform and effective eigenvalues.

APPENDIX A: COMPUTATIONAL ALGORITHM
As previously discussed, the main feature here is to obtain the solutions of the real and imaginary parts of both the modal responses in forward and reversed times from Equations (31) and (32). Herein, therefore, the process of Hilbert transformation is extremely important to obtain the accurate solutions. Furthermore, as stated earlier, these state space equations need to satisfy the boundary conditions at t = −∞ and t = ∞ in which both the displacement responses and forces are zeros. Therefore, for a finite time duration, nt, spans from an initial time, t 0 , to an ending time, t e , a numerical approximation is required to satisfy the boundary conditions. In the present work, it is defined that the force boundary conditions are F a (−∞) ≈ F a (t 0 − nt) = 0 ( A 1 )

FIGURE A1
Computational algorithm and F a (∞) ≈ F a (t e + nt) = 0, where nt = t e − t 0 . This definition also similarly applies to the displacement boundary conditions. Hence, prior to the Hilbert transformation of the original force vector, F, "dummy" initial and ending conditions are added to the original force function. As an example, if t 0 = 0s and t e = 100 s, then the extended force function, F ext , should span from t = −100 s to t = 200 s. This can be expressed as where the time parameters inside the square brackets define the time span of each force vector. The "dummy" F initial and F ending will depend to the shape of the original force function, F. In the present case, if it is an impulse or random excitation force function, then F initial and F ending are assumed to be zeros. However, if the original is a periodic function, that is, harmonic force function, F initial and F ending need to be constructed as decaying functions such that their values at the time boundaries are zeros.
The next step is then to obtain the Hilbert transformation of the extended force function, F ext . Thus, the earlier analytical force function as expressed in Equation (24) becomes and then flipped to obtain the reversed-time force function, that is, if F a (t) spans from −100 to 200 s, thenF a (z) spans from 200 to −100s.
After the force function is constructed, the following steps are to construct the state space forms as given in Equations (31) and (32). The detailed algorithm of the present computational approach is depicted in Figure A1.
In the present work, a built-in Runge-Kutta module in MATLAB, so-called ODE45 module, is used to obtain the solutions of the sets of ordinary differential equations given in Equations (31) and (32). It is important to note that once the solutions of Equation (32) is obtained, it is still in the reversed time, z. Hence, the solutions needs to be reversed back to the original time, t, as depicted in Figure A1.
Once the modal responses vectors, q k a (t) and q l a (t), are computed, then the displacement-velocity vector, x a (t), can be calculated via Equation (23). The solutions are then isolated for the finite duration time, nt. The real parts of the solutions are the actual responses, x(t), .
x(t), of the system due to the force F(t).