Multi-Objective Design Optimization of Multiple Tuned Mass Dampers for Attenuation of Dynamic Aeroelastic Response of Aerospace Structures

Abstract

This paper proposes a design procedure to determine the optimal configuration of multi-degrees of freedom (MDOF) multiple tuned mass dampers (MTMD) to mitigate the global dynamic aeroelastic response of aerospace structures. The computation of the aerodynamic excitations is performed considering two models of atmospheric disturbances, namely, the Power Spectral Density (PSD) modelled with the Davenport Spectrum (DS) and the Tuned Discrete Gust (TDG) with the one-minus cosine profile. In order to determine the optimum sets of MTMD, a Multi-objective design Optimization considering Genetic Algorithm (MOGA) is implemented, where the selected fitness functions for the analysis are the minimization of the total mass of the resonators as well as the concurrent minimization of the peak displacements of a specified structural node in all translational degrees of freedom. A case study is presented to demonstrate the proposed methodology, where the optimal sets of MTMD are determined for the concurrent minimization of the pointing error of a truss-like antenna structure as well as the mass of the considered MTMD. It is found that the placement of the MTMD in the primary reflector of the antenna structure provided a maximum reduction in the pointing error of 62.0% and 39.2%, considering the PSD and the TDG models, respectively. Finally, this paper presents an advanced framework to estimate optimal parameters of MTMD control devices under convoluted loading cases as an initial step towards the use of such passive systems in applications that commonly employ active or semi-active solutions.

1. Introduction

Aerospace structures are subjected to a wide range of dynamic loadings. When excited close to resonance frequencies, the dynamic loads can cause excessive vibrations in the structures, resulting in loss of accuracy [1,2], structural instabilities [3,4], and material fatigue [5,6], among others, causing potential catastrophic failures. Additionally, it is always found that high levels of vibration responses normally occur in aerospace structures when subjected to aerodynamic loadings due to their inherent slenderness and lightweight design considerations [7,8]. On the other hand, in the last few decades, there has been a significant increase in the stringent design and certification requirements for aerospace structures, such that more efficient design solutions related to material-to-strength ratio and vibration control systems are implemented, aiming to improve the operation quality of those structures [9,10,11,12].

In this context, one of the major goals in the design of aerospace structures is to increase their efficiency regarding structural fatigue and dynamic stability. A wide range of vibration control devices have been proposed in the literature, characterized by passive, adaptive, semi-active [13,14], or active nature [15,16], and can be effectively implemented to achieve those design requirements. Particularly, passive vibration attenuation devices are based on the principle that a non-active system is installed in a structure, being capable of transferring or absorbing part of the kinetic energy of global mode shapes and not requiring an external power source to perform its function. A specific type of passive vibration attenuation device is the tuned mass damper (TMD), which essentially consists of an auxiliary mass attached to the main structure. TMDs can be classified into four categories, namely, the tuned liquid column dampers [17], pendular TMDs [18], bidirectional and homogeneous TMDs [19], and conventional TMDs. Due to the ease of their implementation and robust design methodologies, they have been extensively used in civil engineering structures, such as high-rise buildings [20,21,22], telecommunication towers [23], wind turbines [24,25], or bridges [26,27] in order to mitigate seismic, hydrodynamic, and aerodynamic vibrations. TMDs are also used in aircraft wing structures for aeroelastic response reduction [28]. Recently, Lee et al. [29] designed a conventional TMD system to mitigate the elastic response of a large retractable truss-like frame under seismic load excitations, achieving a 90% reduction in the vibration levels with the inclusion of 8 TMDs when compared to the uncontrolled structure. Liu et al. [30] used a single degree of freedom (SDOF) TMD to mitigate the vibration response of force-aft wind turbine towers, comparing the analytical formulation with a radial basis function neural network to achieve the optimum solution.

The conventional design of a single TMD, similarly to a vibration neutralizer, consists of tuning the device to the fundamental natural frequency of the structure, i.e., the first mode of vibration [31,32]. This approach may not be entirely suitable for aerospace structures due to the complex interaction between the aerodynamic loads and the elastic structure, which could return different optimal values of design parameters according to the established operational conditions [33]. Recently, Ma et al. [34] investigated the efficiency of MDOF TDM applied in beam structures, showing that with such a configuration, it is possible to attenuate multiple modes of interest while requiring less mounting space than conventional SDOF TMD systems. Additionally, the use of MDOF MTMD devices is an approach that overcomes the limitations inherent to the use of a single TMD, due to the fact that the auxiliary masses are distributed throughout the target structure at designated locations, therefore permitting the control of upper vibration modes in complex shape structures [35,36]. The performance of MTMD control devices depends intrinsically on the design parameters of mass, stiffness, and damping, as well as the number and location of the auxiliary masses, thus defining the covering frequency range of response attenuation [37]. In order to determine the best sets of design parameters, it is important to implement an optimization algorithm that could minimize physical objective functions such as displacement, velocity, or acceleration responses. In the case of aerospace structures, another critical factor that needs to be considered by the designer is the minimization of the total mass inclusion when using the MTMD as the main control devices. This way, it is pertinent to develop a reliable optimization framework for determining optimal design parameters of MTMD applied to aerospace structures, under the consideration of realistic aerodynamic load scenarios, such that the performance of the control devices can be effectively quantified.

Hence, this work proposes a multi-objective design optimization framework to find configurations of MDOF MTMD that are capable of attenuating the dynamic aeroelastic response in aerospace structures under aerodynamic gust excitations while concurrently minimizing the total mass inclusion of the vibration control devices. For this, we developed an in-house routine that integrates the auxiliary masses with a dynamic aeroelastic response analysis simulation that considers two models of aerodynamic excitation profiles, the PSD and the TDG, with DS and one-minus-cosine profiles, respectively. As for the optimization algorithm selection, due to the inherent complexity of the problem, we used the metaheuristic MOGA [38] to determine the best MTMD individuals. Intending to test the proposed method, we performed a thorough analysis using a case study of a truss-like antenna structure under different operational gust cases. For all cases analyzed, three configurations of MTMD placement in the antenna were considered: in the primary reflector truss, in the boom arms of the secondary reflector, and in both locations combined.

