A Comparative Analysis of the Response-Tracking Techniques in Aerospace Dynamic Systems Using Modal Participation Factors

Abstract

Mechanical structural systems are subject to multiple dynamic disturbances during service. While several possible scenarios can be examined to determine their design loading conditions, only a relatively small set of such scenarios is considered critical. Therefore, only such particular deterministic set of critical load cases is commonly employed for the structural design and optimization. Nevertheless, during the design and optimization stages, the mass and stiffness distributions of such assemblies vary, and, in consequence, their dynamic response also varies. Thus, it is important to consider the variations in the dynamic loading conditions during the design-and-optimization cycles. This paper studies the modal participation factors at length and proposes an alternative to the current point-wise treatment of the dynamic equations of motion of flexible bodies during design optimization. First, the most relevant-to-structural-dynamics definitions available in the literature are reviewed in depth. Second, the analysis of those definitions that have the potential to be adopted as point-wise constraint equations during structural optimization is extended. Finally, a proof of concept is presented to demonstrate the usability of each definition, followed by a case study in which the potential advantages of the proposed extended analysis are shown.

1. Introduction

Complex mechanical systems often encounter dynamic disturbances that can affect their performance. These disturbances, whether random or systematic, must be considered at various stages of the design process to ensure the structure’s safety and efficiency. During structural optimization, changes to mass and stiffness distributions can alter the dynamic response of the system. Therefore, to optimize a structure effectively, designers must account for variations in the dynamic loads throughout the design and optimization process. Various dynamic optimization techniques have been developed for specific design and structural applications. For example, in seismic design, a method called the sum of modal compliances (SMC) [1] is used for high-rise structures. SMC involves solving design-independent optimization problems and periodically updating load vectors based on modal analysis.

The estimation of design loads can become both costly and time-intensive, depending on the complexity of the numerical model and the number of simulations needed for the design exploration. In such cases, it may be necessary to conduct load estimation independently of the overall optimization process. Even when dealing with a single scenario, simplifying the time-dependent constraint equations (dynamic equations of motion) becomes crucial, achieved either through point-wise constraints or by transforming them into a single equivalent functional equation [2,3]. In large-scale system analyses involving frequency sweeps, the interpolatory model order reduction technique [4] has been demonstrated to be effective. It approximates response functions, reducing the substantial computational workload in applications like vibro-acoustics and large-scale structural dynamics. Additionally, several iterative sub-structuring methods [5] have been devised to assess the influence of higher modes in the reduced eigenvalue equation. These approaches enable the calculation of eigenvalues and their derivatives with computational efficiency, maintaining numerical precision compared to global structural analysis.

In one of the above-mentioned approaches, developed to reduce computational expense, the point-wise treatment of the equations of motion in flexible bodies is thoroughly considered. The idea is to replace the time dependency by identifying the worst-case point (i.e., the time step and the magnitude) of the solution. However, one of the challenges posed by this treatment is that the time step at which the maximum dynamic response occurs during an iteration will change in the subsequent iteration [6,7,8,9]. To solve this issue, a simpler but more expensive approach is to consider a coarser time grid for the solution of the equation of motion [10,11,12,13]. Therefore, the development of a point-wise methodology that allows the fast and accurate re-analysis of the worst-case point could lead to significant computational savings and an improved overall optimization cycle. This is especially significant when transient loads play a fundamental role in the design of a mechanical assembly.

A search within the open literature shows that the concept of the “modal participation factor” is somewhat ambiguous, as several variants exist for various applications; each variant is used to identify specific metrics that are convenient to each discipline. For example, in the context of electrical engineering, the term “modal participation” refers to the product of the right and the left eigenvectors to measure the sensitivities of voltage to reactive power variations [14,15,16,17,18].

Considering various spheres of study, such as mechanical structural response, vibrational analysis, and structural dynamics, several definitions of modal participation factors have been developed. They are described below. In mechanical engineering, further examples of the usage of the term “modal participation factor” can be found. For example, Wallrapp and Wiedemann [19] used a quasi-static solution, from which they derived a modal participation factor to assess the importance of a mode within that solution. Likewise, Chung-Band and Bahng [20] used the concept of modal strain energy to define a modal participation factor and identify dominant modes. A similar study was carried out by Zhou [21,22] for laminated thin composite plates in aerospace applications involving severe dynamic loading conditions, in which modal participation factors were utilized as influential parameters in the stress modal analysis (SMA) and evaluation of the critical modes in the numerical finite element (FE) study. The reliability of modal participation factors in the field of structural health monitoring, using guided waves, has also been validated [23] via a semi-analytical finite element method, in which modal control was studied in excitation and sensing distributions of waveguides for composites of an irregular cross-section.

The use of the term “modal participation factor” is also a common topic in experimental modal analysis, in which different definitions have been proposed [24,25]. For example, Van Lagenhove and Brughmans [26] used the concept to identify the optimal location for sensors [27] to correlate experimental and numerical data. Using modal participation factors, the free response for the free vibration of a multilayered beam structure was investigated [28], and the influence of hole configuration parameters, such as the array, geometry, and filling ratio, was analyzed using numerical FE methods. These steps formed an integral part of the design-optimization scheme for perforated beam configuration and allowed the design decision-making process to be more focused and objective-oriented.

A comparative study of vibration modal analysis was also carried out [29] for various constituent materials, such as Kevlar and carbon-based composites, to characterize mechanical properties such as stiffness, storage modulus, and tensile strength for applications in the dental field.

Within the context of structural dynamics, the term usually refers to the effective modal mass participation factor, where the physical solution is recovered by the number of modes [30,31,32]. Within the same context, some studies have utilized the same definition to determine the eigenvector matrix [33,34], whereas, in other instances, researchers defined a more convenient mathematical form and made an indistinct use of the term [35,36,37,38,39,40]. However, it is to be noted that the definition of modal participation has an associated dichotomy [17] with it, especially in time-invariant linear systems, where it is observed that there is interchangeability between the measurement of participation of the state in modes [41] and that of participation of modes in states. This is solved using an averaging technique of the initial conditions [17] with a symmetric uncertainty. A similar approach was applied to nonlinear systems [41], where a similar dichotomy was expected and provided a precise definition for the modal participation factors for an autonomous stem with a smooth nonlinear nature.

In the field of civil engineering, an interesting formal definition, from a structural dynamics perspective, is the static modal participation factor (SMPF) definition [10,42,43,44,45,46,47,48,49]. According to the conventional definition, the steady-state dynamic response of simple numerical models can be recovered without directly solving the differential equation of motion. This idea was later extended to the case of a cantilevered wing modeled using finite elements [50]. A modal participation factor-based procedure has also been developed and verified for the model updating technique for the computational model [42] of beam structures. The model accuracy was improved by analyzing the error functions compared to the experimental measurements and the contribution of modes for analyzing structural response in both static and dynamic loading conditions.

