Optimization of Multiple Tuned Mass Damper (MTMD) Parameters for a Primary System Reduced to a Single Degree of Freedom (SDOF) through the Modal Approach

Abstract

The research paper presents a novel approach toward constructing motion equations for structures with attached MTMDs (multiple tuned mass dampers). A primary system with MDOF (multiple dynamic degrees of freedom) was reduced to an equivalent system with a SDOF (single degree of freedom) through the modal approach, and equations from additional MTMDs were added to a thus-created system. Optimization based on ℌ₂ and ℌ_∞ for the transfer function associated with the generalized displacement of an SDOF system was applied. The research work utilized GA (genetic algorithms) and SA (simulated annealing method) optimization algorithms to determine the stiffness and damping parameters for individual TMDs. The effect of damping and stiffness (MTMD tuning) distribution depending on the number of TMDs was also analyzed. The paper also reviews the impact of primary system mass change on the efficiency of optimized MTMDs, as well as confirms the results of other authors involving greater MTMD effectiveness relative to a single TMD.

1. Introduction

The first general formulation of the research problem involving the parameter optimization in a tuned mass dampers (TMD) attached to a single degree of freedom (SDOF) structure was suggested by Den Hartog [1]. He determined known formulas for TMD tuning and a critical damping ratio in consideration of a TDOF (two degrees of freedom) system with harmonic excitation. Recommendations in terms of the optimal TMD parameter set can also be found in the work by Warburton [2] for a deterministic case, that is, when the main mass of a SDOF system undergoes harmonic excitation. The response of a linear system to broadband response is different in the case of a harmonic excitation, because the first one occurs at system’s natural frequencies. In cases where broadband excitation has an almost constant spectrum within the natural frequency range, it is convenient to replace it with white noise. In his work [3], Warburton provided formulas for optimal TMD parameters at exactly such excitation. Similar discussions were conducted by Bakre and Jangid in [4,5]. A different approach to the issue of determining TMD parameter values was proposed by Krenk [6], who suggested aligning the ordinates of three points A, B (just like Den Hartog), and one central, between the previous two. He derived a new optimal damping factor value, which is 15% higher than the classic result.

The most recent work that needs mentioning includes the one by Batou and Adhikari [7], which focused on an SDOF system with an attached viscoelastic damper. A standard rheological model as a TMD model was used to obtain analytically optimal TMD parameters. Classic results for a TMD with viscous damping can be obtained as a specific case for this damper type. It was also demonstrated that using a TMD with viscoelastic damping enables obtaining better vibration absorption compared to an equivalent TMD with viscous damping. Another new approach to a well-known TMD attached to an SDOF system was reviewed in the work by Marian and Giaralis [8]. The authors proposed a tuned mass damper inerter (TMDI), which combines a classic tuned mass damper (TMD) with a grounded inerter, which counteracts relative acceleration through additional inertia. The work utilized a modified “two point” approach proposed by Den Hartog [1], which enabled determining optimal TMDI parameters. TMD analyses are conducted in various fields of science. The study by Peterka et al. [9] involved modeling the vibrations of a drilling mandrel with a vibration damper linked with the mandrel via a viscoelastic coupling. A system with three degrees of freedom with an additionally installed TMD was modeled and the efficiency of such a solution demonstrated.

Multiples TMDs (MTMDs) can be introduced to a system in order to improve structural parameters. The efficiency of a spatially distributed MTMD was originally discussed by Bergman et al. [10]. The presented results indicated that structural behavior could be improved when an MTMD is tuned for appropriate structural frequencies, maintaining the total mass ratio at the same level as in the case of a TMD. It was demonstrated that an increased number of dampers resulted in the flattening of the response curve over a wider frequency range (Kareem and Kline [11]). The performance of TMD and MTMD was compared, and it was shown that the optimal damping values of an MTMD damper were lower than in the case of TMD. The early stages of studies concerning MTMD configurations under simplified and constrained conditions primarily involved to reducing the number of associated design variables. For example, in [12], Xu and Igusa studied MTMD with evenly spaced natural frequencies and an equal constant damping, and based on the asymptotic analysis, they demonstrated that such an MTMD was efficient in reducing primary system response. For a finite number of MTMDs with similar constraints, Joshi and Jangid [13] and Jangid [14] presented optimal MTMD parameters for an undamped and damped primary system, respectively. MTMDs with equal damping factors and evenly spaced natural frequencies were also studied by such authors as Yamaguchi and Harnpornchai [15].

A number of algorithms for selecting optimal MTMD parameter values have been proposed over the years (see Zuo and Nayfeh [16]). The study focused on developing optimization algorithms based on the ℌ₂ and ℌ_∞ norms and aimed at minimizing the response to random and harmonic excitations. ℌ_∞ optimization was also utilized in the work by Zhao [17], for determining optimal ATID (active tuned inerter damper) parameters. The article by Aggumus and Guclu [18] describes a study involving a single tuned mass damper (STMD) using a model of a building with multiple degrees of freedom (MDOF). It applied an optimization method based on the H_∞ norm in the STMD control process. Asami, in the work [19], presents an exact ℌ_∞ optimal solution for a serial TMD with a double mass connected to a damped primary system. In the optimization of ℌ_∞ of the transfer function, a very precise numerical solution was successfully obtained by solving the sixth-order algebraic equation.

MTMD parameters were studied by Zuo and Nayfeh [20], who utilized ℌ₂ optimization. They demonstrated that optimal structures did not have uniformly distributed tuning frequencies or identical damping factors for individual TMDs, and the applied optimization of individual parameters within an MTMD system provided significant improvement of the entire SDOF system with MTMD. In response to the aforementioned study, Li and Ni [21] presented an optimization method for MTMDs with non-uniformly distributed mass, based on a gradient with linear search of the objective function. By solving the problem of multi-target optimization based on the objective function in the mode of maximum displacement or frequency response, the authors obtained optimized, unevenly distributed MTMDs. In work [22], Zuo and Nayfeh applied a method utilizing subgradients of a non-smooth objective function. Besides the gradient methods, also genetic algorithms are used in terms of MTMD optimization, for example, by Ok, Song, and Park [23]. Heuristic methods, such as the particle swarm optimization (PSO), were also used. The PSO method for the determination of optimal MTMD parameters was used, among others, in the research by Leung and Zhang [24] and Zhang et al. [25]. Other methods used for optimizing the objective function in MTMD issues include simulated annealing (SA) applied by Aydin et al. [26] for optimizing the parameters in a multi-story building.

Recent works on MTMD include the one by Kim and Lee [27]. The article discusses the analysis of linear MTMDs with numerous practical configurations, attached to an SDOF structure, subjected to white noise input. Six practical configurations were developed and analyzed comparatively, each of which was linearly restricted with distributed tuning factors, mass ratios, damping factors, and their combinations. In [28], Stanikzai et al. analyzed TMDs and MTMDs distributed at various degrees of freedom of an MDOF system in a building with base isolation (BI). The analysis involved 40 earthquake ground motions for the adopted pattern of an MDOF structure with attached TMD and MTMDs and indicated the efficiency of the latter. In [29], Yin et al. analyzed a new type of TMD system named pounding tuned mass damper (PTMD). The coupled equations were created by combining the equations of motion of both the bridge and moving vehicles. In order to compare the damping performance, a parametric study of the various numbers and locations, mass ratio, and stiffness of the MPTMDs were investigated.

The article presents an alternative method for developing motion equations for a MTMD structure. The method enables adding single TMDs or MTMD groups to completely different degrees of freedom of the primary system. The system of equations allows for easy MTMD tuning for complex vibrations modes, with MTMDs located in local maxima of these vibrations modes, while still analyzing the SDOF system with attached MTMDs. MTMD parameter optimization based on ℌ₂ and ℌ_∞ norms for the transfer function associated with the modal coordinate of an SDOF system was proposed. The research work utilized GA (genetic algorithms) and SA (simulated annealing method) optimization algorithms to determine the stiffness and damping parameters for individual TMDs. The comparison and discussion were based on previously developed models of a primary system with an attached MTMD.

2. Methods

2.1. Original Proposal of Motion Equations for a Primary System with Attached MTMDs, Reduced to an SDOF System through the Modal Approach

In order to be able to model an equivalent system with a single degree of freedom (SDOF) for a system with multiple degrees of freedom using the modal approach, several conditions need to be satisfied, such as the spatial excitation force distribution must be relatively uniform and, in the case of concentrated load, it should be located near the highest ordinate of the first vibration mode and natural frequencies higher than ${\omega }_i$ are not close to this frequency.

Using this relationship enables obtaining an equation, which describes an equivalent motion of a SDOF system, corresponding to a system with N degrees of freedom. The primary system is reduced through the modal approach. The main assumption is the orthogonality of the main structure C damping matrix and a main structure is treated as a linear time-invariant system (LTI system). It is a system that produces an output from any input signal subject to linearity and time constraints. In general, for N degrees of primary system, the motion equation solution:

$$M\ddot{q}(t)+C\dot{q}(t)+Kq(t)=p(t)$$

(1)

expands into a series of eigenvectors

$$q(t)={\displaystyle \sum _{i=1}^N{a}_i{\psi }_i(t)}=W\psi (t),$$

(2)

where $\psi (t)$ is the modal coordinate vector and $a_i$ is i-th eigenvector, while W is the matrix corresponding mode shapes in its columns.

If we include the $N_i$ mode of the vibrations, which have the greatest share in the vibration system, the above equation can be expressed as:

$$q(t)={\displaystyle \sum _{i=1}^{N_i}{a}_i{\psi }_i(t)}.$$

(3)

If we normalize the eigenvectors in the manner shown below, we get the following motion equations:

$${\tilde{m}}_i{\ddot{\psi }}_i(t)+{\tilde{c}}_i{\dot{\psi }}_i(t)+{\tilde{k}}_i{\psi }_i(t)={\tilde{p}}_i(t),$$

(4)

where $i=1,2,\dots N_i$, ${\tilde{m}}_i=a_i^{T}M{a}_i=1$, ${\tilde{c}}_i=a_i^{T}C{a}_i=2{\zeta }_i{\omega }_i$, ${\tilde{k}}_i=a_i^{T}K{a}_i={\omega }_i^2$, and ${\tilde{p}}_i(t)=a_i^{T}p(t)$.

Designations as in Figure 1 were adopted for the analyzed situation. The diagram shows a beam structure with eliminated degrees of rotation and degrees of freedom along the X-axis direction of the global coordinate system.

The following designations for the diagram above were adopted: j = 1,2…N—system degrees of freedom; I = 1,2…N_i– number of vibration modes taken into account (further vibration eigenmodes not necessarily included in the analysis)—i index means any subsequent mode number included; and q_j(t)—displacement of j–th degree of freedom.

Assuming that natural frequencies corresponding to the modes, which are decisive for structural vibrations, are separated from each other, the dynamic response in the case of excitations with a selected vibration mode can be determined using the formula:

$$q_j(t)={\displaystyle \sum _{i=1}^{N_i}{a}_{ji}{\psi }_i(t)},$$

(5)

where ${\psi }_i(t)$ is the modal coordinate associated with the i-th vibration eigenmode, $a_{ji}$ is the coordinate of the i-th vibration mode of the j–th degree of freedom, and N_i is the number of vibration eigenmodes taken into account.

The S point to which the system is reduced is called the equivalent system mass reduction point. It is adopted at the point of the highest coordinate of the $a_i$ mode shape or the point of application of concentrated load within the system, including, e.g., force originating from a mechanical vibration damper.

After bilaterally dividing the modal motion Equation (4) by the coordinate of the i-th mode shape of the equivalent system mass reduction point $a_{ji}$ and after normalizing eigenvectors in the form of ${\tilde{m}}_i=a_i^{T}M{a}_i=1$, we get an equivalent system motion equation:

$$\frac{{\tilde{m}}_i}{a_{ji}{}^2}{a}_{ji}{\ddot{\psi }}_i(t)+\frac{{\tilde{c}}_i}{a_{ji}{}^2}{a}_{ji}{\dot{\psi }}_i(t)+\frac{{\tilde{k}}_i}{a_{ji}{}^2}{a}_{ji}{\psi }_i(t)=\frac{1}{a_{ji}}{\tilde{p}}_i(t),$$

(6)

or

$${\overline{m}}_{ji}{\ddot{\overline{q}}}_{ji}(t)+{\overline{c}}_{ji}{\dot{\overline{q}}}_{ji}(t)+{\overline{k}}_{ji}{\overline{q}}_{ji}(t)={\overline{p}}_{ji}(t),$$

(7)

where $i=1,2,..,N_i$, ${\overline{q}}_{ji}(t)=a_{ji}{\psi }_i(t)$ is the equivalent system displacement, ${\overline{m}}_{ji}={\tilde{m}}_i/a_{ji}{}^2=1/a_{ji}{}^2$ is the equivalent mass, ${\overline{c}}_{ji}={\tilde{c}}_i/a_{ji}{}^2=2{\zeta }_i{\omega }_i/a_{ji}{}^2$ is the equivalent damping of system, ${\overline{k}}_{ji}={\tilde{k}}_i/a_{ji}{}^2={\omega }_i{}^2/a_{ji}{}^2$ is the equivalent stiffness of system, and ${\overline{p}}_{ji}(t)={\tilde{p}}_i(t)/a_{ji}$ is the equivalent force of system.

A new tuned mass damper (TMD) or multiple TMDs (MTMDs) can be attached to the new system. The well-known Kelvin–Voigt model was adopted for the description of each TMD. In general, attached TMDs do not have to be located at the S system reduction point. If the S system reduction point is an attachment point for a single TMD, we get a 2DOF system, which is a well-known issue addressed by Den Hartog in [1], in the case of no damping of the primary system. Of course, TMD is tuned to a frequency near ${\omega }_i$, for which the equivalent system was determined. The case of using MTMDs located in various degrees of freedom of the primary system requires a separate discussion. MTMD arrangement diagram is shown in Figure 2. Of course, each TMD shall be tuned to near frequency ${\omega }_i$.

The following designations for the above diagram were adopted: k = 1,2…N_k is the TMD degrees of freedom; $m_k^t,c_k^t,k_k^t$ are the mass, damping, and stiffness of the k-th damper; and ${\xi }_k(t),q_k^t(t)$ is the relative and absolute displacement of the k-th TMD. If we add an $N_k$ number of tuned mass dampers to the created system, the new system will have $N_k+1$ degrees of freedom (each TMD is an additional degree of freedom).

If we introduce location vectors for each MTMD, with the value 1 present on the degree of freedom to which the tuned mass damper is attached, we can write the following structural displacement equation for a degree of freedom to which the k-th TMD is attached:

$${\tilde{q}}_k^t(t)=e_k{}^T{a}_i{\psi }_i(t)={\tilde{a}}_{ki}{\psi }_i(t),$$

(8)

where $e_k{}^T=\left[0,0,1,\dots ,0\right]$ is the location vector of the k–th damper. The value 1 is present on the degree of freedom to which the k-th TMD is attached and ${\tilde{a}}_{ki}=e_k{}^T{a}_i$ is the ordinate of the i-th eigenvector at the TMD attachment point.

We also have to introduce the value of the displacements, velocity, and relative acceleration, which describe the motion of additional MTMDs. Generalized displacements, velocities, and relative accelerations for the dampers can be expressed with the following formulas:

$${\xi }_k^{}=q_k^t-{\tilde{q}}_k^t=q_k^t-{\tilde{a}}_{ki}{\psi }_i=q_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\overline{q}}_{ji},$$

(9)

$${\dot{\xi }}_k^{}={\dot{q}}_k^t-{\dot{\tilde{q}}}_k^t={\dot{q}}_k^t-{\tilde{a}}_{ki}{\dot{\psi }}_i={\dot{q}}_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\dot{\overline{q}}}_{ji},$$

(10)

$${\ddot{\xi }}_k^{}={\ddot{q}}_k^t-{\ddot{\tilde{q}}}_k^t={\ddot{q}}_k^t-{\tilde{a}}_{ki}{\ddot{\psi }}_i={\ddot{q}}_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\ddot{\overline{q}}}_{ji}.$$

(11)

A system of motion equations for such a case, with the number of degrees of freedom of $N_k+1$, is shown by the formula below:

$$\begin{array}{l}{\overline{m}}_{ji}{\ddot{\overline{q}}}_{ji}(t)+{\overline{c}}_{ji}{\dot{\overline{q}}}_{ji}(t)+{\overline{k}}_{ji}{\overline{q}}_{ji}(t)={\overline{p}}_{ji}(t)+T_{ji}\\ m_k^t{\ddot{q}}_k^t+c_k^t\left({\dot{q}}_k^t-{\dot{\tilde{q}}}_k^t\right)+k_k^t\left(q_k^t-{\tilde{q}}_k^t\right)=0\end{array}$$

(12)

The force from each damper, which needs to be attached to the degree of freedom associated with the k–th TMD, taking into account absolute values, can be expressed with the following formulas:

$$T_k=-m_k^t{\ddot{q}}_k^t,$$

(13)

$$T_k=c_k^t\left({\dot{q}}_k^t-{\dot{\tilde{q}}}_k^t\right)+k_k^t\left(q_k^t-{\tilde{q}}_k^t\right)=c_k^t\left({\dot{q}}_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\dot{\overline{q}}}_{ji}\right)+k_k^t\left(q_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\overline{q}}_{ji}\right).$$

(14)