The paper is organized as follows: After this introduction, the methodology developed is presented, where the computational model of the MTMD is developed and discussed, followed by the dynamic aeroelastic analysis and the MOGA framework. In sequence, the application case study is investigated, where the details of the geometry are defined alongside the pointing accuracy formulation, and the gust operational cases considered. Then, the results for the inclusion of MTMD for the different placement configurations are presented and discussed. The last section presents the conclusions of the study.

2. Materials and Methods

In this section, the formulation of the MTMD is presented in the global equation of motion of the structural FE model, considering the mass, stiffness, and damping matrices. In sequence, the aeroelastic equation of motion is derived with the modified FE constitutive matrices, considering the gust models given by the DS and the TDG. Then, the general formulation of the Time Consistent Loads and Displacements (TCLD) is proposed, such that different degrees of freedom (DOF) can be tuned in phase to account for time-consistent displacements when analyzing a response signal in the frequency domain. Finally, the optimization procedure based on a MOGA is introduced, where the chosen design variables, constraints, and fitness functions are discussed.

2.1. Structural Model with Multiple-Tuned Mass Dampers

Consider the FE model of a generalized structure with N number of DOF, where each degree of freedom can receive up to K number of attached TMDs in the form of an auxiliary MTMD configuration, as depicted in Figure 1. Here, each auxiliary control device included can be modelled as local inclusions of mass, stiffness, and damping values at the generalized structural nodes, which modifies the dynamic global constitutive matrices of the FE model. The equation of motion of the system shown in Figure 1 can then be written as [39]:

$$\tilde{M}\ddot{u}\left(t\right)+\tilde{C}\dot{u}\left(t\right)+\tilde{K}\mathbf{u}\left(t\right)=\tilde{\mathbf{P}}\left(t\right)$$

(1)

where $\tilde{\mathbf{M}}$, $\tilde{\mathbf{C}}$, and $\tilde{\mathbf{K}}$ are, respectively, the constitutive global FE matrices of mass, damping, and stiffness of the system modified by the inclusion of the MTMD, and $u\left(t\right)=\mathrm{col}\left(x\left(t\right),y\left(t\right)\right)$ is the nodal displacement column vector that concatenates the displacement DOF of the structure’s nodes, $x\left(t\right)$, with those of the MTMD, $y\left(t\right)$. The notations $\left(\dot{}\right)$ and $\left(\ddot{}\right)$ indicate the first and second derivatives with respect to time, respectively, hence denoting the velocity and acceleration vectors for each DOF in the system. Moreover, $\tilde{\mathbf{P}}\left(t\right)=\mathrm{col}\left(P\left(t\right),0\right)$ is the column vector that defines external forces applied to the structure’s nodes. As it will be discussed further in the text, the variable of external forces $\tilde{\mathbf{P}}\left(t\right)$ is separated into two categories, the first being force functions containing non-aerodynamic excitations, and the second being the aerodynamic disturbances applied to the structure.

It is of interest to understand how the MTMD modifies the global FE matrices. The general form of the $\tilde{\mathbf{M}}$ matrix, from Equation (1), is described by Equation (2a,b) as,

$$\tilde{M}=\left[\begin{array}{cc}M& 0\\ 0& m\end{array}\right]$$

(2a)

$$m=\left[\begin{array}{cccccccccc}{m}_{11}& 0& & & & \dots & & & & 0\\ 0& \ddots & & & & & & & & \\ & & m_{1K}& & & & & & & \\ & & & m_{21}& & & & & & \\ & & & & \ddots & & & & & ⋮\\ ⋮& & & & & m_{2K}& & & & \\ & & & & & & \ddots & & & \\ & & & & & & & m_{N1}& & \\ & & & & & & & & \ddots & 0\\ 0& & & & \dots & & & & 0& m_{NK}\end{array}\right]$$

(2b)

where $\tilde{\mathbf{M}}$ is the block diagonal matrix assembled with $M$, which is the N × N material mass matrix of the global system, and $m$, represents the N.K × N.K mass inclusion of the MTMD considered.

The modified global stiffness matrix, $\tilde{\mathbf{K}}$, can be written in Equation (3a–d) as,

$$\tilde{\mathbf{K}}=\left[\begin{array}{cc}K+\overline{k}& \stackrel{=}{k}\\ {\stackrel{=}{k}}^T& k\end{array}\right]$$

(3a)

$$\overline{k}=\left[\begin{array}{ccc}{k}_{11}+k_{12}\cdots +k_{1K}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & k_{N1}+k_{N2}\cdots +k_{NK}\end{array}\right]$$

(3b)

$$\stackrel{=}{k}=\left[\begin{array}{cccccccccccc}-k_{11}& -k_{12}& \cdots & -k_{1K}& 0& 0& 0& 0& & \cdots & & 0\\ 0& 0& 0& 0& -k_{21}& -k_{22}& \cdots & -k_{2K}& & \cdots & & 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \ddots & \cdots & & 0\\ ⋮& & & & & & & & & & & ⋮\\ 0& & & & \cdots & & & 0& -k_{N1}& -k_{N2}& \cdots & -k_{NK}\end{array}\right]$$

(3c)

$$k=\left[\begin{array}{cccccccccc}{k}_{11}& 0& & & & \dots & & & & 0\\ 0& \ddots & & & & & & & & \\ & & k_{1K}& & & & & & & \\ & & & k_{21}& & & & & & \\ & & & & \ddots & & & & & ⋮\\ ⋮& & & & & k_{2K}& & & & \\ & & & & & & \ddots & & & \\ & & & & & & & k_{N1}& & \\ & & & & & & & & \ddots & 0\\ 0& & & & \dots & & & & 0& k_{NK}\end{array}\right]$$

(3d)

where $K$ is the global material stiffness matrix of the structure, which is assembled with the stiffness compliance matrices of the MTMD inclusion, defined as $\overline{k}$, $\stackrel{=}{k}$ and $k$. The individual values of stiffness for each TMD are defined as $k_{ij}$, where $i=1, \cdots , N$ and $j=1, \cdots , K$. The notation $\left({}^T\right)$ indicates the transposition operation of a matrix.