In this paper, the concept of modal participation factors is investigated as a computationally inexpensive alternative to the current point-wise treatment of the dynamic equations of motion of flexible bodies. Taking into account the variety of definitions listed above, proposing a unified concept and a simple classification of the different definitions relevant to structural dynamics which can be found in the open literature seems convenient. A proof of concept is demonstrated by applying the developed concept to a simple system with 5 degrees of freedom (DoFs). Each selected definition explains the interpretation of the results obtained. We further extend the analysis of those definitions that can potentially replace the differential equations of motion of flexible bodies as point-wise constraint equations during structural optimization. Finally, a case study involving an aircraft wing structure is presented to demonstrate the capabilities of the developed analysis framework.

2. A Unified Concept

A unified concept for the definitions of modal participation and the modal contribution factors can be established. A modal participation factor quantifies the contribution of a mode to the structural response without the requirement of solving the dynamic equations of motion. On the other hand, a modal contribution factor indicates the percentage of the response contributed by the i-th mode to the n-th degree of freedom in the time-domain solution.

An essential difference between these two concepts is that to estimate the modal contribution, it is necessary to solve the equations of motion in the domains of space and time. In contrast, modal participation factors can be estimated without solving the equations of motion. Another important difference is that the modal contribution has a single mathematical definition. In contrast, various mathematical definitions exist that fall in the proposed concept of the modal participation factor. Below, the mathematical formulations of both concepts are presented.

2.1. Modal Contribution

The classical expression for the dynamic response of an undamped $n$-DoF system is given below:

$$\mathit{M}\ddot{\mathit{u}}\left(t\right)+\mathit{K}\mathit{u}\left(t\right)=\mathit{f}p\left(t\right)$$

(1)

where $\mathit{K}$ and $\mathit{M}$ are the stiffness and mass matrices of size $n\times n$; $n$ is the total number of degrees of freedom (DoFs) in the system; $\ddot{\mathit{u}}\left(t\right)$ and $\mathit{u}\left(t\right)$ are the vectors of nodal accelerations and displacements of size $n\times 1$; $\mathit{f}$ is the nodal external force vector of size $n\times 1$; and $p\left(t\right)$ is the time-domain function of the external force.

To uncouple the system of equations, the following change of variable is introduced:

$$\mathit{u}\left(t\right)=\mathit{\Phi }\mathit{q}\left(t\right)$$

(2)

where $\mathit{q}\left(t\right)$ is the vector of simple harmonic responses of size $N\times 1$, herein referred to as modal responses, and $\mathit{\Phi }$ is the matrix of eigenvectors of size $n\times N$. $N$ is the total number of mode shapes retained in the system.

Substituting Equation (2) into Equation (1) results in the undamped free-vibration equation of motion under the assumption of simple harmonic motion, which yields the characteristic equation

$$\left(\mathit{K}-\mathit{\lambda }\mathit{M}\right)\mathit{\Phi }=0$$

(3)

where $\mathit{\lambda }$ is a diagonal matrix of the square of the natural frequencies, ${\omega }_n$, of the system. In this context, the natural frequencies,${\omega }_n$, are those values that satisfy the condition of simple harmonic motion while the matrix of eigenvectors, $\mathit{\Phi }$, provides the relationship between the amplitudes of the modal responses, $\mathit{q}\left(t\right)$, for each DoF.

Equation (2) is in a compact form and can be expressed alternatively in component form as

$$u_n\left(t\right)=\sum _{i=1}^N{\phi }_{n,i}{q}_i\left(t\right)={\phi }_{n,1}{q}_1\left(t\right)+{\phi }_{n,2}{q}_2\left(t\right)+\cdots +{\phi }_{n,N}{q}_N\left(t\right)$$

(4)

where i is the mode number.

${\gamma }_{ni}$, indicating the modal contribution, can be mathematically defined, from Equation (4), as follows [51,52,53,54]:

$${\gamma }_{n,i}\left(t\right)={\phi }_{n,i}{q}_i\left(t\right)$$

(5)

2.2. Modal Participation Factors

Though many definitions of modal participation factors exist in the literature [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39], these can be generally classified based on their application into two classes, namely, static modal participation factors (SMPFs) and modal relative effectiveness factors (MREFs). The former aims to approximate the system modal response, ${\mathit{q}}_{\mathit{i}}\left(t\right)$, while the latter utilizes the eigenvector matrix, $\mathit{\Phi }$, to obtain an assessment of the system behavior without necessarily quantifying that response.

The SMPF class can potentially be used as point-wise constraint equations during design optimization. They are capable of determining the peak modal response and subsequently approximating the solution in the physical domain (as shown in Equation (5)) without the need for a numerical solution of the equation of motion.

The second class of modal participation factors has a more restricted use. However, for the sake of completeness, the most widely used mathematical definitions in this classification are also presented below.

2.2.1. Class I: Static Modal Participation Factors (SMPFs)

In the modal response, the amplitude is indicated as the product of a time-dependent term, $d_i\left(t\right)$ [44,50,53,54], and a constant term, ${\Gamma }_i$.

$$q_i\left(t\right)={\Gamma }_i{d}_i\left(t\right)$$

(6)

In compact notation, each ${\Gamma }_i$ constant is a component of the diagonal matrix $\mathit{\Gamma }$ of size $N\times N$, which is defined as

$$\mathit{\Gamma }\equiv diag\left({\mathit{\Phi }}^T\mathit{f}\right)$$

(7)

while $d_i\left(t\right)$ is the dynamic response of the i-th mode for a unitary forcing function, as given by Equation (8), and a component of the vector $\mathit{d}\left(t\right)$ of size $N\times 1$ in the compact notation.

$${\ddot{d}}_i\left(t\right)+{\omega }_{n_i}^2{d}_i\left(t\right)=p\left(t\right)$$

(8)

For structural design and optimization, the focus is on the peak magnitude of the response, denoted by ${\mathit{u}}^0$, rather than in its time history $\mathit{u}\left(t\right)$. Introducing Equation (6) in its compact form into Equation (2) yields Equation (9), where it is evident that the vector of peak nodal responses is directly proportional to the peak value of each entry of the vector $\mathit{d}\left(t\right)$, namely, ${\mathit{d}}^{\mathbf{0}}$.

$${\mathit{u}}^{\mathbf{0}}=\mathit{\Phi }\mathit{\Gamma }{\mathit{d}}^{\mathbf{0}}$$

(9)

Depending on the form of the unitary excitation $p\left(t\right)$, a specific method can be selected to find an analytical expression for each component of the vector ${\mathit{d}}^{\mathbf{0}}$, namely $d_i^{\mathbf{0}}$, without the need to solve Equation (8) numerically for each i mode.

The Steady-State Participation Factor

The peak amplitude of the modal response, $d_i^0$ for a unitary forcing function, can be expressed as the product of a static response $d_i^{st}$ and a dynamic amplification factor $v_i$ [42]:

$$d_i^0=d_i^{st}{v}_i$$

(10)

$$d_i^0=\frac{1}{\left({\omega }_{n_i}^2\left({\omega }_{n_i}^2-{\omega }^2\right)\right)}$$

(11)

To compare this expression with other discussed definitions, we define a steady-state modal participation fraction as the ratio between the absolute value of the steady-state response of the i-th mode and the sum of all the absolute peak modal responses:

$$L_i^{st}=\frac{\left|{\Gamma }_i{d}_i^0\right|}{{\sum }_{i=1}^N\left|{\Gamma }_i{d}_i^0\right|}=\frac{\left|\frac{{\Gamma }_i}{\left({\omega }_{n_i}^2\left({\omega }_{n_i}^2-{\omega }^2\right)\right)}\right|}{{\sum }_{i=1}^N\left|\frac{{\Gamma }_i}{\left({\omega }_{n_i}^2\left({\omega }_{n_i}^2-{\omega }^2\right)\right)}\right|}$$

(12)

The Transient I Participation Factor: Maximum Modal Response

Evaluating the steady participation factors does not consider the transient components of the time-domain solution. An alternative to approximating the peak nodal response is to consider the peak modal responses:

$$d_i^{max}=\mathit{max}\left(d_i\left(t\right)\right)$$

(13)

Solving Equation (13) requires a closed-form solution or an optimization method to locate the global maximum. However, classical gradient-based algorithms struggle to locate the global maxima because Equation (8) is inherently multimodal.

A transient modal participation fraction, $L_i^{tsI}$, is determined as the ratio between the absolute magnitude of the peak modal response estimated using $d_i^{max}$ and the sum of the absolute magnitude of the peak modal responses. This is analogous to the steady participation fraction discussed earlier.

$$L_i^{tsI}=\frac{\left|{\Gamma }_i{d}_i^{max}\right|}{{\sum }_{i=1}^N\left|{\Gamma }_i{d}_i^{max}\right|}$$

(14)

The Transient II Participation Factor: Time-Consistent Modal Response

It should be noted that Equation (13) makes an implicit assumption: the maximum response in the physical domain occurs when all the modes act in phase. This assumption is not always valid, and as a result, Equation (13) may tend to under- or overestimate the peak nodal response.

To address these limitations, we introduce a definition for the modal participation factor based on time consistency. If a mode’s contribution significantly outweighs that of the remaining modes, that mode effectively “dominates” the solution in the physical domain. Consequently, the response in the physical domain closely aligns with the behavior of its dominant mode.

Continuing from the previous discussion, to approximate the peak nodal response, we simply need to identify the time instant, $t_{dm}^{max}$, when the dominant mode exhibits its peak response and estimate the values of the remaining modes at that specific time. Mathematically, the peak nodal response can be expressed as

$$d_i^{tsII}=d_i\left(t_{dom}^{max}\right)$$

(15)

Similar to Equations (12) and (14), a transient modal participation fraction, $L_i^{tsII}$, is then defined as

$$L_i^{\mathrm{tsII}}=\frac{\left|{\Gamma }_i{d}_i^{tsII}\right|}{\sum _i^N\left|{\Gamma }_i{d}_i^{tsII}\right|}$$

(16)

Internal Load Participation Factor

The time-dependent elemental internal loads, ${\mathit{l}}_{\mathit{e}}\left(t\right)$, induced in the structure by the dynamic response are described by the following equation, in compact notation:

$${\mathit{l}}_{\mathit{e}}\left(t\right)={\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{e}}{\mathit{u}}_{\mathit{e}}\left(t\right)={\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{e}}{\mathit{\Phi }}_{\mathit{e}}\mathit{\Gamma }\mathit{d}\left(t\right)$$

(17)

where ${\mathit{K}}_{\mathit{e}}$ and ${\mathit{u}}_{\mathit{e}}\left(t\right)$ are the element stiffness matrix and the vector of nodal displacements associated with the element in the global coordinate system, respectively, and ${\mathit{C}}^{\mathit{T}}$ is the transformation matrix used to translate the internal forces to the element coordinate system.

If the term ${\mathit{u}}_{\mathit{e}}\left(t\right)$ is further expanded into its modal components, it becomes evident that four of its components ${\mathit{C}}^{\mathit{T}}$, ${\mathit{K}}_{\mathit{e}}$, ${\mathit{\Phi }}_{\mathit{e}}$, and $\mathit{\Gamma }$ are time-invariant. Among these, ${\mathit{K}}_{\mathit{e}}$, ${\mathit{\Phi }}_{\mathit{e}}$, and $\mathit{\Gamma }$ are strongly dependent on the mass and stiffness distribution of the structures. The product ${\mathsf{\Sigma }}_{\mathit{e}}={\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{e}}{\mathit{\Phi }}_{\mathit{e}}$ determines the impact of the i-th mode on the r-th component of the load for each element, where $rϵ\left\{1:6\right\}$ and $\mathit{\Gamma }\mathit{d}\left(t\right)$ determines the amplitude at which the i-th mode oscillates.

To compare this definition to the second class of participation factors, we introduce the following modal participation fraction:

$$L_{r,i}^{load}=\frac{\left|{{\Sigma }_e}_{r,i}{\Gamma }_i{d}_i\right|}{{\sum }_{i=1}^N\left|{{\Sigma }_e}_{r,i}{\Gamma }_i{d}_i\right|}$$

(18)

The term $d_i$ is written without an upper index in Equation (18) to indicate that any mathematical form used to solve for $d_i$ may be employed. It should be noted that the purpose behind defining the modal participation fractions is to compare the outcomes of various mathematical formulations of “modal participation factors”.

2.2.2. Class II: Modal Relative Effectiveness Factors (MREFs)

The second class of modal participation factors aims to identify mode shapes that are easily excited. Many of these factors employ equations to determine a constant parameter that offers insight into the system’s behavior without the need to solve any differential equations of motion.

The modal effective mass participation factor, a prominent member of the second-class classification, plays a crucial role in selecting the optimal number of modes needed to represent the physical domain [31,55]. Various approaches have been proposed to define these factors, as discussed in the literature. For instance, Wilson [38] introduced a modal participation factor akin to Equation (7) to identify influential modes for advanced analyses.

The effective mass participation factor has found applications in diverse fields, including the study of vibration control in air thrusters [56]. In one application, a cantilever beam was analyzed using a finite element model (FEM), with model reduction via the mode displacement method for continuous control assessment [57]. The versatility of FEM techniques has been demonstrated in numerous mechanical analysis studies. For instance, it was applied to investigate CNT-reinforced composite plate and shell structures, encompassing both static [58] and dynamic [59] analyses with a director FEM model approach. These structures exhibited various functional gradations, and the accuracy of results was validated against literature values.