Changing to generalized force, for absolute values we get:

$$T_{ji}=\frac{1}{a_{ji}}{\displaystyle \sum _{k=1}^{N_k}{\tilde{a}}_{ki}{T}_k}=-{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}}{a_{ji}}{m}_k^t{\ddot{q}}_k^t},$$

(15)

$$T_{ji}=\frac{1}{a_{ji}}{\displaystyle \sum _{k=1}^{N_k}{\tilde{a}}_{ki}{T}_k}={\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}}{a_{ji}}{c}_k^t}\left({\dot{q}}_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\dot{\overline{q}}}_{ji}\right)+{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}}{a_{ji}}{k}_k^t}\left(q_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\overline{q}}_{ji}\right).$$

(16)

By substituting the additional damping force values to motion equations, and by using the relationships for absolute values in Equations (15) and (16), we get:

$$\begin{array}{l}{\overline{m}}_{ji}{\ddot{\overline{q}}}_{ji}(t)+{\overline{c}}_{ji}{\dot{\overline{q}}}_{ji}(t)+{\overline{k}}_{ji}{\overline{q}}_{ji}(t)+{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}{}^2}{a_{ji}{}^2}{c}_k^t}{\dot{\overline{q}}}_{ji}-{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}}{a_{ji}}{c}_k^t}{\dot{q}}_k^t+{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}{}^2}{a_{ji}{}^2}{k}_k^t}{\overline{q}}_{ji}-{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}}{a_{ji}}{k}_k^t}{q}_k^t={\overline{p}}_{ji}\\ m_k^t{\ddot{q}}_k^t+c_k^t\left({\dot{q}}_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\dot{\overline{q}}}_{ji}\right)+k_k^t\left(q_k^t-\frac{{\tilde{a}}_{ki}}{a_{ji}}{\overline{q}}_{ji}\right)=0\end{array}$$

(17)

Motion equations in matrix form, for the most general case, with absolute values, are shown by the relationship below:

$$\overline{M}\ddot{\overline{q}}+\overline{C}\dot{\overline{q}}+\overline{K}\overline{q}=\overline{p},$$

(18)

where M, C, and K are positive-defined mass, damping, and stiffness matrices, respectively (it is impossible to assume negative values of stiffness and damping of individual TMDs).

From the generalized displacement and excitation forces’ vector let us separate vector blocks associated with the equivalent displacement modal ordinate and actual TMD displacements:

$$\ddot{\overline{q}}=\left[\begin{array}{c}{\ddot{\overline{q}}}_1\\ {\ddot{\overline{q}}}_2\end{array}\right],\dot{\overline{q}}=\left[\begin{array}{c}{\dot{\overline{q}}}_1\\ {\dot{\overline{q}}}_2\end{array}\right],\overline{q}=\left[\begin{array}{c}{\overline{q}}_1\\ {\overline{q}}_2\end{array}\right],\overline{p}=\left[\begin{array}{c}{\overline{p}}_1\\ {\overline{p}}_2\end{array}\right],$$

(19)

where:

$${\ddot{\overline{q}}}_1=\left[{\ddot{\overline{q}}}_{ji}\right],{\dot{\overline{q}}}_1=\left[{\dot{\overline{q}}}_{ji}\right],\overline{q}=\left[{\overline{q}}_{ji}\right],\overline{p}=\left[{\overline{p}}_{ji}\right],$$

(20)

$${\ddot{\overline{q}}}_2=\left[\begin{array}{c}{\ddot{q}}_1^t\\ {\ddot{q}}_2^t\\ ⋮\\ {\ddot{q}}_k^t\\ ⋮\\ {\ddot{q}}_{N_k}^t\end{array}\right],{\dot{\overline{q}}}_2=\left[\begin{array}{c}{\dot{q}}_1^t\\ {\dot{q}}_2^t\\ ⋮\\ {\dot{q}}_k^t\\ ⋮\\ {\dot{q}}_{N_k}^t\end{array}\right],{\overline{q}}_2=\left[\begin{array}{c}{q}_1^t\\ q_2^t\\ ⋮\\ q_k^t\\ ⋮\\ q_{N_k}^t\end{array}\right],{\overline{p}}_2=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ ⋮\\ 0\end{array}\right].$$

(21)

Similarly, the blocks associated with the discussed degrees of freedom shall also be separated from the matrix $\overline{M}{,}_{}\overline{C}{,}_{}\overline{K}$:

$$\overline{M}=\left[\begin{array}{cc}{\overline{M}}_{11}& {\overline{M}}_{12}\\ {\overline{M}}_{21}& {\overline{M}}_{22}\end{array}\right],\overline{C}=\left[\begin{array}{cc}{\overline{C}}_{11}& {\overline{C}}_{12}\\ {\overline{C}}_{21}& {\overline{C}}_{22}\end{array}\right],\overline{K}=\left[\begin{array}{cc}{\overline{K}}_{11}& {\overline{K}}_{12}\\ {\overline{K}}_{21}& {\overline{K}}_{22}\end{array}\right].$$

(22)

The size of blocks ${\overline{M}}_{11}$, ${\overline{C}}_{11}$, ${\overline{K}}_{11}$ is $1\times 1$; ${\overline{M}}_{12}$, ${\overline{C}}_{12}$, ${\overline{K}}_{12}$ is $1\times N_k$; ${\overline{M}}_{21}$, ${\overline{C}}_{21}$, ${\overline{K}}_{21}$ is $N_k\times 1$; and of ${\overline{M}}_{22}$, ${\overline{C}}_{22}$, ${\overline{K}}_{22}$ is $N_k\times N_k$. Below you can find the forms of individual blocks of matrices $\overline{M}{,}_{}\overline{C}{,}_{}\overline{K}$.

$${\overline{M}}_{11}=\left[{\overline{m}}_{ji}\right],{\overline{M}}_{22}=\left[\begin{array}{cccccc}{m}_1^t& 0& \cdots & 0& \cdots & 0\\ 0& m_2^t& \cdots & 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & m_k^t& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& \cdots & m_{N_k}^t\end{array}\right],$$

(23)

$${\overline{M}}_{12}={\overline{M}}_{21}^T=\left[\begin{array}{cccccc}0& 0& \cdots & 0& \cdots & 0\end{array}\right],$$

(24)

$${\overline{C}}_{11}=\left[{\overline{c}}_{ji}+{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}{}^2}{a_{ji}{}^2}}{c}_k^t\right],{\overline{C}}_{22}=\left[\begin{array}{cccccc}{c}_1^t& 0& \cdots & 0& \cdots & 0\\ 0& c_2^t& \cdots & 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & c_k^t& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& \cdots & c_{N_k}^t\end{array}\right],$$

(25)

$${\overline{C}}_{12}={\overline{C}}^T{}_{21}=\left[\begin{array}{cccccc}-\frac{{\tilde{a}}_{1i}}{a_{ji}}{c}_1^t& -\frac{{\tilde{a}}_{2i}}{a_{ji}}{c}_2^t& \cdots & -\frac{{\tilde{a}}_{ki}}{a_{ji}}{c}_k^t& \cdots & -\frac{{\tilde{a}}_{N_{k}i}}{a_{ji}}{c}_{N_k}^t\end{array}\right],$$

(26)

$${\overline{K}}_{11}=\left[{\overline{k}}_{ji}+{\displaystyle \sum _{k=1}^{N_k}\frac{{\tilde{a}}_{ki}{}^2}{a_{ji}{}^2}}{k}_k^t\right],{\overline{K}}_{22}=\left[\begin{array}{cccccc}{k}_1^t& 0& \cdots & 0& \cdots & 0\\ 0& k_2^t& \cdots & 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & k_k^t& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& \cdots & k_{N_k}^t\end{array}\right],$$

(27)

$${\overline{K}}_{21}={\overline{K}}^T{}_{12}=\left[\begin{array}{cccccc}-\frac{{\tilde{a}}_{1i}}{a_{ji}}{k}_1^t& -\frac{{\tilde{a}}_{2i}}{a_{ji}}{k}_2^t& \cdots & -\frac{{\tilde{a}}_{ki}}{a_{ji}}{k}_k^t& \cdots & -\frac{{\tilde{a}}_{N_{k}i}}{a_{ji}}{k}_{N_k}^t\end{array}\right].$$

(28)

The new system suggested by the authors significantly limits the time consumption of the calculations and also enables the optimization of MTMD parameters discussed below. In general, $N_i$ equivalent systems can be created (as many as natural frequencies are taken into account). Each of these systems, with the number i, is a system with $N_k+1$ degrees of freedom. Achieving full displacement at the S reduction point requires determining the response for each of the created systems and then applying the formula:

$${\overline{q}}_j(t)={\displaystyle \sum _{i=1}^{N_i}{\overline{q}}_{ji}(t)}.$$

(29)