Analogous to the stiffness matrix-assembling process, the modified global damping matrix $\tilde{\mathbf{C}}$ can be formulated in Equation (4a–d) as,

$$\tilde{\mathbf{C}}=\left[\begin{array}{cc}C+\overline{c}& \stackrel{=}{c}\\ {\stackrel{=}{c}}^T& c\end{array}\right]$$

(4a)

$$\overline{c}=\left[\begin{array}{ccc}{c}_{11}+c_{12}\cdots +c_{1K}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & c_{N1}+c_{N2}\cdots +c_{NK}\end{array}\right]$$

(4b)

$$\stackrel{=}{c}=\left[\begin{array}{cccccccccccc}-c_{11}& -c_{12}& \cdots & -c_{1K}& 0& 0& 0& 0& & \cdots & & 0\\ 0& 0& 0& 0& -c_{21}& -c_{22}& \cdots & -c_{2K}& & \cdots & & 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \ddots & \cdots & & 0\\ ⋮& & & & & & & & & & & ⋮\\ 0& & & & \cdots & & & 0& -c_{N1}& -c_{N2}& \cdots & -c_{NK}\end{array}\right]$$

(4c)

$$c=\left[\begin{array}{cccccccccc}{c}_{11}& 0& & & & \dots & & & & 0\\ 0& \ddots & & & & & & & & \\ & & c_{1K}& & & & & & & \\ & & & c_{21}& & & & & & \\ & & & & \ddots & & & & & ⋮\\ ⋮& & & & & c_{2K}& & & & \\ & & & & & & \ddots & & & \\ & & & & & & & c_{N1}& & \\ & & & & & & & & \ddots & 0\\ 0& & & & \dots & & & & 0& c_{NK}\end{array}\right]$$

(4d)

where $C$ is the global material damping matrix of the structure, which is assembled with the damping compliance matrices of the MTMD inclusion, defined as $\overline{c}$, $\stackrel{=}{c}$, and $c$. The individual values of damping for each TMD are defined as $c_{ij}$, where $i=1, \cdots , N$ and $j=1, \cdots , K$.

From Equation (1), the response $u\left(t\right)$ can be expressed as a linear combination of harmonic motion functions, $q_i\left(t\right)$, with distinct deformation patterns, named mode shapes ${\mathsf{\phi }}_i$, that are associated with discrete values of frequency, also known as natural frequencies ${\omega }_{n_i}$ [40]. The mode shapes and natural frequencies of the structure with MTMD can be found by means of the formulation of an eigenvalue problem considering the undamped free vibrating system. These relations are given by Equation (5a,b),

$$\left(\tilde{\mathbf{K}}-{\omega }_{n_i}\tilde{\mathbf{M}}\right){\mathsf{\phi }}_i=0$$

(5a)

$$u\left(t\right)={\displaystyle \sum }_{i=1}^R{\mathsf{\phi }}_i{q}_i\left(t\right), R=N+K$$

(5b)

FE models of aerospace structures usually contain thousands of DOF, which often causes a high computational cost to solve the numerical problem. In order to simplify the equation of motion for faster computation of solutions, given by Equation (1), the property of orthogonality between the mode shapes can be exploited in order to partially uncouple the DOF. The projection of the mass and stiffness matrices into the modal space results in diagonal matrices. While an undamped or classically damped system (e.g., a proportionally damped multi-degree-of-freedom system) can be entirely decoupled, the MTMD added to the system introduces viscous damping, creating a non-classically damped system with an un-diagonalized symmetric modal-damping matrix. The resulting mass, damping, and stiffness matrices defined in the modal, or generalized, coordinates are presented in Equation (6) as,