Additionally, a similar model based on nonlinear director FEM was developed for studying nonlinear bending [60,61] and large deflection [62] in CNT- and graphene-reinforced nanocomposite shell structures. Parametric studies were conducted to explore the design aspects of reinforcements and their influence on structural responses across various case studies, including beams, panels, and shells.

Modal Effective Mass Participation Factor

In the initial stages, the modal effective mass was calculated for a single-DoF system subjected to a base acceleration. It was observed that the effective mass of the system varies directly with the reaction force, s. Thus, it can be inferred that a significant inertial force is present when a forced acceleration $\ddot{\mathrm{a}}$ acts on a modal effective mass of considerable magnitude.

A Craig–Brampton transformation [55,63] can be utilized to extend the above theory to multi-DoF systems which produce Equation (19) in compact notation.

$$\mathit{s}={\mathit{M}}^{\mathit{e}\mathit{f}\mathit{f}}\left[\mathbf{1}+{\left(\frac{\omega }{{\omega }_n}\right)}^2\mathit{H}\left(\frac{\omega }{{\omega }_n}\right)\right]\ddot{\mathrm{a}}$$

(19)

where a transfer function parameter is represented by $\mathit{H}\left(\frac{\omega }{{\omega }_n}\right)$ and ${\mathit{M}}^{\mathit{e}\mathit{f}\mathit{f}}$ represents the matrix of the modal effective masses. Each component of the matrix ${\mathit{M}}^{\mathit{e}\mathit{f}\mathit{f}}$ is found using the following expression:

$$m_{p,i}^{eff}=\frac{{\Gamma }_{p,i}^{mm}\times {\Gamma }_{i,p}^{mm}}{{m_q}_i}$$

(20)

This quantity represents the effective mass of an elastic mode when coupled with a rigid-body mode. The mass of each rigid-body degree of freedom, denoted as p, is equal to the sum of masses associated with each elastic mode, i, for that rigid-body mode, p.

The constant ${\mathit{\Gamma }}^{\mathit{m}\mathit{m}}$ denotes the matrix of modal effective mass participation factors of size $\left(N-p\right)\times p$(elastic-modes $\times$ rigid-modes), and ${m_q}_i$ is the modal mass of the i-th eigenvector in Equation (20). The former relates the rigid-body modes ${\mathit{\Phi }}^{\mathbf{r}}$ of size $n\times p$(nodes $\times$ rigid-modes) with the matrix ${\mathit{\Phi }}^{\mathit{e}}$ indicating elastic modes and of size $n\times (N-p)$ (nodes $\times$ elastic-modes) in compact notation as follows:

$${{\mathit{\Gamma }}^{\mathit{m}\mathit{m}}}^T={{\mathit{\Phi }}^{\mathit{r}}}^T\mathit{M}{\mathit{\Phi }}^{\mathit{e}}$$

(21)

An important consideration in defining the modal effective mass, even in a free–free multi-DoF system, is to examine the elastic degrees of freedom in relation to a specific degree of freedom (DoF) known as the boundary DoF. The analysis requires constraining this boundary DoF, a vital step in determining the modal characteristics of the system. Ignoring this aspect can lead to inaccurate results when using Equation (20) because, for free–free systems, the contribution of elastic DoFs tends to approach zero.

Within the context of the Craig–Bampton transformation, this important point is intricately connected. For more details, readers are referred to [31,55]. Lastly, the effective modal mass fraction is defined as the ratio of the effective modal mass to the total system mass, expressed in component form, as follows:

$$L_{p,i}^{mm}=\frac{m_{p,i}^{eff}}{{\sum }_{j=1}^n{m}_n}$$

(22)

Free–Free Modal Participation Factors

Based on the concept of effective interface mass (EIM), Kammer et al. [64] recently introduced a new definition to assess the relative importance of each mode for free–free systems.

By introducing the variable change as per Equation (2) into Equation (1) in compact notation and solving for the physical acceleration, the following expression is obtained:

$$\ddot{\mathit{u}}\left(t\right)=\mathit{\Phi }\ddot{\mathit{q}}\left(t\right)=\mathit{\Phi }{\mathit{\Phi }}^{\mathit{T}}\mathit{f}p\left(t\right)-\mathit{\Phi }\mathit{\lambda }\mathit{q}\left(t\right)$$

(23)

The product $\mathit{\Phi }{\mathit{\Phi }}^{\mathit{T}}$ is conventionally considered as the connection between the response of the j-th node when the n-th node is exposed to an excitation in structural dynamics. By taking into account the product $\mathit{\Phi }{\mathit{\Phi }}^{\mathit{T}}\mathit{f}$, this interpretation can be strengthened, and as a result, it is anticipated that the modes that are excited by the external force vector will significantly add to the response of the system.

The trace of the product between two eigenvectors and the free–free modal participation factor assumes the shape of Equation (24) in component form to evaluate the above-discussed interactive contribution [64].

$${\Gamma }_i^{free}=tr\left({\mathit{\phi }}_{\mathit{i}}{{\mathit{\phi }}_{\mathit{i}}}^{\mathit{T}}\right)$$

(24)

A free–free modal participation fraction is the ratio between the free–free participation factors and the trace of the elastic degrees of freedom, ${\Gamma }^{elas}$, which is defined using the output of Equation (24) and aids in forming a comparative analysis with the other formulations discussed.

$$L_i^{free}=\frac{{\Gamma }_i^{free}}{{\Gamma }^{elas}}=\frac{tr\left({\mathit{\phi }}_{\mathit{i}}{{\mathit{\phi }}_{\mathit{i}}}^{\mathit{T}}\right)}{tr\left({\mathit{\Phi }}^{\mathit{e}}{{\mathit{\Phi }}^{\mathit{e}}}^{\mathit{T}}\right)}$$

(25)

Modal Strain Energy Participation Factor

Here, Equation (26) is the widely adopted formula for strain energy SE(t).

$$SE(t)=\frac{1}{2}{\mathit{u}\left(t\right)}^{\mathit{T}}\mathit{K}\mathit{u}(t)$$

(26)

The term ${\mathit{K}}_{\mathit{q}}$, which denotes the matrix of modal stiffnesses, emerges when the strain energy time history is decomposed into its modal components. The diagonal eigenvalue matrix and the modal stiffness matrix are identical if the eigenvectors are mass-normalized. A matrix of peak strain energies ${\mathit{S}\mathit{E}}^{\mathit{m}\mathit{a}\mathit{x}}$, of dimension $n\times N$, may be constructed as follows, assuming that the peak value is the object of interest:

$${\mathit{S}\mathit{E}}^{\mathit{m}\mathit{a}\mathit{x}}=\frac{1}{2}{{\mathit{q}}^{\mathit{m}\mathit{a}\mathit{x}}}^{\mathit{T}}{\mathit{\Phi }}^{\mathit{T}}\mathit{K}\mathit{\Phi }{\mathit{q}}^{\mathit{m}\mathit{a}\mathit{x}}=\frac{1}{2}{{\mathit{q}}^{\mathit{m}\mathit{a}\mathit{x}}}^{\mathit{T}}{\mathit{K}}_{\mathit{q}}{\mathit{q}}^{\mathit{m}\mathit{a}\mathit{x}}$$

(27)

Since the highest magnitude of the harmonic modal responses must be solved, the strain energy matrix in Equation (27) can be complex and expensive to assess. Instead, several researchers [20,65,66], and have suggested using the element modal strain energy ${MSE}_{ji}$ as a criterion for evaluation. The global coordinate system is used to describe and make up the global stiffness matrix K, and the ${MSE}_{ji}$ assesses the participation of each mode in each element $j$ stiffness submatrix ${\mathit{K}}_{{\mathit{e}}_{\mathit{j}}}$,

$${MSE}_{ji}=\frac{1}{2}{\mathit{\phi }}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{{\mathit{e}}_{\mathit{j}}}{\mathit{\phi }}_{\mathit{i}}$$

(28)

According to [67], a modal strain participation fraction, $L_{ji}^{MSE}$, can be described as the relationship between kinetic energy and strain energy when a system is in free oscillation. Thus, the following relationship can be deduced:

$$L_{ji}^{MSE}=\frac{{\mathit{\phi }}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{{\mathit{e}}_{\mathit{j}}}{\mathit{\phi }}_{\mathit{i}}}{{\mathit{\phi }}_{\mathit{i}}^{\mathit{T}}\mathit{M}{\mathit{\phi }}_{\mathit{i}}{\omega }_{n_i}^2}$$

(29)

A high value of $L_{ji}^{MSE}$ for the i-th mode in the j-th element would indicate that any variation in the i-th mode shape would produce a significant change in the value of the elemental internal load.

3. Sensitivity of Load via Static Modal Participation Factors

In structural optimization, the dynamic behavior of a system is governed by various parameters. To achieve minimum weight and maximum performance in functional structures, these parameters are often adjusted. Gradient-based optimization methods, commonly used, require the computation of costly derivatives for cost and constraint functions. Typically, finite difference methods [68] are employed for this purpose, but this approach can limit the practical application of optimizers.

Alternatively, the variation in dynamic loads can be approximated using the SMPF method, as discussed in earlier sections. Various finite element models (FEMs) and numerical methods [69,70] have also been developed for applications demanding high-performance design parameters with a flexible structural modeling strategy.

Complex scenarios, such as the mechanical buckling behavior of fractionally graded reinforcement in nanocomposites in the form of curved panels and plates, have been analyzed [71]. This analysis includes determining influential parameters using a discrete finite element shell model, which defines a displacement field and considers the effect of transverse shear deformation. The same model has been extended to include a first-order shear deformation theory for post-buckling analysis of functionally graded nanocomposite cylindrical shells and plates under mechanical [72] and thermo-mechanical [73] load cases.

Techniques like hybrid structural modeling and the momentum-based approach known as direct matrix input at a grid point (DMIG) are efficient methods for load recovery in load-bearing structures. These methods differ from the deformation approach used in this study and have been shown to yield lower values of modal truncation errors compared to deformation-based approaches and other load recovery methods.

Considering the time-invariant system as given by Equation (9), a small perturbation $∆$ can be introduced into expression 17 in compact notation, yielding

$${\mathit{l}}_{\mathit{e}}^{\mathbf{0}}+\mathsf{\Delta }{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}={\mathit{C}}^{\mathit{T}}\left({\mathit{K}}_{\mathit{e}}+\mathsf{\Delta }{\mathit{K}}_{\mathit{e}}\left)\right({\mathit{\Phi }}_{\mathit{e}}+\mathsf{\Delta }{\mathit{\Phi }}_{\mathit{e}}\right)\left(\mathit{\Gamma }+\mathsf{\Delta }\mathit{\Gamma }\left)\right({\mathit{d}}^{\mathbf{0}}+\mathsf{\Delta }{\mathit{d}}^{\mathbf{0}}\right)$$

(30)

Assuming that the optimizer-induced changes are parametric and sufficiently small, we can disregard higher-order terms, resulting in Equation (31). This straightforward expression relates the change in the magnitude of the maximum incremental dynamic load $\mathsf{\Delta }{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}$ due to small perturbations in the stiffness matrix $\mathsf{\Delta }{\mathit{K}}_{\mathit{e}}$, the eigenvectors $\mathsf{\Delta }{\mathit{\Phi }}_{\mathit{e}}$, and the peak dynamic response of the i-th mode for the unitary forcing function $\mathsf{\Delta }{\mathit{d}}^{\mathbf{0}}$.

$$\begin{gathered} \mathit{\Delta }{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}={\mathit{C}}^{\mathit{T}}\left({\mathit{K}}_{\mathit{e}}\mathit{\Delta }{\mathit{\Phi }}_{\mathit{e}}+\mathit{\Delta }{\mathit{K}}_{\mathit{e}}{\mathit{\Phi }}_{\mathit{e}}\right)diag\left({\mathit{\Phi }}^{\mathit{T}}\mathit{f}\right){\mathit{d}}^{\mathbf{0}}+ \\ {\mathit{C}}^{\mathit{T}}{\mathit{K}}_{\mathit{e}}{\mathit{\Phi }}_{\mathit{e}}\left(diag\left({\mathit{\Phi }}^{\mathit{T}}\mathit{f}\right)\mathit{\Delta }{\mathit{d}}^{\mathbf{0}}+diag\left({\mathit{\Delta }\mathit{\Phi }}^{\mathit{T}}\mathit{f}\right){\mathit{d}}^{\mathbf{0}}\right) \end{gathered}$$

(31)

Since Equation (31) allows the quick estimation of the incremental load, this expression can be used directly in the finite difference equation. By neglecting the contribution of higher-order terms of the Taylor series expansion, the following expression is obtained:

$${{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}}^{\mathbf{\prime }}=\frac{{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}\left(z_0+∆z\right)-{\mathit{l}}_{\mathit{e}}^{\mathbf{0}}\left(z_0\right)}{∆z}=\frac{\Delta {\mathit{l}}_{\mathit{e}}^{\mathbf{0}}}{∆z}$$

(32)

where $z_0$ represents an initial design variable and $∆z$ is its corresponding incremental value.

4. Proof of Concept

A proof of concept is demonstrated using a cantilever beam model with 5 degrees of freedom (DoFs), as illustrated in Figure 1. This example aims to help readers understand the outcomes of each definition discussed in the paper and facilitate comparative analysis. The model is subjected to a periodic dynamic load, $p\left(t\right),$ applied to its free end, as defined by Equation (33), where $\omega$ represents the forcing frequency and $t$ denotes time.