This approach is obviously not convenient. However, the method of reducing the system to an SDOF is very useful for structures with simple static diagrams (beams, cantilever structures such as chimneys, and masts), when the lowest vibration natural frequency is taken into account. Next, let us determine the parameters of an SDOF equivalent system. For such a system, using the aforementioned equations, we can create new motion equations with attached single TMDs or MTMDs tuned to a distinguished natural frequency.

2.2. Determination of a Transfer Matrix Used to Determine the Objective Function in Optimization Issues

For a model described with the general motion equation, which is a system with N degrees of freedom, the transfer matrix $H(s)$, the single element of which $H_{ik}(s)$ is called a transfer function, is a Laplace transform of the impulse response $h_{ik}(t)$.

$$H_{ik}(s)={\displaystyle {\int }_{-\infty }^{\infty }{h}_{ik}(t)e^{-st}dt},$$

(30)

where impulse response $h_{ik}(t)$ is the i-th response of the system to the k–th excitation in the form of a unit impulse function ${\delta }_k$ (Dirac delta) applied at the initial moment $t=0$, j means an imaginary unit, and s is complex variable $s=\sigma +j\omega$. Therefore, in order to determine impulse responses in the i points of the system, the excitation shall be adopted in the following form:

$$p(t)={\left[0,0,\dots ,0,{\delta }_k(t),0,\dots ,0\right]}^T.$$

(31)

In the specific case in which $\sigma =0$, the input is a complex integral $e^{-j\omega t}$ at frequency $\omega$ and $H_{ik}(j\omega )$, viewed as function of $\omega$, is known as the frequency response function (FRF) of the system and is given by the Fourier transform. The below Equations (32) and (33) show a pair of Fourier transforms, which can be found in the work by Bendat [30]:

$$H_{ik}(j\omega )={\displaystyle {\int }_{-\infty }^{\infty }{h}_{ik}(t)e^{-j\omega t}dt},$$

(32)

$$h_{ik}(t)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{H}_{ik}(j\omega )e^{j\omega t}d\omega }.$$

(33)

By substituting the excitation force vector (31) to the Equation (1) and after applying the Fourier transform we get:

$$(K-{\omega }^{2}M+j\omega C)H(j\omega )=I.$$

(34)

In general, $H_{ik}(j\omega )$ is the complex quantity and can be presented as:

$$H_{ik}(j\omega )=\left|H_{ik}\right|e^{j{\Theta }_{ik}}=\mathrm{Re}{H}_{ik}+j\mathrm{Im}{H}_{ik},$$

(35)

$$\left|H_{ik}(j\omega )\right|=\sqrt{{(\mathrm{Re}{H}_{ik})}^2+{(\mathrm{Im}{H}_{ik})}^2},\arg H_{ik}={\Theta }_{ik}=\mathrm{arctg}\frac{\mathrm{Re}{H}_{ik}}{\mathrm{Im}{H}_{ik}},$$

(36)

where $\left|H_{ik}(j\omega )\right|,{\Theta }_{ik}$ are the frequency response function module and argument, respectively.

If we assume the denotation of the dynamic stiffness matrix in the form:

$$G(j\omega )=K-{\omega }^{2}M+j\omega C=\mathrm{Re}G+j\mathrm{Im}G,$$

(37)

we get a complex equation:

$$\mathrm{Re}G\mathrm{Re}H-\mathrm{Im}G\mathrm{Im}H+j\left[\mathrm{Im}G\mathrm{Re}H+\mathrm{Re}G\mathrm{Im}H\right]=I,$$

(38)

and after breaking down, two real equations:

$$\mathrm{Re}G\mathrm{Re}H-\mathrm{Im}G\mathrm{Im}H=I,$$

(39)

$$\mathrm{Im}G\mathrm{Re}H+\mathrm{Re}G\mathrm{Im}H=0.$$

(40)

Ultimately, we get a matrix of the real and imaginary section of the full frequency response matrix $H(j\omega )$:

$$\mathrm{Re}H={\left\{\mathrm{Re}G+\mathrm{Im}G{\left[\mathrm{Re}G\right]}^{-1}\mathrm{Im}G\right\}}^{-1}I,$$

(41)

$$\mathrm{Im}H={\left[\mathrm{Re}G\right]}^{-1}\mathrm{Im}G\mathrm{Re}H.$$

(42)

Of course, with an appropriate computational procedure, we can directly determine a complex frequency response matrix from the formula:

$$H(j\omega )={\left(G(j\omega )\right)}^{-1}.$$

(43)

Knowing the frequency response matrix, it is possible to determine the structural response within the frequency domain in the form of:

$$q(j\omega )=H(j\omega )p(j\omega ),$$

(44)

where $q(j\omega ),p(j\omega )$ are the Fournier transforms of the displacement and vector of excitation forces.

The structure’s load Fourier transform is expressed by the formula:

$$p(j\omega )={\displaystyle {\int }_{-\infty }^{\infty }p(t)e^{-j\omega t}}dt.$$

(45)

If we determine the Fourier retransform of the structure’s excitation forces vector, we get the following relationship:

$$q(t)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }H(j\omega )\left({\displaystyle {\int }_{-\infty }^{\infty }p(t)e^{-j\omega t}}dt\right)e^{j\omega t}}d\omega .$$

(46)

The equation for $q(t)$ determination based on the Duhamel integral is often used in the literature. The matrix form of the equation is expressed by the formula:

$$q(t)={\displaystyle {\int }_0^{t}h(t-\tau )}p(\tau )d\tau ,$$

(47)

where $q(t)$ is a response vector, $h(t)$ is a matrix of impulse response, and $p(t)$ is the vector of excitation forces. Since it was assumed that the random excitation is a stationary process, the dynamic response is also a stationary process. Due to this fact, the correlation matrix of dynamic response of structure can be written as:

$$R_q(t_1,t_2)=\mathbb{E}\left[q(t_1)q^{\mathrm{T}}(t_2)\right].$$

(48)

Inserting the Equation (47) into the Equation (48), we assume:

$$R_q(t_1,t_2)={\displaystyle {\int }_0^{t_1}{\displaystyle {\int }_0^{t_2}h(t_1-{\tau }_1)\mathbb{E}\left[p(t_1)p^{\mathrm{T}}(t_2)\right]h(t_2-{\tau }_2)d{\tau }_{1}d{\tau }_2}}.$$

(49)

If we assume that the correlation of excitation forces is:

$$\mathbb{E}\left[p(t)p^{\mathrm{T}}(t-\tau )\right]=R_p(\tau )=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{S}_p(j\omega )e^{j\omega \tau }}d\omega ,$$

(50)

we can write the formula (49) in the form:

$$R_q(\tau )=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }H(j\omega )S_p(j\omega )H^{\mathrm{H}}(j\omega )e^{j\omega \tau }}d\omega ,$$

(51)

where $H^{\mathrm{H}}(j\omega )$ denotes the complex conjugate and transpose of frequency response matrix $H(j\omega )$ (Hermitian matrix) and $S_p(j\omega )$ is a full matrix of spectral density of random excitation. If we determine the values of correlation matrix at the moment $\tau =0$, the dependency for a matrix of variance of a dynamic response of a system is obtained:

$${\sigma }_q^2=R_q(0)=\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }H(j\omega )S_p(j\omega )H^{\mathrm{H}}(j\omega )}d\omega =\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{S}_q(j\omega )}d\omega ,$$

(52)

where $S_q(j\omega )$ is a matrix of double-sided spectral density of responses $S_q(j\omega )=H(j\omega )S_p(j\omega )H^{\mathrm{H}}(j\omega )$.

In the case of acting on the system of excitation of white noise nature $S_p(\omega )=I$, the matrix of spectral densities of response is the product of frequency response matrix $H(j\omega )$ and a matrix conjugated and transposed to $H^{\mathrm{H}}(j\omega )$. The standard deviation of response can be expressed by the formula:

$${\sigma }_q^{}={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }H(j\omega )H^{\mathrm{H}}(j\omega )}d\omega \right)}^{1/2}.$$

(53)

Now consider the p-norm of matrix $A=\left[A_{ij}\right]$. Let p ≥ 1 be a real number. The p-norm (also called $ℒ$_p-norm) is expressed by the formula:

$${\Vert A\Vert }_{{ℒ}_p}={\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{\left|A_{ij}\right|}^p}}\right)}^{1/p}.$$

(54)

For p = 2, we get the $ℒ$₂-norm:

$${\Vert A\Vert }_{{ℒ}_2}={\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{\left|A_{ij}\right|}_{}^2}}\right)}^{1/2}.$$

(55)

For a multivariable (multiple inputs and multiple outputs (MIMO)), LTI system with transfer matrix $H(s)=\left[H_{ij}(s)\right]$, the definition (55) generalizes to ℌ₂-norm:

$$\begin{array}{l}{\Vert H\Vert }_{{\mathscr{H}}_2}={\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{\Vert H_{ij}\Vert }_{{\mathscr{H}}_2}^2}}\right)}^{1/2}={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{\left|H_{ij}(j\omega )\right|}^2}}}d\omega \right)}^{1/2}\\ ={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{H}_{ij}(j\omega )H_{ij}^{}(-j\omega )}}}d\omega \right)}^{1/2}={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }\mathrm{tr}[H(j\omega )H^{\mathrm{H}}(j\omega )]}d\omega \right)}^{1/2}\end{array},$$

(56)

The notation ℌ₂ (instead of $ℒ$₂) is due to the fact that the function spaces that, in addition to having finite $ℒ$_p norms on the imaginary axis, are bounded and analytic functions in the right-half plane (RHP) (i.e., with no poles in the RHP) are called Hardy spaces ℌ_p. Thus, stable transfer functions belong to these spaces, provided the associated integral is finite. So we can see that the ℌ₂-norm of the transfer matrix $H(s)$ is the sum of diagonal values of the standard deviations ${\sigma }_q^{}$ matrix for white noise excitation obtained from the formula (53). Most of the authors ([16,20]) often use the expression ${\Vert H\Vert }_{{\mathscr{H}}_2}^2$, which is actually the square of the ${\Vert H\Vert }_{{\mathscr{H}}_2}$ norm, so it is not the ℌ₂-norm by definition.

For a stable single input and single output (SISO) linear system with transfer function $H(s)$, the ℌ₂-norm is defined as:

$${\Vert H\Vert }_{{\mathscr{H}}_2}={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left|H(j\omega )\right|}^2}d\omega \right)}^{1/2}.$$

(57)

In addition to the ℌ₂-norm, which gives a characterization of the average gain of a system, a more fundamental norm for systems is the ℌ_∞-norm, which provides a measure of a worst-case system gain. Consider a stable SISO linear system with transfer function H(s). The ℌ_∞-norm is defined as:

$${\Vert H\Vert }_{{\mathscr{H}}_{\infty }}=\underset{\omega \in ℝ}{\max }\left|H(j\omega )\right|,$$

(58)

or, in the event that the maximum may not exist, more correctly as:

$${\Vert H\Vert }_{{\mathscr{H}}_{\infty }}=\underset{\omega \in ℝ}{\sup }\left|H(j\omega )\right|.$$

(59)

Recall that $\left|H(j\omega )\right|$ is the factor by which the amplitude of a sinusoidal input with angular frequency $\omega$ is magnified by the system.

For multivariable MIMO systems, the ℌ_∞-norm is defined in an analogous way:

$${\Vert H\Vert }_{{\mathscr{H}}_{\infty }}=\underset{\omega \in ℝ}{\sup }\overline{\sigma }\left(H(j\omega )\right),$$

(60)

where $\overline{\sigma }$ is the maximum singular value $\overline{\sigma }\left(H(j\omega )\right)$ of the matrix $H(j\omega )$.

The system represented by Equation (18) can be considered as a SISO system with additional forces from installed TMDs. The form of the objective function proposed for the ℌ₂ optimization issue will have the following form:

$$\min J1={\Vert H\Vert }_{{\mathscr{H}}_2}={\left(\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{\infty }{\left|\overline{H}{}_{jj}(j\omega )\right|}^2}d\omega \right)}^{1/2}.$$

(61)

An objective function associated with the ℌ_∞ optimization issue was also proposed. In this case, the optimization criterion is the minimization of the maximum values of the FRF $\left|H(j\omega )\right|$ (see Figure 3), which can be described by the formula:

$$\min J2={\Vert H\Vert }_{{\mathscr{H}}_{\infty }}=\underset{\omega \in ℝ}{\sup }\left|\overline{H}{}_{jj}(j\omega )\right|.$$

(62)

Both cases of the FRF form and optimization included the FRF of the equivalent system $\left|\overline{H}{}_{jj}(j\omega )\right|$, associated with the j degree of freedom, which was adopted as the equivalent system reduction point, and with the i-th mode of structural natural vibrations and the corresponding value of modal displacement, associated with the primary ordinate ${\overline{q}}_{ji}(t)=a_{ji}{\psi }_i(t)$.

2.3. The Issue of Optimization

The analysis used two optimization methods in the form of a standard genetic algorithm (GA), and a simulated annealing (SA) algorithm. Ready libraries for the DELPHI software environment were used.

The procedure for the numerical determination of the optimal parameters can be summed up in the following steps:

• Calculations of the eigenvalues and eigenvectors in Finite Element (FEM) software, adopting baseline frequency ${\omega }_i$, to which MTMDs will be tuned;
• Adopting a point for reducing a system to an equivalent SDOF system, normalization of eigenvectors and determination of the parameters for the Equation (7);
• Adopting a TMD number and initial parameters $m_k^t,c_k^t,k_k^t$ for the TMDs;
• Adopting variables, which will be subject to the optimization process, e.g., $c_k^t,k_k^t$ for a constant value of $m_k^t$ (different constant parameter options possible);
• Selection of the optimization method (GA or SA) and the optimization issue, calculations of the equivalent system FRF $\left|\overline{H}{}_{jj}(j\omega )\right|$, and the determination of values using the Equations (61) or (62) at each optimization step;
• After reaching the desired calculation accuracy (change of the J1 or J2) value, saving the obtained MTMD parameters, and calculations of the equivalent system FRF $\left|\overline{H}{}_{jj}(j\omega )\right|$ associated with the j degree of freedom, which was adopted as the equivalent system reduction point.

3. Numerical Example 1

3.1. Input Data

The analysis adopted a cantilever structure in the form of a chimney made of S235 steel, with a height of h=160 m. The chimney parameters were selected so as to obtain the first circular natural frequency equal to ${\omega }_1=1\mathrm{rad}/\mathrm{s}$. The chimney diameter and its wall thickness were adopted as a constant, along its entire height. These values equal to, respectively, D = 4010 mm, g = 25.5 mm. The cantilever structure was divided into 40 beam elements with the following parameters: Young modulus E = 210 GPa, Poisson coefficient ν = 0.3, and density ρ = 7850 kg/m³. First three natural vibration modes and their corresponding eigenvalues were calculated. FE software was used to import a normalized eigenvector for the first mode of the natural vibrations, which was then used to construct motion equations for the structure reduced to SDOF system with attached TMD.

Figure 4 shows a FEM model with the numberings of the bars and modes, and a presentation of the first three mode shapes and normalized eigenvectors.

The analysis adopted the issue of tuning TMD and MTMDs to the first natural frequency ${\omega }_1=1\mathrm{rad}/\mathrm{s}$. For this frequency, the adopted structural damping was in the form of a damping ratio ${\zeta }_1=0.02$. The chimney tip (node 41) was selected as the S point for reducing the structure to an SDOF equivalent system. The problem involved finding optimal parameters for a single TMD and MTMDs in the form of 2, 4, 8, and 20 TMDs located also at the place of the highest ordinate of the first vibration mode, which was node 41. New $\overline{M},\overline{C},\overline{K}$ matrices for the equivalent system were determined based on the presented modal reduction equations. The equivalent mass determined from the formula ${\overline{m}}_{ji}={\tilde{m}}_i/a_{ji}{}^2=1/a_{ji}{}^2$ (where j is the system point of reduction, namely, node 41, whereas i = 1 corresponds to the first natural vibration mode) amounted to ${\overline{m}}_{41,1}={\tilde{m}}_1/a_{41,1}{}^2={10}^5\mathrm{kg}$, equivalent stiffness was ${\overline{k}}_{41,1}={\tilde{k}}_1/a_{41,1}{}^2=100\mathrm{kN}/\mathrm{m}$, and equivalent damping was ${\overline{c}}_{41,1}={\tilde{c}}_1/a_{41,1}{}^2=2{\zeta }_1{\omega }_1/a_{41,1}{}^2=4\times {10}^3\mathrm{kg}/\mathrm{s}$. Numerical optimization was conducted for all cases. For each case, the total MTMD mass was equal to the mass of the single TMD, which was adopted as equal to $m_1{}^t=2\times {10}^3\mathrm{kg}$. Therefore, the mass factor was $\mu =m_1{}^t/{\overline{m}}_{41,1}=0.02$. For all other cases, TMD masses were determined from the formula $m_k{}^t=\mu \cdot {\overline{m}}_{41,1}/N_k$. The natural circular frequency for a single TMD was determined from the formula ${\omega }_k{}^t=\sqrt{k_k{}^t/m_k{}^t}$, whereas the damping factor values for individual TMDs were determined from the formula $c_k{}^t=2{\zeta }_k{}^t{m}_k{}^t{\omega }_k$. Stiffness $k_k{}^t$ was used to determine the dimensionless tuning ratio for a single TMD ${\beta }_k={\omega }_k{}^t/{\omega }_1$. Initial tuning adopted for all dampers within the optimization amounted to ${\beta }_k=1$. This constituted a base to adopt the starting values $k_k{}^t$.