$${\tilde{\mathbf{M}}}_q={\mathsf{\phi }}_i{}^T\tilde{\mathbf{M}}{\mathsf{\phi }}_j=\{\begin{array}{c}0, i\ne j\\ m_{q_j}, i=j\end{array}$$

(6a)

$${\tilde{\mathbf{C}}}_q={\mathsf{\phi }}_i{}^T\tilde{\mathbf{C}}{\mathsf{\phi }}_j$$

(6b)

$${\tilde{\mathbf{K}}}_q={\mathsf{\phi }}_i{}^T\tilde{\mathbf{K}}{\mathsf{\phi }}_j=\{\begin{array}{c}0, i\ne j\\ k_{q_j}, i=j\end{array}$$

(6c)

2.2. Aerodynamic Model

When investigating the dynamic behavior of aerospace structures, it is pivotal to consider the appropriate aerodynamic model to determine the external excitations on the structure. In the literature, there are a plethora of methods that allow one to achieve excellent approximations for the computation of the distribution of aerodynamic forces, such as the strip theory [41], the vortex-lattice method [42], or the doublet-lattice method (DLM) [43]. For this work, the DLM is utilized to compute the unsteady aerodynamic influence coefficients (AIC), defined in the frequency domain.

Atmospheric Gust Disturbances

Aerospace structures are constantly subjected to induced atmospheric disturbances and turbulence while in operation. Gusts are local disturbances that affect the pressure distribution around aerodynamic surfaces, hence generating additional loads that can potentially result in instability or excessive vibrations. Vertical or lateral gusts are created when a gust-induced velocity gradient is normal to the flow path, while head-on gusts are defined when the gust-induced velocity gradient is parallel to the flow path. The force generated by the gust is found from the turbulence wind-gust velocity, $\Delta \mathit{v}$, that is obtained by means of the formulation of a chosen gust model. The total wind velocity, $v$, is given by a combination of the mean velocity, ${\mathit{v}}_m$, and the turbulence wind-gust velocity, $\Delta \mathit{v}$, as in Equation (7).

$$\mathit{v}={\mathit{v}}_m+\Delta \mathit{v}$$

(7)

One type of disturbance profile is the discrete gust [44,45], in which an isolated vertical disturbance creates a pulse-type profile in the excitation, such as the one-minus cosine model. The expression that models the spatial behavior of the one-minus cosine gust velocity, $\Delta \mathit{v}$, is given by Equation (8a,b),

$$\Delta \mathit{v}\left({\mathit{x}}_g\right)=\frac{v_{g0}}{2}\left(1-\cos \frac{2\pi {\mathit{x}}_g}{L_g}\right), 0\le x_g\le L_g$$

(8a)

$${\mathit{x}}_g={\mathit{v}}_{m}t$$

(8b)

where $v_{g0}$ is the peak gust disturbance or design gust velocity, $L_g$ is the gust length, ${\mathit{x}}_g$ is the position of the structure in the spatial coordinates of the gust with respect to a fixed origin, and t is the respective measured time from the position ${\mathit{x}}_g=0$. By combining both Equations (8a,b), one obtains the aerospace structure response in the time domain.

In order to model continuous turbulence within the atmospheric boundary layer, a good representative model is given by the DS due to the fact that it takes into account the ground effects [46,47]. In this model, the wind gust is represented by a uniformly spread force acting on the lifting surface from a certain direction. The disturbance gust component is characterized by a random process defined by the DS formulation, which depends on the mean wind speed (${\mathit{v}}_m$) and terrain roughness parameters. It can be expressed by Equation (9a,b),

$$S_D\left(\omega \right)=4800{\mathit{v}}_m\kappa \frac{x}{{\left(1+x^2\right)}^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$$

(9a)

$$\kappa =\frac{1}{{\left[2.5\ln \left(\frac{z}{z_0}\right)\right]}^2}$$

(9b)

where $x=1200 \omega /{\mathit{v}}_m$ ($\omega /{\mathit{v}}_m$ is the wave number or inverse wavelength in cycles per meter), $\kappa$ is the surface drag coefficient, $z$ is the height from the ground up to the center of the structure, and $z_0$ is height of the terrain roughness, which is measured between $z_0=3.94$ to $11.81$ in at Goldstone, California [48]. The wind-gust velocity, $\Delta \mathit{v}$, is found by applying a white-noise input to a fourth-order filter that approximates the DS model [49,50].

2.3. Time-Consistent Loads and Displacements

The methodology of a Time-Consistent Load and Displacement (TCLD) is employed to compute the structural pointing error due to aerodynamic disturbances that are computed in the frequency domain, namely, PSD gust disturbances. Here, as multiple degrees of freedom are employed to define the pointing error (discussed in further detail in Section 2.4), it is very important to use the response for those degrees of freedom at the same instant, i.e., the time-consistent response. Here, while the worst response for a specific degree of freedom can be recovered in the frequency domain at a certain frequency, the response of other degrees of freedom, while happening at the same frequency, may be occurring at a different time instant.

To achieve time-consistent loads and displacements in the frequency domain, two parameters must be defined. The first is the critical frequency, ${\omega }_{crit}$, which defines the circular frequency of a certain DOF, measured in degrees per second. The second is the critical phase angle, ${\phi }_{crit}$, which is the correspondent phase angle of the chosen critical frequency at the same DOF, measured in degrees. The critical frequency is usually selected according to the maximum amplitude of the frequency-response function, given a selected frequency range, of a specific DOF considered in the model. With the values of both parameters determined, the critical time with respect to the reference DOF, $t_{crit}$, is then defined as in Equation (10).

$$t_{crit}=\frac{{\phi }_{crit}}{{\omega }_{crit}}$$

(10)

The parameter $t_{crit}$ represents the reference instance that will allow all DOFs to be time consistent with each other. Then, if we consider an arbitrary DOF quantity to become time consistent and name it the target station, the conversion can be expressed by Equation (11),

$$U_{stn,TCLD}=U_{stn}\cos \left({\phi }_{crit}-{\phi }_{stn}\right)$$

(11)

where $U_{stn,TCLD}$ is the time-consistent load or displacement for a specific station, U_stn is the out-of-phase load or displacement of the same station, and ${\phi }_{stn}$ is the phase angle of the load or displacement determined from the frequency-response function. The variables $U_{stn,TCLD}$ can now be operated together for different DOFs due to the application of the time-consistent formulation. In the case of displacement responses, for example, the application of the TCLD formulation could be of interest to determine a net vector that represents the magnitude of the structure’s deformation when excited with the critical frequency. As it will be presented in the next section, this concept will be part of the proposed optimization procedure.

2.4. Multi-Objective Design Optimization Procedure

In this section, it is proposed that the multi-objective optimization problem that determines the best values for the MTMD attached to the structural nodes of a given aerospace structure intends to concurrently minimize the aeroelastic response and the total mass inclusion due to the placement of the control devices. The optimization problem is formulated as follows:

Fitness function:	Minimize ${\displaystyle \sum }_{i=1}^{n_{TMD}}{m}_{TMD}\left(i\right)$
	Minimize $U_{fitness}$
Design variables:	$\mathsf{\xi }=\left\{m_{TMD},c_{TMD},k_{TMD}\right\}$
Constraints:	$0\le {\displaystyle {\displaystyle \sum }_{i=1}^{n_{TMD}}}{\mathrm{m}}_{TMD}\left(i\right)\le \left(0.25\times M\right)$
	$c_{TMD}^{min}\le {\mathrm{c}}_{TMD}\left(i\right)\le c_{TMD}^{max}$
	$k_{TMD}^{min}\le {\mathrm{k}}_{TMD}\left(i\right)\le k_{TMD}^{max}$

where $\mathsf{\xi }=\left\{m_{TMD},c_{TMD},k_{TMD}\right\}$ is the vector that contains values for mass, damping, and stiffness of each TMD considered, $n_{TMD}$ is the total number of TMD that were attached to structural grid points, and $U_{fitness}$ is a generic net output response of the system, computed with the root sum of squares (RSSQ) of the DOF. The position of each TMD could be similarly implemented using the proposed optimization framework; however, it was considered to be out-of-scope for this paper. In the case of a frequency-domain analysis, the output response of the DOF is calculated using the TCLD formulation for a determined critical frequency on the frequency-response function. $M$ is the total mass of the structure analyzed. It is important to mention that, in order to simplify the variability of the masses of the TMD, all were considered the same, and the total inclusion should be less than or equal to a chosen value of 25% of the uncontrolled structure’s mass. The variables $c_{TMD}^{min}$, $c_{TMD}^{max}$, $k_{TMD}^{min}$, and $k_{TMD}^{max}$ represent the minimum and maximum bounds of the design space of the damping and stiffness parameters from the MTMD set analyzed.

The optimization procedure is carried out by means of an elitist MOGA [51,52,53] that computes the MTMD optimal parameters for the structural model, which is subjected to the atmospheric wind disturbances as modeled in the previous sections. In general, the MOGA begins by creating an initial population of a predefined size that satisfies the linear constraints and bounds imposed on the optimization problem. Then, the “parents individuals” from the current population are selected by means of a binary tournament, and the children are created by crossover and mutation from the selection tournament. The “children individuals” are evaluated using the fitness-function outputs and feasibility, creating a score for each individual. The obtained “children individuals” are then combined, in matrix form, with the current population evaluated, generating the extended population. Within this new set, the rank and crowding distance are calculated for all individuals. The rank is essentially a parameter that measures the chance of selection and is related to the dominance of other individuals. When an individual is not dominated by others, i.e., is classified as independent, it returns lower rank values, which represent the highest chances of feasibility and selection. The crowding distance is a measurement of the closeness of individuals of the same rank to their nearest surrounding neighbors. Higher crowding distances among individuals of the same rank have a higher chance of selection. Finally, the last step of the MOGA is to trim the extended population generated to the actual predefined population size by keeping a certain number of individuals with similar rank values.

Regarding the MSC NASTRAN^® modelling of each TMD attached to a target structure, it can be performed in a variety of ways. The modelling performed in this paper consists of concurrently using three different elements: the CELAS2, which defines a scalar spring element; the CDAMP2, which defines a scalar viscous damp element; and the CMASS2, which defines a scalar concentrated mass element [54]. The referred elements support multiple DOFs of actuation, resulting in a multiple DOF TMD. In order to attach the control devices to a structural grid node, it is necessary to create an extra dummy node that is used to place the CMASS2 element and interconnect the CELAS2 and CDAMP2 elements, as depicted in Figure 2. With this modelling it is possible to insert multiple TMD in a structure and evaluate the effect of changing the mass, stiffness, and damping of each element that potentially minimizes the amplitude of response under dynamic excitations by manipulating the equation of motion given by Equation (1).

Figure 3 depicts the steps of the optimization procedure implemented. Initially, the inputs are defined for the simulation setup file for the modelling of the dynamic aeroelastic analysis in MSC NASTRAN^® (SOL146), which encompasses the geometry, element connection, and properties of the structure, along with the aerodynamic DLM panels and gust excitation inputs. The latter is computed in a separate script for different initial conditions and atmospheric conditions. The MTMD structural grid points of placement, the optimization bounds ($m_{TMD},c_{TMD},k_{TMD}$ values) and the GA setup parameters are also defined. In sequence, the iterative design optimization procedure begins by evaluating the fitness function considering the current individual from the defined population. The total mass inclusion of the MTMD is first calculated from the design variable inputs, followed by $U_{fitness}$. To determine the latter, a sequence of sub-steps is necessary. First, two files for input to MSC NASTRAN^® are generated: the extra dummy grid points created for the placement of the MTMD according to the predefined location and $n_{TMD}$ used in the analysis, and the actual elements that compose the MTMD given the values from the design variables. Those are represented by the CELAS2, CDAMP2, and CMASS2, which are the spring, viscous damper, and concentrated mass elements, respectively. Second, MSC NASTRAN^® is started, and the dynamic aeroelastic analysis is computed with the given input files previously described. Third, the output response file is generated in the PCH format and loaded in MATLAB^® and is defined in terms of displacement of the DOF from a selected node. Fourth, if the frequency-domain analysis is being performed, the displacement output response is adjusted with the TCLD formulation, where the peak response is determined for a critical frequency, and then $U_{fitness}$ is calculated by means of a RSSQ of the DOF. If the time-domain analysis is being performed, $U_{RSSQ}$ is calculated directly with the RSSQ of the displacement DOF considered and integrated over the time length of the signal analyzed, such that the total area under the $U_{RSSQ}\times t$ curve is minimized, therefore defining the value of $U_{fitness}$ for this case. The optimization procedure continues the iterative process until a stopping criterion is reached, such as the chosen maximum number of generations or the average change in the spread of the Pareto front.

The multi-objective optimization method determines several sets of feasible optimal solutions according to the constraints defined; however, only a single solution can be implemented at a time for a structure. Therefore, considering all optimal solutions, the best of them can be selected in a variety of ways. In this work, to quantitatively determine the best design variables set among the Pareto front points, the utopian point method was used [55]. The utopian point can be found by means of the coordinates of the extremities of the Pareto front, thus defined by values that simultaneously minimize both objective functions. The distance from the utopian point to the Pareto front can be measured by means of the weighted Euclidean metric method, and the minimum distance overall defines the optimal set of design variables, representing the best trade-off between the objectives of the problem. The described method is summarized in Equation (12),

$$d\left(\mathsf{\xi }\right)=\underset{\mathsf{\xi }}{\min }{\left\{{\displaystyle \sum }_{i=1}^N{\left|w_i\left(f_i-f_i^0\right)\right|}^2\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}, i=1,\cdots ,n$$

(12)

where $d=d\left(\mathsf{\xi }\right)$ is the calculated minimum utopian metric between the utopian point, $F^0=\left(f_1^0,f_2^0,f_3^0,\dots ,f_N^0\right)=f_i^0$, and the Pareto front point, $F\left(\mathsf{\xi }\right)=\left(f_1,f_2,f_3,\cdots ,f_N\right)=f_i$, with a total of $N$ objectives. The weights $w_i, i\in \left\{1,\cdots ,n\right\}$ are defined as individual contributions originating from the objectives, representing a criterion of importance defined by the user [55]. For simplicity of solution, we used the weights set as unitary values in this work, such that the trade-off between objectives is equally balanced.

3. Case Study

The methodology presented in the previous sections will now be applied to a truss-like antenna structure that is subjected to atmospheric wind disturbances and turbulence. The pointing accuracy of the antennas is a parameter of pivotal importance for obtaining high-resolution signals. Due to the large dimensions of those aerospace structures, they are constantly subjected to atmospheric wind disturbances; therefore, minimizing the aeroelastic response is of great interest to allow the minimization of the pointing error. The implementation of MTMD devices is an interesting approach for this application since it is highly efficient in attenuating dynamic modal responses and is also a passive vibration control technique. The performance of this approach has not yet been investigated in the literature and will therefore be presented in this paper.

3.1. FE Model of a Truss-Like Antenna Structure

The case study structure is essentially composed of a Kagome-like truss inner geometry with a hexagonal outer shape and four boom arms that support the secondary reflector with an offset to the primary reflector, as depicted in Figure 4. Each member of the truss has 25 in of length with a cross-section area of 1 in², manufactured with aluminum 6061-T6 [56] (modulus of elasticity of 10,000 ksi and density of 0.0975 $\mathrm{lb}/{\mathrm{in}}^3$), resulting in a total estimated mass of approximately 1183 lb. The structure is modeled with a structural damping coefficient of 0.03, a value commonly used in the aerospace industry for continuous metal structures [57]. The choice of constitutive element to represent the antenna structure in MSC NASTRAN^® is a beam element (CBAR) with two nodes and six DOFs per node. The frame is supported by 12 nodes located in the backside truss of the primary reflector, as indicated in Figure 4.

In order to evaluate the pointing accuracy of the antenna, two nodes of importance need to be defined: the reference and target nodes. The first is located at the center of the frontside truss, while the second is located at the junction of the four boom arms, also shown in Figure 4. The pointing accuracy is modelled as the deviation angle between the target and reference nodes using Equation (13a–d),

$$\Theta ={\tan }^{-1}\frac{\Delta x}{\Delta L}$$

(13a)

$$\Phi ={\tan }^{-1}\frac{\Delta y}{\Delta L}$$

(13b)

$$\Psi ={\tan }^{-1}\frac{\Delta z}{\Delta L}$$

(13c)

$$U_{RSSQ}=\sqrt{{\Theta }^2+{\Phi }^2+{\Psi }^2}$$

(13d)

where $\Theta$, $\Phi ,$ and $\Psi$ are the calculated deviation angles, respectively, in the x, y, and z axes. The variables $\Delta x$, $\Delta y$, and $\Delta z$ are the time-consistent displacements originated from the dynamic aeroelastic response in the x, y, and z directions. $\Delta L$ is the fixed offset distance between the target and reference nodes, and $U_{RSSQ}$ is the net deviation angle of the structure, thus defining the pointing accuracy.

3.2. Aeroelastic Response of the Uncontrolled Model

The primary reflector of the antenna frame is considered the main aerodynamic surface that will result in relevant values of disturbance caused by the wind-gust excitations, as modelled with the DLM panels. It is approximated as a flat surface due to the low thickness to chordwise length ratio, evaluated at approximately 6%. The aerodynamic approximation is implemented in MSC NASTRAN^® with two main aero panels containing 100 boxes each, where the wind-gust disturbance is defined as a vector perpendicular to the plane of the primary reflector surface and the mean wind speed vector being parallel to it, as depicted in Figure 5. Introducing MTMD into the system does not alter the force-aerodynamic influence coefficient matrix, as the defined splining matrix used to transfer forces between the aerodynamic and structural sets was consistent between uncontrolled and optimized models.

The DLM is commonly used as a reference tool in the aerospace industry in order to determine the response of structures under unsteady aerodynamics. Its main advantage relies on the fact that it has a the low computational cost and is integrated with commercial software such as MSC NASTRAN^®. However, this method presents limitations in terms of its approximations due to the use of linearized potential equations, therefore not describing essential effects for aerodynamic analysis, such as geometry thickness, fluid viscosity, or shockwave formation, often resulting in overprediction of lift in aerodynamic surfaces [58,59,60]. Considering the DLM limitations in estimating aerodynamic coefficients, it is a common procedure in the industry to match the steady results at zero frequency with a higher-order model by using correction factors, such as those computed with CFD analysis or validated with steady-state wind tunnel data, in order to guarantee a good compromise between performance and accuracy in the results [61].

Two operational cases were considered for the analysis of the proposed antenna model, in which the range of actuation was calculated for TDG and PSD formulations, respectively, using Equations (8a,b) and (9a,b). The gust scale length utilized in the analysis was based on the values given by Nieto et al. [1], defined as $L_g=47,244.094$ in. In addition, the mean wind speed and the disturbance speed of the scenarios utilized are summarized in

Table 1. Operational case summary [62].

Case Identification	Disturbance Speed $\Delta \mathit{v}$ (in/s)	Mean Wind Speed ${\mathit{v}}_{\mathit{m}}$ (in/s)	Mach Number	Dynamic Pressure (lb/in3)
Primary	109.10	346.95	0.02615	0.005686
Secondary	339.20	1093.61	0.08173	0.055524

and were given by the industrial partner of the project and the National Aeronautics and Space Administration (NASA) based on real operational conditions estimated at installation sites [62]. The wind-gust models considering the referred operational scenarios are depicted in Figure 6.

Figure 7 presents the frequency and time-domain aeroelastic deviation angle responses for both primary and secondary operational cases of the uncontrolled structure, measured at the target node, considering the PSD and TDG models, respectively. For the TDG, after 10 s of excitation, the magnitude of the response was negligible and therefore truncated. It can be observed that the pointing error in the direction of $\Theta$, even though close to the $\Phi$ response, is the greatest overall for both gust models and operational scenarios considered. This behavior is expected due to the fact that the gust disturbance input was essentially given for the x direction, thus resulting in a response more significant in the referred direction when compared to the others.

Table 2. Evaluation of the fitness function for the uncontrolled antenna structure.

Case ID	DS ${\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ (deg)	TDG ${\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ (deg × s)
Primary	5.58 × 10−5	2.76 × 10−4
Secondary	3.78 × 10−4	8.39 × 10−3

summarizes the evaluation of the fitness function of the optimization algorithm for the uncontrolled structure in order to be taken as reference values for the results section.

The frequency response output from Figure 7a,b presents an accentuated peak, estimated at about 2.8 Hz, which is coincident with the first mode shape of the antenna structure, as indicated in Figure 8. The other frequencies in the spectrum showed magnitudes considerably lower than the one related to the first mode shape, demonstrating their insignificant participation in the overall response. Therefore, the optimization algorithm should minimize the amplitude of response near that frequency by means of the inclusion of MTMD with optimum parameters. In addition, from Figure 7c,d, the time domain response is determined by a summation of a truncated finite number of mode shapes of the structure, and, in this case, it mainly represents the response of the first vibration mode due to its greater dominance. Hence, the focus now is to study the passive control of the aeroelastic response by including MTMD devices symmetrically within the antenna structure, investigating optimal design parameters, and comparing different placements by using the optimization framework developed in the previous sections.

3.3. Aeroelastic Response Attenuation Using Optimal MTMD

For the investigation of pointing error minimization of the antenna structure with MTMD control devices, three different configurations of auxiliary masses are proposed. Configuration 1 places the MTMD in the structural grid points of the larger hexagonal shape of the backside truss, totaling an inclusion of 61 devices. Configuration 2 places the MTMD in the structural grid points of the boom arms, totaling an inclusion of 29 devices. Finally, Configuration 3 merges both of the previous configurations, therefore totaling an inclusion of 90 devices. As it was previously mentioned, all MTMD devices were defined with the same mass, which varies from 0% to a maximum of 25% inclusion with respect to the total antenna structure’s mass.

Table 3. Design variable bounds for each configuration investigated.

Configuration ID	${\mathit{n}}_{\mathit{T}\mathit{M}\mathit{D}}$	${\mathit{m}}_{\mathit{T}\mathit{M}\mathit{D}}^{\mathit{m}\mathit{a}\mathit{x}}$ (lb)	${\mathit{c}}_{\mathit{T}\mathit{M}\mathit{D}}^{\mathit{m}\mathit{i}\mathit{n}}$ (lbf.s/in)	${\mathit{c}}_{\mathit{T}\mathit{M}\mathit{D}}^{\mathit{m}\mathit{a}\mathit{x}}$ (lbf.s/in)	${\mathit{k}}_{\mathit{T}\mathit{M}\mathit{D}}^{\mathit{m}\mathit{i}\mathit{n}}$ (lbf/in)	${\mathit{k}}_{\mathit{T}\mathit{M}\mathit{D}}^{\mathit{m}\mathit{a}\mathit{x}}$ (lbf/in)
1	61	5.0	0.29	285.50	2.86	2855
2	29	10.0	0.29	285.50	2.86	2855
3	90	3.3	0.29	285.50	2.86	2855

summarizes the upper and lower boundaries chosen for each individual TMD parameter. A complete description of the MTMD device nomenclature with respect to the structural grid nodes is presented in the Appendix of the paper.

Due to the considerable number of design variables, the parameters of the MOGA must be wisely chosen in order to reduce the convergence time. Details about the meaning of the parameters can be found in Kalyanmoy [38]. The following parameters were used in the simulations:

• Maximum number of generations: 30;
• Population size: 300;
• Crossover fraction: 80%;
• Elitism probability: 2%;
• Migration factor: 20%;
• Migration interval: 20;
• Pareto fraction: 35%.

With those parameters, a total of 9000 analyses are carried out within the MOGA optimization procedure, per configuration. The results for the optimization framework are now presented, considering the operational cases and both models of atmospheric disturbances. The results are given in terms of Pareto fronts, output response functions for the target node, and charts that indicate the optimum values of the MTMD found.

The optimization procedure presented in Section 2.4 was implemented for the gust excitation modelled with the PSD. Representatively, an example of output from the framework is shown in Figure 9, where the pointing error and mass inclusion are minimized for the primary operational case. In Figure 9a, it is possible to observe that the decrease in the pointing error is associated with a direct increase in the mass inclusion in the antenna, therefore establishing a trade-off problem. The utopian point indicated encompasses the extremes of the Pareto front, hence considering all individuals obtained at the convergence of the optimization algorithm. According to the demands of the project, it is possible to select different best individuals by means of the weight function previously described. When they are taken as unitary values, both objectives have the same relevance in defining the best individual, which is the case in this analysis. Hence, the best individual found for the referred operational case and configuration selected allows a reduction in the pointing error to 62%, while requiring a total mass inclusion of 5.9%. In Figure 9b, all individuals from the Pareto front were plotted, and it is observable that the maximum mass inclusion provided a significant attenuation of the first resonance frequency, reaching a reduction in the pointing error of 87.6% with a total mass inclusion of 23.9%. Another relevant observation is that all angles had their response attenuated by the inclusion of the MTMD, therefore demonstrating the capability of MDOF aeroelastic response mitigation of those devices. Such a fact was also observed in all other operational cases and configurations investigated.

In order to compare the Pareto fronts obtained with different configurations under both operational scenarios, Figure 10 is presented. It is possible to conclude that Configuration 1 presented the overall best performance for both operational cases, providing the best trade-off between mass inclusion and pointing error minimization. Coincidently, such an observation was also valid for the TDG simulations, as it will be shown in the next session. A summary of the Pareto front values for the best individuals in Figure 10 is presented in

Table 4. Pareto front values for the best individuals under the PSD gust excitation model.

Configuration ID	Primary Case		Secondary Case
	% Mass Inclusion (lb)	$\% {\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ Reduction (Deg)	% Mass Inclusion (lb)	$\% {\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ Reduction (Deg)
1	5.9	62.0	4.6	51.6
2	2.9	43.0	3.2	50.1
3	4.9	60.0	4.0	39.5

.

Figure 11 depicts the output response for the dominant DOF angle, $\Theta$, for all best individuals of the different configurations under the two operational scenarios, compared with the uncontrolled response of the antenna, hence demonstrating the better attenuation capability of Configuration 1. It is known that the placement of MTMD is more effective in regions of higher displacement, therefore providing better attenuation capabilities [21]. The placement of the MTMD in Configuration 1 and Configuration 3 presented advantages with respect to Configuration 2 due to the fact that they are mainly located in grid points of maximum displacement of the first resonance mode shape from the antenna structure, according to Figure 8, therefore resulting in superior attenuating capabilities. Additionally, the amplitude of the frequency response is proportional to Pareto front values for pointing error reduction. The difference between the amplitude of the frequency response for Configuration 2 and Configuration 3 within Figure 11a,b is primarily driven by the reduction in pointing error previously displayed in Table 4. Since a trade-off is made between the two design objectives, the best design for Configuration 2 has a larger point-error reduction when comparing the primary operational case to the secondary case. Meanwhile, Configuration 3 has a decrease in pointing error reduction. The optimized configurations for the primary and secondary operational cases cannot be directly compared as they represent different design selections from optimization algorithm convergence. The best design was selected to minimize pointing error and mass inclusion with equal relevance, where the optimization results are dependent on the input gust model.

Now, the two operational cases are applied to the antenna structure considering the TDG model, and the performances of each MTMD configuration are compared. Similarly, as in the previous gust model analysis, a representative output from the framework is presented, where the concurrent minimization of mass inclusion and area underneath the time response function, for Configuration 1, considering the secondary operational case, are depicted in Figure 12a. The choice of the second objective is based on the fact that a smaller area under the curve would provide faster attenuation to the gust signal, which presents a behavior, such as an impulse signal due to its short duration. Figure 12b displays the deviation angle in terms of the RSSQ, therefore providing the net displacement of the DOF angles. As it is presented, an alike behavior of direct proportionality between the mass inclusion and the reduction in the U_fitness was characterized, as in the case of the PSD gust model, achieving an overall maximum of 58.5% U_fitness minimization with 23.6% added mass in the antenna structure.

By applying the utopian point method to the optimal individuals generated from the MOGA framework, one can determine the best sets from the Pareto frontiers, considering all configurations and operational cases, as presented in Figure 13. Observing the graphs depicted, it is possible to conclude that while maintaining a relatively similar mass inclusion percentage, Configuration 1 presented the highest reduction in $U_{fitness}$ for both configurations, and operational cases investigated. The objectives are summarized in

Table 5. Pareto front values for the best individuals under the TDG gust excitation model.

Configuration ID	Primary Case		Secondary Case
	% Mass Inclusion (lb)	$\% {\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ Reduction (Deg × s)	% Mass Inclusion (lb)	$\% {\mathit{U}}_{\mathit{f}\mathit{i}\mathit{t}\mathit{n}\mathit{e}\mathit{s}\mathit{s}}$ Reduction (Deg × s)
1	3.6	39.2	4.7	36.8
2	3.0	29.4	4.4	33.0
3	3.7	24.2	6.4	22.8

. Moreover, the performance of the configurations is displayed in Figure 14, where the time response for the TDG excitations is compared. It can be seen that the inclusion of the MTMD devices provides a considerable attenuation of the uncontrolled signal between the time intervals presented, where Configuration 1 resulted in the best deviation angle reduction performance.

From the analysis performed in this section, it was found that the overall best performance was achieved with the MTMD in Configuration 1, therefore the optimal values for the mass, damping, and stiffness are given, respectively, in

Table 6. Summary of the best values for the MTMD masses for Configuration 1 under the operational scenarios modelled with the PSD and TDG.

Operational Case	Optimal Values for ${\mathbf{m}}_{\mathit{T}\mathit{M}\mathit{D}}$ (lb)
DS—Primary	1.14
TDG—Primary	0.69
DS—Secondary	0.90
TDG—Secondary	0.90

and Figure 15. As it is expected, the multitude of boundary conditions explored returned non-unique optimal values for the MTMD among the cases analyzed. This is an inherent characteristic of this type of passive control device, such that the optimal parameters strongly depend on the wind disturbance velocity for different operational conditions [63].

4. Conclusions

In this paper, it was proposed that an optimization framework to concurrently attenuate the dynamic aeroelastic response and minimize mass inclusion, applied to aerospace structures, by means of the use of MTMD, which belongs to the category of passive control devices. We introduced an algorithm that integrates routines from MATLAB^®/MSC NASTRAN^®, using the dynamic aeroelastic solution 146 and an elitist MOGA to determine the optimal parameters of the control devices to minimize the objective functions, considering different models of atmospheric disturbance excitations in both the frequency and time domains. One of the main characteristics of the developed framework is that it can be essentially implemented for any structure or boundary condition modelled within the MSC NASTRAN^® environment, therefore being extremely versatile and applicable to other fields of engineering than aerospace structures. To investigate its performance, an FE model of a truss-like antenna was created, and its aeroelastic response was assessed in terms of pointing accuracy. Furthermore, the use of MTMD has not yet been investigated in the literature applied to this type of structure, thereby consolidating one important novelty of the work. Given the evaluation of the mode shapes in the peak frequency response, three different configurations of MTMD placement were investigated and compared. From the optimization results, the best individuals were determined according to the implementation of the Utopian point method considering the weighted Euclidean metric for different operational load cases. It was found that, for both atmospheric disturbance models, the overall best individual was given by Configuration 1, which presents the MTMD placed in the backside truss of the antenna. This was an expected result due to the fact that the placement of MTMD is more effective in regions of higher displacements, therefore providing better vibration attenuation capabilities. By utilizing the MTMD in Configuration 1, it was possible to characterize a reduction in the pointing error of the antenna by 62.0% and 51.6% for the primary and secondary operational cases considering the DS model, respectively, and 39.2% and 36.8% for the primary and secondary operational cases considering the TDG model, respectively, demonstrating the excellent performance of the implementation of such a passive control method. Therefore, this paper contributed to proposing advanced and useful tools to determine the optimal parameters of MTMD passive control devices under complex loading scenarios, as a primordial step towards the use of such systems in applications that commonly employ active or semi-active solutions, which require sophisticated management and calibration.