$$p\left(t\right)=sin\left(\omega t\right)$$

(33)

The modal response of the system, subjected to a load as described in Equation (3) and assuming an under-damping damping ratio,$\xi$, of 5% for each mode, is numerically computed using the fourth-order Runge–Kutta method [74]. The resulting time history is presented in Figure 2a, clearly demonstrating the dominance of the first harmonic mode. In Figure 2b, 2c, and 2d, we illustrate the modal contributions to the physical response at degrees of freedom 1, 3, and 5, respectively, as calculated based on Equation (3). These figures also include the superimposed total nodal response. It is important to note that the contribution of the i-th mode varies for each j-th nodal response.

To estimate the modal participation fraction of the system at hand using Equations (14) and (16), an appropriate expression for $d_i\left(t\right)$ must be derived. This can be easily accomplished by the classic method of undetermined parameters, where the resulting expression becomes

$$d_i\left(t\right)=\frac{1}{\left[{\omega }_{n_i}^4-{\omega }^4+{\left(2\xi {\omega }_{n_i}\omega \right)}^2\right]}\left\{\begin{array}{c}{e}^{-\xi {\omega }_{n_i}t}\left[2\xi {\omega }_{n_i}\omega Cos\left({\omega }_{d}t\right)-\frac{\omega }{{\omega }_d}\left({\omega }_{n_i}^2-2{\xi }^2{\omega }_{n_i}+{\omega }^2\right)Sin\left({\omega }_{d}t\right)\right]\\ +\left({\omega }_{n_i}^2+{\omega }^2\right)Sin\left(\omega t\right)-2\xi {\omega }_{n_i}\omega Cos\left(\omega t\right)\end{array}\right\}$$

(34)

where ${\omega }_d$ is the damped frequency and is expressed as

$${\omega }_d={\omega }_{n_i}\sqrt{\left|{\xi }^2-1\right|}$$

(35)

To determine the peak amplitude of $d_i$, the saddle points of Equation (34) are located by satisfying the first-order necessary condition [75], which yields the following expression:

$$\begin{array}{c}{e}^{-\xi {\omega }_{n_i}t}\left[\frac{\xi {\omega }_{n_i}}{{\omega }_d}Sin\left({\omega }_{d}t\right)\left({\omega }_{n_i}^2-2{\xi }^2{\omega }_{n_i}+{\omega }^2-2{\omega }_d^2\right)+Cos\left({\omega }_{d}t\right)\left(2{\xi }^2{\omega }_{n_i}-{\omega }^2-{\omega }_{n_i}^2\left(1+2{\xi }^2\right)\right)\right]\\ =\dots \dots \left({\omega }_{n_i}^2+{\omega }^2\right)Cos\left(\omega t\right)-2\xi {\omega }_{n_i}\omega Sin\left(\omega t\right)\end{array}$$

(36)

The peak amplitude of the multimodal Equation (34) can be determined by locating the roots of Equation (36). It is worth noting that the peak response occurs during the transient period, so only the roots within this period are relevant. Finally, by applying Equations (14) and (16) in conjunction with Equations (34)–(36), the participation fraction contribution associated with each mode during the transient period can be computed.

The estimated modal participation fractions using modes 12, 14, 16, 22, and 25 are presented in

Table 1. Participation factors for an under-damped system subjected to a sinusoidal forcing function.

Mode	${\mathit{L}}_{\mathit{i}}^{\mathit{s}\mathit{t}}$ (12)	${\mathit{L}}_{\mathit{i}}^{\mathit{t}\mathit{s}\mathit{I}}$ (14)	${\mathit{L}}_{\mathit{i}}^{\mathit{t}\mathit{s}\mathit{I}\mathit{I}}$ (16)	${\mathit{L}}_{\mathit{i}}^{\mathit{m}\mathit{m}}$ (22)	${\mathit{L}}_{\mathit{i}}^{\mathit{f}\mathit{r}\mathit{e}\mathit{e}}$ (25)
1	0.8558	0.8627	0.8658	0.8819	0.1993
2	0.0934	0.0890	0.0872	0.0851	0.1996
3	0.0318	0.0302	0.0294	0.0225	0.2029
4	0.0136	0.0129	0.0126	0.0088	0.1987
5	0.0054	0.0051	0.0050	0.0016	0.1995

. The steady-state $L_i^{st}$, the transition $L_i^{tsI}$, and the time-consistent $L_i^{tsII}$ participation fractions are estimates of the amplitude of the modal response. As Figure 2a illustrates, these percentages of participation align well with the observed trend: the first mode exhibits a significantly larger amplitude compared to the fifth mode at the given excitation frequency.

The results in Figure 2a also align with the modal effective mass participation factor. A high modal effective mass indicates that a particular mode possesses a substantial modal inertial component, rendering it susceptible to excitation. However, it is worth noting that this metric remains constant for a given system and is independent of the excitation frequency. Consequently, while the modal effective mass is informative, it may not be the most suitable indicator for assessing modal response values, as these responses are influenced by the ratio of natural to excitation frequencies.

The various free–free participation fractions discussed differ from the remainder of the computed participation factors simply because the results imply that the contribution of all the modes is similar. This is because the modal residue is comparable to the trace product of the two equal eigenvectors. In other words, it depicts the reaction of the j-th node to an applied force on the r-th node. As a result, this definition does not offer a true evaluation of a mode’s relevance but instead evaluates the reaction in a node, j, when a force is applied.

The maximum value of the elemental shear force can be calculated based on Equation (17) and the prior values of the steady-state $L_i^{st}$, the transition $L_i^{tsI}$, and the time-consistent $L_i^{tsII}$ participation percentages.

Table 2. Estimation of the peak shear force for an under-damped system subjected to a sinusoidal forcing function.

Element	Actual Peak Load (N)	$\mathbf{Load}{\mathit{L}}_{\mathit{i}}^{\mathit{s}\mathit{t}}$ (N)	% Error	$\mathbf{Load}{\mathit{L}}_{\mathit{i}}^{\mathit{t}\mathit{s}\mathit{I}}$ (N)	% Error	$\mathbf{Load}{\mathit{L}}_{\mathit{i}}^{\mathit{t}\mathit{s}\mathit{I}\mathit{I}}$ (N)	% Error
1	53.1225	48.9718	−7.8134%	57.9858	9.1550%	53.0976	−0.0467%
2	53.7885	49.9985	−7.0461%	45.8251	−14.8049%	53.8574	0.1282%
3	51.5479	48.6080	−5.7032%	49.1496	−4.6525%	51.5973	0.0959%
4	51.3453	49.3883	−3.8114%	65.4091	27.3906%	51.3164	−0.0563%
5	51.2290	50.5989	−1.2300%	38.6055	−24.6413%	51.0221	−0.4039%

presents the findings of the analysis, where the most accurate results were recorded using the time-consistent $L_i^{tsII}$.