Table 1. Mass values for individual TMDs and $k_k{}^t$ and ${\zeta }_k{}^t$ starting values.

TMD Number	${\mathit{m}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{k}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\zeta }}_{\mathit{k}}{}^{\mathit{t}}$
	[kg]	[kN/m]	[–]
1 TMD	2000	2.00	0.05
2 TMDs	1000	1.00	0.05
4 TMDs	500	0.50	0.05
8 TMDs	250	0.25	0.05
20 TMDs	100	0.10	0.05

shows the adopted constant masses of a single TMD and the adopted starting values of the $k_k{}^t$ and ${\zeta }_k{}^t$ MTMD parameters (values subject to optimization), depending on the TMD number.

The calculation process involved using the aforementioned GA and SA optimization methods and the objective functions based on ℌ₂ and ℌ_∞ norms in the form of J1 and J2. The frequency range when determining the objective function was (0.0 rad/s:π rad/s).

3.2. Numerical Optimization Results—Comparison of the GA and SA Heuristic Methods

Both heuristic methods do not guarantee finding an exact solution, which is the consequence of the very specificity of the methods. The optimization procedures GA and SA were validated on the 1TMD, 2TMD, and 4TMD models. The parameters of stiffness and damping ratio obtained from both GA and SA optimization methods did not differ. Additionally, both methods were compared for the case of 10TMD (Zuo and Nayfeh [20]) attached to an undamped primary system. Figure 5 shows the optimal TMD parameter values obtained using both methods. Compared to SA, the GA method did not provide results convergent with the results presented in the work by Zuo and Nayfeh [20], while the obtained values of TMDs’ parameters form SA were completely convergent. Unfortunately, with a greater number of optimization variables, the GA algorithm was less stable than the SA algorithm. The obtained results are discussed in greater detail in Section 5. Only the SA optimization method was applied in the further part of the optimization calculations.

3.3. Numerical Optimization Results—SA Method

Optimal MTMD parameters, which are discussed below, were obtained as a result of computations.

Table 2. The value of the J1 and J2 functions after completed optimization with the SA method.

TMD Number	J1 Value	J1 Change	J2 Value	J2 Change
	[m/(kN*s^0.5)]	[%]	[m/kN]	[%]
SDOF	2.4910 × 10−2		2.5005 × 10−1
1 TMD	1.5515 × 10−2	0.00%	7.4579 × 10−2	0.00%
2 TMDs	1.5203 × 10−2	−2.01%	6.8248 × 10−2	−8.49%
4 TMDs	1.5001 × 10−2	−3.31%	6.4091 × 10−2	−14.06%
8 TMDs	1.4875 × 10−2	−4.12%	6.1620 × 10−2	−17.38%
20 TMDs	1.4783 × 10−2	−4.72%	6.0202 × 10−2	−19.28%

shows the end values of the J1 and J2 objective functions, after the completion of the SA optimization process, using the J1 and J2 objective functions. In addition, it also specifies the percentage change of the objective function value, relative to a 1TMD system. Computation accuracy was adopted at a level of 10⁻¹⁰.

Table 3. Optimal $k_k{}^t,{\zeta }_k{}^t,{\beta }_k$ values obtained for the J1 objective function.

MTMD	No.	${\mathit{m}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{k}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\zeta }}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\beta }}_{\mathit{k}}$
		[kg]	[kN/m]	[–]	[–]
1 TMD	1	2000	1.960905	0.068462	0.990178
2 TMDs	1	1000	0.895498	0.041605	0.946307
	2	1000	1.096614	0.044381	1.047193
4 TMDs	1	500	0.419478	0.025349	0.915945
	2	500	0.470454	0.025186	0.970004
	3	500	0.526416	0.025998	1.026076
	4	500	0.595890	0.028082	1.091687
8 TMDs	1	250	0.200353	0.016536	0.895215
	2	250	0.214738	0.014990	0.926796
	3	250	0.228431	0.014593	0.955889
	4	250	0.242203	0.014669	0.984282
	5	250	0.256858	0.015086	1.013624
	6	250	0.273079	0.015356	1.045139
	7	250	0.291974	0.015877	1.080691
	8	250	0.316299	0.017789	1.124809

and

Table 4. Optimal $k_k{}^t,{\zeta }_k{}^t,{\beta }_k$ values obtained for the J2 objective function.

MTMD	No.	${\mathit{m}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{k}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\zeta }}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\beta }}_{\mathit{k}}$
		[kg]	[kN/m]	[–]	[–]
1 TMD	1	2000	1.905820	0.089169	0.976171
2 TMDs	1	1000	0.878697	0.056278	0.937388
	2	1000	1.046257	0.062921	1.022867
4 TMDs	1	500	0.413107	0.036608	0.908964
	2	500	0.456409	0.038155	0.955415
	3	500	0.503223	0.039621	1.003218
	4	500	0.560626	0.041295	1.058892
8 TMDs	1	250	0.198173	0.023339	0.890332
	2	250	0.210676	0.024023	0.917988
	3	250	0.222606	0.025988	0.943624
	4	250	0.233932	0.028452	0.967329
	5	250	0.245215	0.028899	0.990383
	6	250	0.258339	0.027081	1.016542
	7	250	0.274253	0.026689	1.047383
	8	250	0.293970	0.026947	1.084380

show optimal MTMD parameters obtained with optimized J1 and J2 objective functions. A tabular presentation of the 20 TMDs was omitted due to the high number of the results. Results for these cases are shown as drawings in the next section. Figure 6a shows a graph of the FRF module for an equivalent system $\left|{\overline{H}}_{41,41}\right|$, for a different number of TMDs, and with applied J1 objective function, while Figure 6b shows the same for the J2 objective function.

3.4. Impact of Primary System Mistuning

The analysis also covered the impact by the changing mass of the equivalent structure ${\overline{m}}_{41,1}$, on the value of the equivalent structure FRF $\left|{\overline{H}}_{41,41}\right|$ for a structure with attached MTMDs with optimal parameters, predetermined based on the J1 and J2 objective functions. A 10% reduction and increase in the equivalent structure mass was assumed, which resulted in the change of the natural circular frequencies to a value of ${\omega }_1{}^{-10\%}=1.05409 \mathrm{rad}/\mathrm{s}$ and ${\omega }_1{}^{+10\%}=0.95346 \mathrm{rad}/\mathrm{s}$, respectively. Figure 7a shows a graph of the FRF module for an equivalent system $\left|{\overline{H}}_{41,41}\right|$, for a mass reduction of 10% and a different number of TMDs, and with applied J1 objective function, while Figure 7b shows the same for the J2 objective function. Figure 8a,b shows analogous graphs but for a mass increase of 10% (designation +10% and −10% means an increase and reduction of the ${\overline{m}}_{41,1}$ equivalent mass, respectively).

3.5. Discussion

The analyses presented in Section 3 were used as a base to compare the application of the SA optimization method and the J1 and J2 objective functions with known analytical solutions in terms of a single TMD.

Table 5. Optimal parameters and J1 and J2 objective function values for 1TMD, based on various authors (values for $\mu =0.02$).

	$\mathbf{Tuning}$ ${\mathit{\beta }}_{\mathit{o}\mathit{p}\mathit{t}}$	$\mathbf{Damping}$ ${\mathit{\zeta }}_{\mathit{o}\mathit{p}\mathit{t}}$	J1 Value ${\mathit{\zeta }}_1=0$	J2 Value ${\mathit{\zeta }}_1=0$	J1 Value ${\mathit{\zeta }}_1=0.02$	J2 Value ${\mathit{\zeta }}_1=0.02$
			[m/(kN*s0.5]	[m/kN]	[m/(kN*s0.5]	[m/kN]
Den Hartog	$\frac{1}{1+\mu }$ 0.9803922	$\sqrt{\frac{3\mu }{8(1+\mu )}}$ 0.0870388	1.886 × 10−2	1.005 × 10−1	1.559 × 10−2	7.676 × 10−2
Warburton (1981)	$\frac{1}{1+\mu }\sqrt{\frac{2+\mu }{2}}$ 0.9852819	$\sqrt{\frac{\mu (4+3\mu )}{8(1+\mu )(2+\mu )}}$ 0.0701871	1.866 × 10−2	1.091 × 10−1	1.550 × 10−2	8.176 × 10−2
Warburton (1982)	$\frac{1}{1+\mu }\sqrt{\frac{2-\mu }{2}}$ 0.9754779	$\sqrt{\frac{3\mu }{8(1+\mu )(1-\mu /2)}}$ 0.0861813	1.891 × 10−2	1.057 × 10−1	1.560 × 10−2	7.497 × 10−2
Ren	$\sqrt{\frac{1}{1-\mu }}$ 1.0101525	$\sqrt{\frac{3\mu }{8(1-\mu /2)}}$ 0.0870388	1.941 × 10−2	1.366 × 10−1	-	-
Nishihara and Matsuhisa	$\frac{1}{1+\mu }\left(1-{\zeta }_1\sqrt{\frac{\mu }{1+\mu -{\zeta }_1{}^2}}\right)$ 0.9776460	$\frac{{\zeta }_1+\sqrt{\mu \left(1+\mu -{\zeta }_1{}^2\right)}}{1+\mu }$ 0.1596084	-	-	1.677 × 10−2	9.865 × 10−2
Proposed J1 ${\zeta }_1=0$	- 0.9890635	- 0.0682140	1.885 × 10−2	1.148 × 10−2	-	-
Proposed J2 ${\zeta }_1=0$	- 0.9803647	- 0.0853719	1.868 × 10−2	1.006 × 10−2	-	-
Proposed J1 ${\zeta }_1=0.02$	- 0.9901779	- 0.0684623	-	-	1.563 × 10−2	8.521 × 10−2
Proposed J2 ${\zeta }_1=0.02$	- 0.9761710	- 0.0891687	-	-	1.552 × 10−2	7.458 × 10−2

shows known formulas for the determination of the optimal parameters for TMD usually determined from a TDOF system, with stated values of the optimal parameters. The table also presents original results for 1 TMD, taking into account natural damping of the primary system, as well as without it. Figure 8 shows an equivalent system FRF module graph for optimal TMD parameters. When analyzing the obtained solution (see Table 5 and Figure 9a), in the case of no damping of the primary system ${\zeta }_1=0$ and in comparison with the Den Hartog solution [1], which has been known for years, one can state very good conformity of the optimal tuning parameters ${\beta }_1$ and the TMD damping factor ${\zeta }_1{}^t$, whereas optimal TMD parameters proposed by Warburton [2], upon a harmonic excitation, are close to the parameters obtained through calculations of the J2 objective function optimization. The solution proposed by Ren [31] suits a different kind of TMD, where the damping element is not connected with the primary mass but with the substrate (other motion equations); hence, it is impossible to directly apply the optimal parameter formulas proposed in this study as parameters for a standard TDOF system. Furthermore, please note that optimal TMD parameters based on the ℌ₂ norm (white noise input), i.e., determined from the J1 objective function, provided similar TMD efficiency, as suggested by Warburton in [3].

In the case of an analysis covering a structure including natural damping of the primary system, we can observe that FRF module graphs, developed based on optimal Den Hartog and Warburton parameters, are no longer consistent with the graphs based on the parameters obtained through the ℌ₂ and ℌ_∞ methods, that is, the J1 and J2 objective functions. This, of course, stems from the fact that these formulas do not include natural damping ${\zeta }_1$, whereas, when comparing to the Matsuhis solution [32], where the optimal parameters were obtained based on the stability criterion, taking into account ${\zeta }_1$, it should be concluded that his solution was optimal neither for harmonic excitations nor for white noise input.

In [20], Zuo and Nayfeh conducted a free analysis of the ${\beta }_k,{\zeta }_k{}^t$ parameters, with a constant value of $m_k{}^t$, for 10 and 100 TMDs, with not damping of the primary system ${\zeta }_1=0$ and a $\mu =0.05$ mass ratio. They applied optimization based on the ℌ₂ norm and used state space equations for an SDOF system with attached MTMDs. The work also presents tuning increment graphs for individual TMDs in the form of ${\beta }_k-{\beta }_{k-1}$. The authors of that article also conducted such simulations for n = 10 (number of TMD), ${\zeta }_1=0$, and $\mu =0.05$. Figure 10 shows tuning ${\beta }_k$ and damping ratio ${\zeta }_k$ distribution based on the number of TMDs, while Figure 11 shows tuning increments and tuning factor for individual TMDs in the form of ${\beta }_k-{\beta }_{k-1}$ and ${\zeta }_k-{\zeta }_{k-1}$. The obtained ${\beta }_k,{\zeta }_k{}^t$ and ${\beta }_k-{\beta }_{k-1}$ parameters exhibited very good conformity with the results of Zuo and Nayfeh [20].

In the course of analyzing the results shown in Section 3, which concern the optimization of MTMDs located at the tip of a h = 160-m-high chimney, higher effectiveness of MTMDs relative to a single TMD can be confirmed (see Table 2 and Figure 6). This is particularly clear to optimal parameters determined with the J2 objective function, where the location of the flat section of graph $\left|{\overline{H}}_{41,41}\right|$ in the area of ${\omega }_1$ can be seen increasingly lower (Figure 6b). Another issue that needs to be stressed are the graphs $\left|{\overline{H}}_{41,41}\right|$ (Figure 6a) for TMD tuning in the form of white noise ℌ₂ optimization of the J1 function exhibit a broader damper impact range on the side of higher input frequencies.

Furthermore, Figure 12 shows a function variability graph depending on the number of TMDs. Of course, higher MTMD efficiency relative to TMDs can be observed for both objective functions, J1 and J2, applied within the numerical optimization. Furthermore, Figure 12 shows the values of individual objective functions for the optimal parameters obtained using both optimization methods, namely, ℌ₂ (J1 function) and ℌ_∞ (J2 function). It can be seen (see Figure 12a) that in the case of random excitation, which corresponds to the value of the J1 function and TMDs tuned to harmonic excitation (optimal parameters based on the optimization of function J2), we get a lower MTMD efficiency. Figure 12b indicates an inverse relationship.

When analyzing the results concerning the primary system mass change (see Figure 7 and Figure 8), for an approximately 5% difference in the value of ${\omega }_1$ (“structural mistuning”), MTMDs are more effective than a single TMD. A lower effectiveness with a higher number of TMDs tuned to harmonic excitations can be observed only in the case of TMDs tuned based on J2 optimization and reducing the primary system mass increase in frequency ${\omega }_1{}^{-10\%}=1.05409 \mathrm{rad}/\mathrm{s}$. In the event of such a primary system mistuning, better effectiveness was exhibited by MTMDs tuned based on J1, which resulted from the aforementioned broadband excitation on the side of frequencies higher than ${\omega }_1$. It can also be concluded that a reduction in the primary system mass had a more adverse influence on the structural response and MTMD operation than mass increase. These relationships are shown in Figure 13. It can also be assumed that with a low difference (below 5%) between the calculated natural frequencies of the structure and the ones determined for a real structure, it is advisable to use MTMDs instead of a single TMD. The aforementioned analyses confirm the general conclusions included in the works by Zuo and Nayfeh [20] as well as Li and Ni [21] concerning the higher effectiveness of MTMDs relative to TMDs.

Additionally, for analysis purposes, it was decided to present the results of optimal MTMD parameters in the form of optimal tuning ratio ${\beta }_k$ and optimal damping ratio ${\zeta }_k{}^t$ graphs depending on the number of TMDs. The tuning and damping factor increments for individual TMDs in the form of ${\beta }_k-{\beta }_{k-1}$ and ${\zeta }_k{}^t-{\zeta }_{k-1}^t$ were also presented. After a thorough review of the graphs (see Figure 14, Figure 15 and Figure 16) presented for cases with 4, 8, and 20 TMDs, conclusions can be drawn regarding the application of the J1 objective function based on the ℌ₂ norm and the J2 function based on the ℌ_∞ norm. As far as the tuning distribution is concerned, we can observe increased factor ${\beta }_k$ for TMDs tuned above and below the ${\omega }_1$ frequency. The area of TMDs tuned to values beyond ${\omega }_1$ exhibited a significantly higher growth of the value ${\beta }_k$, especially for the parameters obtained from the J1 objective function, which is particularly well illustrated by the graphs ${\beta }_k-{\beta }_{k-1}$. If we consider the optimal tuning ratio graphs for individual ${\zeta }_k{}^t$, we can observe a derivative characteristic, but only for the J1 function, although the increments ${\zeta }_k{}^t-{\zeta }_{k-1}^t$ no longer show a similar dependency as ${\beta }_k-{\beta }_{k-1}$. On the other hand, the optimal parameters ${\zeta }_k{}^t$ obtained based on the J2 function look completely different. In this case, the graph shows a significant drop in the ${\zeta }_k{}^t$ values for “central” TMDs tuned on ${\omega }_1$.

4. Numerical Example 2

4.1. Input Data

The next analysis is an example of the optimization of MTMD mounted on a complex, reinforced concrete structure in the form of a tall building. The analysis included Gray Office “A” located in the city of Lublin in Poland. The building is set on a foundation slab with two underground and 13 aboveground stories. The heights of the underground stories are 3.5 m and the overground stories are 3.9 m. The total height of the building is 56.55 m. Figure 17a shows the architectural visualization of the building, while Figure 17b shows the structural plans of floors 3 and 13. The main load-bearing system is the reinforced concrete core of the building (modified for the present analysis by adding walls in the E-H axes), with wall thicknesses from 20 cm to 65 cm. The remaining load-bearing structure is a reinforced concrete slab ceiling (ceilings thickness from 25 cm to 45 cm), based on reinforced concrete pillars of various dimensions from 60–100 cm.

Figure 18 shows the FEM model of the system. The used beam and shell elements with concrete parameters C25/C30 are Young modulus E = 32 GPa, Poisson coefficient ν = 0.2, and density ρ = 2870 kg/m3 (increasing the density to accommodate the flooring layers). The number of nodes in the model was 11,894, and the total number of equations in the modal analysis was 69,912.

The first 11 corresponding natural vibration modes and their eigenvalues were calculated.

Table 6. The compilation of dynamic characteristics of structure.

Item No.	Frequency [Hz]	Period of Vibrations [s]	Form of Vibration
1	0.8431	1.18607	flexural in the Y direction
2	1.3476	0.74207	torsional
3	1.4488	0.69022	torsional
4	2.5197	0.39688	local flexural form of cantilevers
5	2.5222	0.39647	local flexural form of cantilevers
6	2.9290	0.34141	local flexural form of cantilevers
7	2.9361	0.34059	local flexural form of cantilevers
8	2.9372	0.34046	local flexural form of cantilevers
9	4.0457	0.24717	flexural in the Y direction
10	4.3516	0.22980	torsional
11	4.8888	0.20455	flexural in the Z direction

presents the description of particular modes of vibration below the value of f < 5 Hz, while Figure 19 shows selected modes’ shapes, disregarding the local natural frequencies (f₄ to f₈) of the cantilever beams located on the last floor of the building. The frequencies of natural vibrations were separated. The first flexural form of natural vibrations in the Y direction corresponded to f₁ = 0.8431 Hz, and the next bending form of vibrations occurred at the frequency f₂ = 1.3476 Hz. FE software was used to import a normalized eigenvector for the first mode of the natural vibrations, which was then used to construct motion equations for the structure reduced to SDOF system with attached TMDs.

The analysis adopted the issue of tuning MTMD to the first natural frequency f₁ = 0.8431 Hz. The damping value was assumed in the form of a logarithmic damping decrement δ = 0.08. Node 11,110 was accepted was selected as the S point for reducing the structure to an SDOF equivalent system. The problem involved finding optimal parameters for MTMDs in the form of 2, 4, and 6 TMDs located on the ceiling of the top floor. There were three cases of TMDs attachment: 2 TMDs, at nodes 11,110 and 11,111; 4 TMDs, at nodes 11,243, 10,977, 11,222, and 11,004; 6 TMDs, at nodes 11,243, 11,110, 10,977, 11,222, 11,111, and 11004. The Y direction was adopted as the operating direction of all TMDs. The designations of the nodes are shown in Figure 20.

New $\overline{M},\overline{C},\overline{K}$ matrices for the equivalent system were determined based on the presented modal reduction equations. The equivalent mass determined from the formula ${\overline{m}}_{ji}={\tilde{m}}_i/a_{ji}{}^2=1/a_{ji}{}^2$ (where j is the system point of reduction, namely, node 11,110, whereas i = 1 corresponds to the first natural vibration mode) amounted to ${\overline{m}}_{11110,1}={\tilde{m}}_1/a_{11110,1}^2=6.7449\times {10}^6\mathrm{kg}$, equivalent stiffness was ${\overline{k}}_{11110,1}={\tilde{k}}_1/a_{11110,1}^2={\omega }_1{}^2/a_{11110,1}^2=1.8909\times {10}^5\mathrm{kN}/\mathrm{m}$, and equivalent damping was ${\overline{c}}_{11110,1}^{}={\tilde{c}}_1/a_{11110,1}^2=2{\zeta }_1{\omega }_1/a_{11110,1}^2=9.09\times {10}^2\mathrm{kg}/\mathrm{s}$. Numerical optimization was conducted for all cases. For each case, the total MTMD stiffness was equal to the stiffness of the 1 TMD, which was adopted as equal to $k_1{}^t=120\mathrm{kN}/\mathrm{m}$. For all other cases, TMDs’ stiffness was determined from the formula $k_k{}^t=k_1{}^t/N_k$. In the optimization process, the variables were the mass $m_k{}^t$ and the damping ratio ${\zeta }_k{}^t$.

The calculation process involved using the aforementioned SA optimization methods and the objective functions based on ℌ₂ and ℌ_∞ norms in the form of J1 and J2.

4.2. Numerical Optimization Results and Discussion

Optimal MTMD parameters were obtained as a result of computations.

Table 7. Optimal $m_k{}^t,{\zeta }_k{}^t,{\beta }_k$ values obtained for the J1 objective function.

MTMD	No.	${\mathit{m}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{k}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\zeta }}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\beta }}_{\mathit{k}}$
		[103 kg]	[kN/m]	[–]	[–]