The error in the calculation of $L_i^{st}$ is owing to the exclusion of the transient components, and in a real scenario, the load recovered corresponds to only the actual steady-state shear force. The error value related to $L_i^{tsI}$ on the other hand, was observed to be much higher. There is an over/underestimation in the calculation of the maximum response value, and therefore, the peak load is due to some fundamental assumptions such as the assumption that the modes are in phase.

Finally, by combining Equation (18) with Equations (11), (13) and (15), we computed the participation fractions associated with individual modes in response to the internal load. Additionally, Equation (29) was employed to calculate the modal strain energy. To facilitate comparison, all datasets were normalized relative to their largest magnitude component. The results are presented in Figure 3.

Figure 3 displays the normalized load participation factors for elements 1–5 of the cantilever beam, as shown in Figure 3a–e, respectively. These results confirm that the estimated participation factors obtained using various approaches described in Equations (11), (13) and (15) in conjunction with Equation (18) exhibit excellent agreement. However, it is noteworthy that the mode participation calculated through the modal strain energy method exhibits considerable variation.

Equation (17) is a product of two terms: first, ${\Gamma }_i{d}_i$ is a measure of the amplitude of the modal response; second, ${{\Sigma }_e}_{r,i}$ reflects the difference in the magnitude of the relative amplitudes of the nodes associated with a certain mode. The smaller the disparity in relative amplitudes, the less the mode contribution to the entire spectrum. To illustrate this, examine the load distribution in element 1, as shown in Figure 3, which is influenced by both the first and second modes. In contrast, element 5 still exhibits the dominance of the first mode, albeit to a lesser extent, while other modes become more influential. Consequently, any structural modification that impacts the behavior of the first mode will exert a substantial influence on the load experienced by element 1, but its effect on element 5 will be comparatively smaller.

In Figure 3, we observe the modal strain energy, which indicates that higher-order modes make more substantial contributions to the overall response. It is important to note that the definition in Equation (29) represents the ratio of elemental modal strain energy to the modal kinetic energy of the system. As a result, a lower value implies that the kinetic energy, representing the inertial component, has a greater impact on the response compared to the quasi-static elemental modal strain energy.

For example, element 5 contains a significant third-mode modal strain energy. Consequently, any change to the third mode shape (${\mathit{\phi }}_{\mathbf{3}}$) would significantly affect the response of element 5. Conversely, any change to the shape of the first mode (${\mathit{\phi }}_{\mathbf{1}}$) would not record a significant influence on its response. However, this would not be the case if a variation was made to the natural frequency (${\omega }_{n_1}$); a substantial change in the system’s kinetic energy and, therefore, an effect on the load-acting element 5 can be predicted.

Approximate Derivative of Element Internal Load via SMPF

A series of mass and stiffness parametric variations were performed on the 5-DoF model. The loads were recovered after each modification using the classical solution method (i.e., solving the differential equations of motion in the modal domain) and the SMPF method with the time-consistent $L_i^{tsII}$ approach.

The partial derivative of the peak load concerning the variation in the first eigenvalue was estimated using the traditional finite difference equation with truncated error. Alternatively, it was computed using Equation (32) in conjunction with Equation (15).

The numerical results are presented in

Table 3. Comparison of the partial derivatives obtained using the traditional finite difference equation versus those computed via SMPF.

It	Element 1			Element 2			Element 3			Element 4			Element 5
	Finite Diff	${{\mathit{l}}_{\mathit{e}}^{\mathit{m}\mathit{a}\mathit{x}}}^{\prime }$	Error	Finite Diff	${{\mathit{l}}_{\mathit{e}}^{\mathit{max}}}^{’}$	Error	Finite Diff	${{\mathit{l}}_{\mathit{e}}^{\mathit{m}\mathit{a}\mathit{x}}}^{\prime }$	Error	Finite Diff	${{\mathit{l}}_{\mathit{e}}^{\mathit{max}}}^{’}$	Error	Finite Diff	${{\mathit{l}}_{\mathit{e}}^{\mathit{max}}}^{’}$	Error
1	−1.730	−1.730	0.000	−2.820	−2.820	0.000	−1.627	−1.627	0.000	−2.662	−2.662	0.000	−1.646	−1.646	0.000
2	−1.616	−1.730	0.114	−2.880	−2.820	−0.060	−1.692	−1.627	−0.065	−2.685	−2.662	−0.022	−1.888	−1.646	−0.242
3	−1.853	−1.734	−0.119	−2.645	−2.846	0.201	−1.683	−1.627	−0.056	−2.566	−2.684	0.118	−1.889	−1.644	−0.245
4	−1.873	−1.734	−0.139	−3.116	−2.869	−0.247	−1.630	−1.624	−0.006	−2.963	−2.703	−0.260	−1.894	−1.641	−0.253
5	−1.403	−1.731	0.328	−2.763	−2.891	0.128	−1.332	−1.619	0.287	−2.806	−2.721	−0.086	−1.944	−1.636	−0.308

for all five elements, where good agreement is found between the two methods. The error associated with the use of Equation (32) in combination with Equation (15) is due, firstly, to the accumulated error of the approximated peak response using $L_i^{tsII}$ and, secondly, to the assumption that the changes in the mass and stiffness distribution are sufficiently small to truncate higher-order terms of Equation (30).

5. Case Study

As demonstrated in the initial proof of concept, the second class of participation factors can be employed to rapidly assess the relative modal contributions within the solution. Additionally, the SMPF method could serve as a re-analysis technique for approximating the peak physical response without the need to solve differential equations.

Therefore, to further exploit the latter concept, parametric mass and stiffness variations were performed in an airframe stick model [76], and the after-modification critical loading conditions were found via the SMPF method. In total, 26 parametric variations were considered, and only the time-consistent $L_i^{tsII}$ participation fractions were used due to their previously proven accuracy. The reduced-order stick models are already proven to be consistent with the results of full-order finite element (FE) in several structural load cases for analyzing the responses of specialized cases such as self-elevating units (SEUs) [77]. Several studies have demonstrated techniques in practical multidisciplinary design [4] for networks of aircraft components, including an automated load process with a two-level load approach, including aeroelastic analysis and sensitivity studies of the influential parameters. In certain advanced studies, a discrete director-based finite element (FE) formulation was also shown to be robust for free vibration analysis of various carbon nanotube (CNT)-reinforced structures such as composite shells, cylindrical panels, and plates [78], including cases of functional grading for material distribution and thermal effects [79].

For simplicity, the dynamic excitation vector, $\mathit{f}$, is represented as a linearly distributed force over the wing span, as seen in Figure 4, starting with an applied force of 1000 N at the wing root node (located at a normalized distance $\overline{y}=0.0$) and finishing with a force of 480 N at the wing tip node (located at a normalized distance $\overline{y}=1.0$). The time history, $p\left(t\right)$, of the dynamic excitation is given by Equation (37), which in the context of aeronautics, is representative of a wind gust [80].