2 TMDs	1	2.09840	60	0.007908	1.009423
	2	2.18777	60	0.007888	0.988589
4 TMDs	1	1.03307	30	0.004949	1.017274
	2	1.05885	30	0.004739	1.004811
	3	1.08320	30	0.004698	0.993454
	4	1.11004	30	0.005010	0.981369
6 TMDs	1	0.68397	20	0.003796	1.020795
	2	0.69665	20	0.003459	1.011457
	3	0.70786	20	0.003420	1.003420
	4	0.71903	20	0.003489	0.995597
	5	0.73084	20	0.003547	0.987518
	6	0.74467	20	0.003926	0.978307

and

Table 8. Optimal $m_k{}^t,{\zeta }_k{}^t,{\beta }_k$ values obtained for the J2 objective function.

MTMD	No.	${\mathit{m}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{k}}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\zeta }}_{\mathit{k}}{}^{\mathit{t}}$	${\mathit{\beta }}_{\mathit{k}}$
		[103 kg]	[kN/m]	[–]	[–]
2 TMDs	1	2.10516	60	0.011970	1.007800
	2	2.18180	60	0.011816	0.989942
4 TMDs	1	1.03823	30	0.008284	1.014739
	2	1.06098	30	0.007508	1.003804
	3	1.08231	30	0.007402	0.993861
	4	1.10542	30	0.007627	0.983416
6 TMDs	1	0.68809	20	0.006004	1.017729
	2	0.69927	20	0.005902	1.009561
	3	0.70923	20	0.005837	1.002452
	4	0.71897	20	0.005787	0.995633
	5	0.72918	20	0.005922	0.988639
	6	0.74093	20	0.006122	0.980767

show optimal MTMD parameters obtained with optimized J1 and J2 objective functions. Figure 21 shows plots of the FRF module for node 11,110 for X and Y direction, respectively (for original system with a different number of TMDs and J1 optimization), while Figure 22 shows the same FRF module plots for the J2 objective function. From the analysis of plots of the FRF module (Figure 21 and Figure 22), it can be seen that for the X direction, we did not observe changes in the FRF module value, while for the Y direction, we observed a decrease in the FRF module value in the first natural frequency range f₁ = 0.8431 Hz. Additionally, we can observe that the vibration frequencies f₄ to f₈ are not visible on the FRF module plots.

Figure 23 shows the FRF module in the Y direction for two nodes, 11,110 and 11,111, in the frequency range 0.8 Hz to 0.9 Hz. As in the previous example, the same nature of the FRF module plots for the optimal TMD parameters obtained from the functions J1 and J2 can be observed. However, it should be borne in mind that the optimization process was carried out on the equivalent displacement of the SDOF system, while the FRF module plots were presented for the original system with thousands of DOF. The presented graphs also show the fact of greater efficiency of MTMD with an increasing number of TMDs.

5. Conclusions

The adopted novel method for developing motion equations enables adding single TMDs or MTMDs to completely different degrees of freedom of the primary system. The system of equations allows for easy MTMD tuning for complex vibrations modes, with MTMDs located in local maxima of these vibrations’ modes, while still analyzing the SDOF system with attached MTMDs. Locating MTMDs at the tip of the chimney provided a possibility to compare the results with studies on MTMD tuning. The next stage of the research was to optimize MTMD parameters in a complex structure, but still analyze as an SDOF system. Of course, the thesis of better effectiveness of MTMD over a single TMD was confirmed.

Based on conducted analysis, it was observed that the $k_k{}^t,{\zeta }_k{}^t$ free optimization was of completely different nature in the case of optimizations based on the ℌ₂ and ℌ_∞ norms, especially for the ${\zeta }_k{}^t$ parameters, where the ${\zeta }_k{}^t-{\zeta }_{k-1}^t$ increments showed totally different forms for both J1 and J2 objective functions.