$$p\left(t\right)=\left\{\begin{array}{cc}1-cos\left(\omega t\right)& 0\le t\le \frac{2\pi }{\omega }\\ 0& t\ge \frac{2\pi }{\omega }\end{array}\right.$$

(37)

The model was subjected to a family of “gusts” with various wavelengths, ranging from 12.5 to 31.5 radians, from which the resulting dynamic peak loads had to be identified. In particular, these loads were the peak out-of-plane bending moment (OPB moment), the peak out-of-plane shear force (OPS force), and the peak torsion moment at four monitored elements located at wing spans of $\overline{y}=0.00$, $\overline{y}=0.33$, $\overline{y}=0.66$, and $\overline{y}=1.00$.

To approximate the structures’ peak dynamic response, an appropriate expression for $d_i\left(t\right)$ was derived via the impulse response function, which yields the following equation:

$$d_i\left(t\right)=\left\{\begin{array}{c}\frac{1}{{\omega }_{n_i}^2}\left(1-cos\left({\omega }_{n_i}t\right)\right)-\frac{2}{{\omega }^2-{\omega }_{n_i}^2}sin\left(\frac{\omega +{\omega }_{n_i}}{2}t\right)sin\left(\frac{\omega +{\omega }_{n_i}}{2}t\right)0\le t\le \frac{2\pi }{\omega }\\ \begin{array}{c}\frac{1}{{\omega }_{n_i}^2}\left(cos\left({\omega }_{n_i}\left(t-\tau \right)\right)-cos\left({\omega }_{n_i}t\right)\right)+\frac{2}{2{\omega }_{n_i}\left(\omega -{\omega }_{n_i}\right)}\left(cos\left(\left(\omega -{\omega }_{n_i}\right)\tau +{\omega }_{n_i}t\right)-cos\left({\omega }_{n_i}t\right)\right)-\dots \\ \dots \frac{2}{2{\omega }_{n_i}\left(\omega +{\omega }_{n_i}\right)}\left(cos\left(\left(\omega +{\omega }_{n_i}\right)\tau -{\omega }_{n_i}t\right)-cos\left({\omega }_{n_i}t\right)\right)t\ge \frac{2\pi }{\omega }\end{array}\end{array}\right.$$

(38)

The saddle points of the above equation are located when the first-order necessary condition is satisfied, which provides the following expression:

$$0=\left\{\begin{array}{c}sin\left({\omega }_{n_i}t\right)-\left(\frac{{\omega }_{n_i}}{\omega }\right)sin\left(\omega t\right)\\ \begin{array}{c}\frac{2}{{\omega }_{n_i}}cos\left({\omega }_{n_i}\left(t-\frac{\pi }{\omega }\right)\right)sin\left(\frac{{\omega }_{n_i}}{\omega }\pi \right)+\left(\frac{1}{\left(\omega -{\omega }_{n_i}\right)}\right)cos\left(1+{\omega }_{n_i}\left(t-\frac{1}{\omega }\right)\right)sin\left(\frac{{\omega }_{n_i}}{\omega }\pi -1\right)-\dots \\ \dots \left(\frac{1}{\left(\omega +{\omega }_{n_i}\right)}\right)sin\left({\omega }_{n_i}t+\pi \left(1+\frac{{\omega }_{n_i}}{\omega }\right)\right)cos\left(\pi \left(1+\frac{{\omega }_{n_i}}{\omega }\right)\right)\end{array}t\ge \frac{2\pi }{\omega }\end{array}\right.$$

(39)

To estimate the peak amplitude of $d_i\left(t\right)$, the roots of Equation (39) are found for the dominant mode. Expression (38) is evaluated to find $t_{dom}^{max}$, and the magnitude of remaining modes retained is estimated using Equation (38) at $t_{dom}^{max}$. Finally, Equation (17) is used to recover the element peak loads in the physical domain.

The normalized magnitude of various peak loads obtained using the improved SMPF method is shown in Figure 5, Figure 6, Figure 7 and Figure 8 and compared to results obtained using Nastran’s solution 112. Again, excellent agreement is found between the two numerical results.

The root mean square (RMS) error of the torsional moment at the wing root is notably the highest, at a magnitude of 2.99 percent, primarily due to the significant contribution of high-order modes. It is worth noting that Equation (15) assumes that the peak physical response occurs at $t_{dom}^{max}$ (i.e., the time instant where the peak of the single dominant mode occurs). However, when higher-order modes start to play a more substantial role, the peak physical response may occur close to but not exactly at $t_{dom}^{max}$. Consequently, this assumption may fail to accurately capture the maximum physical response in cases where two or more modes are equally dominant in the solution.

Approximate Derivative of Element Internal Load via SMPF

The approximated sensitivity of the OPB moment, the OPS force, and the torsion moment with respect to the variation in the first eigenvalue are presented in Figure 9, Figure 10, Figure 11 and Figure 12 for the four elements of interest. The RMS error of each plot is presented as a legend at the bottom of each figure. The torsion moment exhibits the highest error among them.

While some outliers are noticeable, there is generally good agreement between the finite difference method and Equation (32). These outliers are likely caused by the omission of higher-order terms in Equation (31) leading to a loss of information.

6. Conclusions

In this paper, the concept of modal participation factors (MPFs) is proposed as an alternative to the current point-wise treatment of the dynamic equations of motion of flexible structures. An in-depth literature review was conducted to determine the MPFs that are most relevant to structural dynamics and can be employed as point-wise constraint equations during structural design. It was found that though the second class of participation factors can be used to evaluate the relative contribution of a mode within the solution inexpensively, it is the SMPF method (which was discussed under the first class of participation factors) that can potentially be utilized in the point-wise constraint equations during the numerical optimization of dynamic systems. The advantage of the method lies in the existence of a closed-form solution and the dominance of a single mode. However, in cases where two or more modes are equally relevant in the solution, each of their contributions must be accounted for before determining the time step at which the peak response (i.e., maximum or minimum) occurs.

The sensitivities of the peak transient loads can be approximated using the SMPF method without the need to solve any differential equations, thus reducing the computational expense. However, the associated truncated higher-order error terms may introduce inaccuracies in the approximation. This error compounds the well-known inaccuracies of the finite difference method when the step size (the mass and stiffness variations, in this case) is not sufficiently small.

Finally, the second class of modal participation factors, especially the modal strain energy, can be considered as a useful pre-analysis tool. These metrics can assist decision making in traditional structural sizing, relying on the engineer’s judgment. During multi-objective optimization cycles, the development of model updating techniques can also be invaluable for analyzing the fluctuations in error within the dynamic characteristics of the studied structure. This approach holds promise, particularly when supported by resources like supercomputing, for analyzing large-scale structures with a substantial number of degrees of freedom and complex computational models.