Vibration and Optimal Control of a Composite Helicopter Rotor Blade

1. Introduction

Vibration is an inherent phenomenon in all dynamic systems and generally manifests as an oscillatory motion arising from the system’s mass and relative component interactions. Regardless of the system complexity, these oscillations act as disturbances that dissipate energy, reduce efficiency, and shorten service life of the system components [1]. In automotive and aerospace vehicles, vibrations contribute significantly to structural fatigue, acoustic noise, and passenger discomfort, ultimately compromising operational reliability. Therefore, effective vibration control is critical to ensure component durability and optimize system performance. The primary objective of vibration control is to mitigate unwanted oscillations induced either by external disturbances or internal structural imbalances and remains challenging to select the most effective control strategy [2]. In helicopter rotor systems, unlike hovering flight, forward flight introduces asymmetric aerodynamic loading [3] owing to variations in the relative airflow velocity between the advancing and retreating blades. This generates periodic aerodynamic loadings at the blade root leaving residual unbalanced components that excite harmonics near the blade-passing frequency, potentially coupling with structural modes and leading to resonance. To address these challenges, researchers have explored effective control methodologies and actuation systems aimed at achieving vibration suppression at desired levels under diverse operating conditions including unsteady aerodynamics and gust loadings [4].

Previous studies have shown the effort of researchers provided for vibration control using various strategies including the application of smart materials. A finite element model [5] was proposed by Lam et al. [6] for piezoelectric composite laminates, where a velocity feedback control algorithm was used. Chuaqui et al. [7] used the local collocation method based on the finite- difference generated radial basis functions to implement a constant-gain velocity feedback active vibration control for cantilever piezoelectric smart beams. The effective use of composite materials with piezoelectric sensors and actuators attached or embedded between composite layers can also be used to suppress vibrations [8,9,10]. All these studies emphasized the need to understand active vibration control in continuous systems through a combination of analytical modeling, state estimation, and control optimization.

Active vibration control is categorized into higher harmonic control (HHC), individual blade control (IBC), and active control of structural response [11,12]. Among them, HHC and IBC involve control of blade pitch or trailing edge flap deflection to generate unsteady, higher harmonic loads that cancel the original vibration. In HHC, actuators are placed below the swash plate; hence, the mechanically applicable control frequencies are limited to three blades [11]. The IBC is a suitable method to overcome the limitation of the HHC [12]. In IBC, the actuators are attached to the rotating frame based on the same HHC algorithm. As each rotor blade is actuated individually, it provides more flexibility to control undesirable dynamic phenomena and hence, with the passage of time, the design of IBC has evolved more with the emergence of smart materials.

Experimental and numerical studies have further advanced active vibration control techniques demonstrating their effectiveness across rotorcraft and other complex dynamic systems. Park et al. [13] investigated the performance of a lift-offset rotor with IBC and found that the rotor hub vibration for high-speed flight could be reduced significantly with multiple harmonic inputs. Experimental studies involving IBC, promising for helicopter vibration attenuation and gust alleviation are highlighted in the research study conducted by Yang et al. [14]. They showed that under the effect of high harmonic input, the vibratory loads on the rotor hub could be reduced by 65% for forward flights. A wind-tunnel study of two rotors without actuators in the rotating frame in a Mach-scaled model revealed that the IBC had potential to reduce the vibration up to 77% in the descent flight [15]. Following this, full-scale rotor blade analysis was also conducted for the four-bladed Bo 105 helicopter rotor to investigate the IBC system at the NASA Ames Research Center showing simultaneous reduction in vibration and blade–vortex interaction noise [16,17]. Chia et al. [18] extended these efforts by developing the HELINOIR computational framework for predicting and controlling low-frequency rotorcraft noise using active flaps, which was validated against the wind-tunnel results.

To apply the control strategy successfully on the rotor blade, the vibratory load on the rotor blade needs to be properly estimated which depends on the corresponding aerodynamic environment for different flight conditions, coupling among the various degrees-of-freedom (DOFs) of motion of the rotor blade, blade structure, and relevant material properties. To estimate the vibratory response of the rotor blade accurately, the vibration frequencies need to be evaluated precisely. Besides the coupling between bending and torsional vibrations caused by noncoincident inertial and elastic axes in composite rotor blades, different ply orientations in the composite layers of the blade also create additional couplings among various DOFs of motion. Wang et al. [19] studied the vibration characteristics of a composite helicopter rotor blade to investigate the effect of moisture absorption on the vibration characteristics using finite element analysis (FEA). Tian et al. [20] conducted the failure analysis of a composite helicopter rotor blade using the VABS tool to explore the local-global design optimization and presented numerical examples to show the capability of the toolset for designing composite blades. In addition to the coupling among the DOFs of motion of the rotor blade generated from various styles of composite lay-ups, the curved tip of the blade offers extra coupling to influence the frequencies of vibration [21], which has the potential to be used as a morphing concept for regulating the vibration. Taymaz [22] presented an optimization study showing that the structural orientation of composite layers in the rotor blade could significantly affect its vibration characteristics and needed to be considered as an important factor for designing a passive vibration control system. All these studies justify that proper prediction of vibration is crucial for estimating the structural deformation of the rotor blade and utilizing that to design the relevant optimal vibration controller.

To unravel the core vibration dynamics in a triply coupled rotor system, a theoretically transparent yet computationally efficient framework is essential, providing a baseline for the development and validation of more advanced control strategies. Majority of the previous research activities are based on identification of the rotor blade model experimentally, or through finite element simulation having a single DOF of motion or considering the blade as a rigid one. However, a fundamental study based on classical mechanics and mathematical modeling for coupled vibration among multiple DOFs of motion of a flexible rotor blade is essential for comprehending the underlying physics and for attenuating the vibration. To identify the research needs, a mathematical model of vibration for an isolated Bo 105 helicopter rotor blade having coupled, three DOFs of motion are derived in this study. An analytical method is developed to estimate the sectional properties of the rotor blade. The frequencies of vibration are estimated from the modified Galerkin method (MGM) [23] and FEA and compared with experimental results. The quasi-steady-state lift, drag, and torsional moment are considered as principal aerodynamic loadings for both hovering and forward flights. The state-space model of the coupled equations of motion is derived by developing the orthogonality condition for triply coupled DOFs of motion. The IBC technique is implemented to reduce the three DOFs of motion rotor blade vibration. A Linear Quadratic Regulator (LQR) controller is designed based on the derived state-space model using MATLAB (MathWorks, 2024) control system toolbox. The controller is tuned to achieve optimal vibration reduction for all the DOFs of the rotor blade. This actuator-free analysis not only avoids hardware-specific constraints but also provides a benchmark for validating numerical or experimental implementations in future studies.

2. Materials and Methods

2.1. Modeling of the Composite Helicopter Rotor Blade

2.1.1. Physical model: The Bo 105 Helicopter Rotor Blade

The physical model of vibration estimation and control incorporates the Messerschmitt-Bölkow-Blohm Bo 105 helicopter rotor blade as the base model [24]. The Bo 105 helicopter rotor blade is flexible and has three, coupled DOFs of motion–the out-of-plane or the flapping deflection, $w$; the in-plane or lead-lag deflection, $v$; and the torsional displacement, $\theta$ of the rotor blade about the elastic axis. A sketch of such a rotor blade is shown in Figure 1(a) with three types of aerodynamic loadings associated with the three aforementioned DOFs of motion–the out-of-plane bending force per unit length, $f_z$ along $z$ axis; in-plane bending force per unit length, $f_y$ along $y$ axis; and pitching or torsional moment per unit length, $m_x$ about the elastic axis or $x$ axis. The distance between the centroidal axis (CA) and the elastic $x$ axis is expressed as $e$. The Bo 105 rotor blade is also considered to have a hingeless rotor system which indicates that the root of the rotor blade is fixed securely to the support [24,25]. A typical hingeless rotor system is described in Figure 1(b).

2.1.2. Cross-Sectional Analysis

The proposed cross-section of the rotor blade for free vibration analysis is based on the geometry of the Bo 105 helicopter rotor blade [24] section. Figure 2 describes a proposed Bo 105 rotor blade section composed of a thin, orthotropic fiberglass-epoxy [25], outer composite shell with the fiber volume fraction of 0.6 and a thick, inner composite core composed of an isotropic polymethacrylimide (PMI) foam called Rohacell [25] and an isotropic Honeycomb structure [25]. The profile of the cross section resembles the NACA 23012 airfoil [26]. The composite shell acts as a balanced laminate composed of four plies with a high strength to weight ratio compared to the inner core of the cross-section and can be considered as a sandwich beam. This cross-sectional configuration of the Bo 105 rotor blade is considered uniform throughout its length from blade root to the blade tip. The blade features a hingeless rotor system which can be treated as a cantilever blade having a fixed end at the root [24,25].

Properties of the materials in the composite section in Figure 2 are given in Table 1 and Table 2 [25] where, ${\rho }_{shell}$, ${\nu }_{12}$, $h_{shell}$, and $G_{12}$ indicate the density, major Poisson’s ratio, thickness, and in-plane shear modulus of the composite shell, respectively; $E_1$ and $E_2$ indicate the elastic modulus of the composite shell in the fiber direction and transverse direction, respectively; and $E_h$ and $E_r$ represent the elastic modulus of the isotropic Honeycomb and Rohacell cores, respectively.

The rotor blade can be considered as a long, thin plate so that the plane stress assumptions hold good. However, the curvature effect of the blade section contour must be considered in the analysis. The key part of the cross-sectional analysis is to find the profile of the reference surface shown in Figure 2, which is a closed curve parallel to the outermost or the innermost profile of the cross section and is located exactly in the middle of the outermost and the innermost profile of the composite shell. It can be treated as the offset curve either from the outermost or the innermost shell profile separated by a distance of $h_{shell}/2$, where $h_{shell}$ is the total thickness of the composite layers. The outermost composite shell profile in Figure 2 is generated from the NACA 23012 data (Airfoil Coordinates Database, UIUC). Both the generated outermost upper and lower profiles of the section are fitted into MATLAB (MathWorks, 2024) by 10^th degree polynomial curves to obtain the corresponding profile equations. Data points on the reference surface are generated from the slope of the outermost surfaces for enough number of points. The whole procedure of generating the coordinates of points on the reference surface from the outermost surface of the composite shell is demonstrated in Figure 3.

From Figure 3, the radial distance between the outermost surface and the reference surface is denoted by $d_r$, which also represents the length of the common normal to the respective tangents at two different points on the two curves. To see the variation of the slope of the outermost surface of the composite shell, the whole cross section is divided into four regions, depending on the directional change of the tangents. Starting from the leading edge of the airfoil, all the four regions show drastic changes of the slope angle, ${\alpha }_v$. Since $c_{h1}$ and $c_{v1}$ are already known from the NACA 23012 airfoil data (Airfoil Coordinates Database, UIUC), only remaining equations to calculate $\mathsf{\Delta }{c}_h$, $\mathsf{\Delta }{c}_v$, and $d_r$ are:

$$\mathsf{\Delta }{c}_h=d_r\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_v$$

(1)

$$\mathsf{\Delta }{c}_v=d_r\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_v$$

(2)

$$d_r=h_{shell}/2$$

(3)

where $\mathsf{\Delta }{c}_h$ is the distance between $c_{h1}$ and $c_{h2}$ along the $c_h$ axis and $\mathsf{\Delta }{c}_v$ is the distance between $c_{v1}$ and $c_{v2}$ along the $c_v$ axis. Once the coordinates of sufficient number of points on the reference surface are known, the curves (upper and lower reference curves) are centered and approximately fitted in MATLAB (MathWorks, 2024) by a 10^th degree polynomial and profile equations of the upper and lower reference surface of the composite shell are obtained.

The major cross-sectional properties required for the coupled, 3-DOFs of motion vibration analysis of the rotor blade are extensional stiffness, bending stiffness (flapping and lead-lag stiffness), and torsional stiffness. Vasiliev and Morozov [27] derived an analytical procedure to find the bending stiffness of the composite shell once the extensional stiffness per unit width of the composite shell, $B_{shell}$ is known. However, calculation of $B_{shell}$ in turn requires estimation of the equivalent extensional modulus of elasticity of the composite shell, $E_{shell}$. In this research, a basic principle that is considered is, $E_{shell}$ will be the same for any arbitrarily curved thin composite shell having the same number of layers, consisting of the same composite material, and having the same length, as described in Figure 4. Using this principle, $E_{shell}$ is estimated based on the classical lamination theory [28] from which $B_{shell}$ is calculated. Following the same procedure, equivalent shear modulus of the composite shell, $G_{shell}$ is also calculated [28].

The theory discussed above is applied to find the natural frequencies of free vibration of a thin-walled beam having a circular, composite section consisting of four layers with ply orientations 45$°$, $-$45$°$, $-$45$°$, and 45$°$ from the inner side to the outer side of the circular cross section [29]. Excellent agreement is found between the natural frequencies obtained from this theory and that given by Hodges et al. [29] which substantiates a good validation of the theory.

Figure 5 shows the schematic of the cross section of the rotor blade in Figure 2 with the presence of the CA and the neutral axis (NA). From Figure 5, the CA and the principal centroidal axis ${\overline{y}}_p$ are not exactly aligned. However, the separation angle is found approximately 1$°$ from SolidWorks (Dassault Systèmes, 2024), which is negligibly small to affect the sectional properties of the shell and the core. The bending stiffness of the whole section is to be calculated with respect to the NA which is not located at the geometric centroid of the section. This is because the cross section is composed of composite materials and the NA is shifted by a distance $d$from the centroid. Due to the small separation angle between the CA and NA, the NA can be assumed horizontal and almost parallel to the CA. The geometric parameters of the cross section ${\overline{d}}_{shell}$, ${\overline{d}}_{core}$, $A_{c_{shell}}$, $A_{c_{core}}$, ${\overline{I}}_{{y_p}_{shell}}$, ${\overline{I}}_{{y_p}_{core}}$, ${\gamma }_{shell}$, ${\gamma }_{core}$, and $e$ are determined by using SolidWorks (Dassault Systèmes, 2024), where $\overline{d}$, $A_c$, ${\overline{I}}_{y_p}$, and $\gamma$ are the distance of the CA from the $c_h$ axis, area of the entire cross section, area moment of inertia of the entire section about the principal centroidal ${\overline{y}}_p$ axis, and torsional constant of the section, respectively. The subscripts $shell$ and $core$ stand for composite shell and composite core of the cross section, respectively.

2.1.3. Extensional Stiffness Per Unit Width

From the lamination theory [28],

$$E_{shell}={\displaystyle \frac{A_{11}{A}_{22}-A_{12}^2}{h_{shell}{A}_{22}}}$$

(4)

where $A_{11}$, $A_{12}$, and $A_{22}$ are the elements of the extensional stiffness matrix of a composite laminate. Then $B_{shell}$ is evaluated as,

$$B_{shell}={\displaystyle \frac{E_{shell}{A}_{c_{shell}}}{s_{uref}+s_{lref}}}$$

(5)

where $s_{uref}$ and $s_{lref}$ are the arc lengths of the upper and the lower reference curves in Figure 5, respectively.

2.1.4. Out-of-Plane Bending (Flapping) and In-Plane Bending (Lead-Lag) Stiffness

Using the value of $B_{shell}$ and including the curvature effect of the cross-sectional profile, bending stiffness of the composite shell is estimated as:

$$E_{shell}{\overline{I}}_{{y_p}_{shell}}=B_{shell}\left[{\int }_{0.001}^{0.2651}{c}_{v_{uref}}^{2}d{s}_{uref}+{\int }_{0.001}^{0.2651}{c}_{v_{lref}}^{2}d{s}_{lref}\right]$$

(6)

where $c_{v_{uref}}$ and $c_{v_{lref}}$ are the ordinates of the points on the upper and the lower reference curves obtained from MATLAB (MathWorks, 2024) curve fitting. Equation (7) below provides an equivalent extensional elastic modulus of the composite core $E_{core}$ when the constituent materials have relatively close values of elastic moduli for isotropic materials.

$$E_{core}={\displaystyle \frac{E_r{A}_{c_r}+E_h{A}_{c_h}}{A_{c_r}+A_{c_h}}}={\displaystyle \frac{E_r{A}_{c_r}+E_h{A}_{c_h}}{A_{c_{core}}}}$$

(7)

where $A_{c_{core}}=A_{c_r}+A_{c_h}$ with the subscripts $r$ and $h$ standing for the Rohacell foam and the Honeycomb structure [25], respectively. Using linear strain assumption, location of the NA of the section is found by solving Equations (9) and (10) as follows:

$$E_{shell}{d}_{shell}{A}_{c_{shell}}+E_{core}{d}_{core}{A}_{c_{core}}=0$$

(8)

$${\Rightarrow E}_{shell}\left(c_{v_{upper}}-h_1-{\overline{d}}_{shell}\right)A_{c_{shell}}+E_{core}\left(c_{v_{upper}}-h_1-{\overline{d}}_{core}\right)A_{c_{core}}=0$$

(9)

$$h_1+h_2=h_{airfoil}$$

(10)

where $c_{v_{upper}}$ is the ordinate of the point on the outermost upper surface measured along the $c_v$ axis from the centroid of the whole section while $h_1$ and $h_2$ represent the distances of the topmost and the bottommost points from the NA of the combined section, respectively, as shown in Figure 5. The term $h_{airfoil}$ indicates the maximum thickness of the NACA 23012 airfoil [26] which is 12% of its total chord length. The out-of-plane bending or flapping stiffness is calculated from Equation (13) below by using Equations (11) and (12).

$${I_{NA}}_{shell}={\overline{I}}_{{y_p}_{shell}}+{A_c}_{shell}{d}_{shell}^2$$

(11)

$${I_{NA}}_{core}={\overline{I}}_{{y_p}_{core}}+{A_c}_{core}{d}_{core}^2$$

(12)

$$D_{bNA}=E_{shell}{I}_{N{A}_{shell}}+E_{core}{I_{NA}}_{core}$$

(13)

where $D_{bNA}$ and $I_{NA}$ indicate the bending stiffness with respect to the NA of the whole cross section and the area moment of inertia of the shell/core areas with respect to the respective NA, respectively. Due to the small angle assumption, the axis perpendicular to the NA of the whole composite section passes through the centroid almost vertically and can be defined approximately with respect to the ${\overline{z}}_p$ axis. Therefore, the corresponding bending stiffness expression can be written as,

$$D_{b{\overline{z}}_p}\cong E_{shell}{\overline{I}}_{{z_p}_{shell}}+E_{core}{\overline{I}}_{{z_p}_{core}}$$

(14)

where $D_{b{\overline{z}}_p}$ and ${\overline{I}}_{z_p}$ are the bending stiffness with respect to the centroidal ${\overline{z}}_p$ axis of the whole cross section and the area moment of inertia of either the composite shell or the core with respect to the respective ${\overline{z}}_p$ axis, respectively. If $d$ and $\overline{d}$ in Figure 5 are very small for an asymmetric cross section, then the NA and CA become very close and $c_h$ or the chord axis can be considered as the axis of symmetry. A small value of $d$ causes the NA to coincide with the CA and at the same time, a small separation angle between the CA and the ${\overline{y}}_p$ axis in Figure 5 makes the cross section approximately symmetric with respect to the ${\overline{y}}_p$ axis, i.e., the $c_h$ axis of Figure 5. $y$ and $z$axes are termed as the global orthogonal axes passing through the shear center of the cross section. If the section is approximately symmetric with respect to the ${\overline{y}}_p$ axis and maintains a small angle with the global y axis, the out-of-plane bending (flapping) stiffness with respect to the global y axis passing through the shear center can be approximated as:

$$D_{by}\cong D_{b{\overline{y}}_p}=D_{bNA}$$

(15)

Similarly, the in-plane bending stiffness with respect to global $z$ axis through the shear center for small value of $e$ can be expressed as:

$$D_{bz}\cong D_{b{\overline{z}}_p}$$

(16)

2.1.5. Torsional Stiffness

The individual shear modulus $G_{shell}$ and $G_{core}$ for the shell and the core, respectively, can be approximated as below [28],

$$G_{shell}={\displaystyle \frac{A_{66}}{h_{shell}}}$$

(17)

$$G_{core}={\displaystyle \frac{E_{core}}{2\left(1+\nu \right)}}$$

(18)

where $A_{66}$ is the element of the extensional stiffness matrix of a composite laminate and $\nu$ is the Poisson’s ratio of isotropic material. Equation (18) is an approximation since both the Rohacell and the Honeycomb structures [25] are considered isotropic and use the same $\nu$. The shear center and centroid of the cross section are relatively close to each other and the total torsional stiffness $D_t$ is written as,

$$D_t\cong G_{shell}{\gamma }_{comp}+G_{core}{\gamma }_{core}$$

(19)

2.1.6. Mass Per Unit Length

The mass per unit length, $m$ of the rotor blade having the proposed composite section in Figure 2 is the sum of the masses per unit length of the individual constituent beams with respective cross sections making up the whole composite section and is expressed as,

$$m={\rho }_{shell}{A}_{c_{shell}}+{\rho }_r{A}_{c_r}+{\rho }_h{A}_{c_h}$$

(20)

where $\rho$ is the material density with appropriate meaning of the already defined subscripts.

2.2. Modeling of Free Vibration for the Composite Helicopter Rotor Blade

2.2.1. Governing Equations of Motion

The Bo 105 helicopter rotor blade model is ideal to be considered as the Euler-Bernoulli beam as the blade can be approximated as a long, thin plate having a solid section for which the effect of sectional warping is ignored due to its relatively minor role [30]. The triply coupled governing equations of motion for the Bo 105 helicopter rotor blade vibration are given by Equations (21)$-$(23) [24]:

$$\begin{gathered} {[D}_{by}(x)w″(x,t)+D_{bzy}v″(x,t)]″-{\left[T{\left(x\right)w}^{\prime }\left(x,t\right)\right]}^{\prime }-{\left[{\mathsf{\Omega }}^{2}me\left(x+e_1\right)\theta \left(x,t\right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha \left(x\right)\right]}^{\prime } \\ +m[\ddot{w}(x,t)+e\ddot{\theta }(x,t)\mathrm{c}\mathrm{o}\mathrm{s}\alpha (x)]=f_z(x,t)+{\left[{\mathsf{\Omega }}^{2}mex\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)\right]}^{\prime } \end{gathered}$$

(21)

$$\begin{gathered} {\left[D_{bz}{\left(x\right)v}^{″}\left(x,t\right)+D_{bzy}{w}^{″}\left(x,t\right)\right]}^{″}-{\left[T{\left(x\right)v}^{\prime }\left(x,t\right)\right]}^{\prime }+{\left[{\mathsf{\Omega }}^{2}me\left(x+e_1\right)\theta \left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)\right]}^{\prime } \\ +{\mathsf{\Omega }}^{2}me\theta \left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)+m\left[\ddot{v}\left(x,t\right)-e\ddot{\theta }\left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)\right]-{\mathsf{\Omega }}^{2}mv\left(x,t\right) \\ =f_y\left(x,t\right)+{[\mathsf{\Omega }}^{2}me\left(x+e_1\right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha (x)]\prime +{\mathsf{\Omega }}^{2}me\mathrm{c}\mathrm{o}\mathrm{s}\alpha (x) \end{gathered}$$

(22)

$$\begin{gathered} -{\left[D_t\left(x\right){\theta }^{\prime }\left(x,t\right)\right]}^{\prime }-{\mathsf{\Omega }}^{2}me\left(x+e_1\right)\left[v^{\prime }\left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)-w^{\prime }\left(x,t\right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha \left(x\right)\right]+{\mathsf{\Omega }}^{2}mev\left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right) \\ +{\mathsf{\Omega }}^{2}m\cos 2\alpha \left(x\right)\left({\kappa }_{m2}^2-{\kappa }_{m1}^2\right)\theta \left(x,t\right)+m{\kappa }_m^2\ddot{\theta }\left(x,t\right)-me\left[\ddot{v}\left(x,t\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)-\ddot{w}\left(x,t\right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha \left(x\right)\right] \\ =m_x\left(x,t\right)-{\mathsf{\Omega }}^{2}m\left[\left({\kappa }_{m2}^2-{\kappa }_{m1}^2\right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \left(x\right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha \left(x\right)\right] \end{gathered}$$

(23)

with

$$T\left(x\right)={\int }_x^l{\mathsf{\Omega }}^{2}m\left(x+e_1\right)dx$$

(24)

where ($\prime$) denotes differentiation with respect to $x$, and ($̇$) indicates differentiation with respect to time, $t$; $D_{by}$ and $D_{bz}$ are the bending stiffnesses with respect to y and z axes, respectively; $D_{bzy}$ is bending stiffness associated with respect to $yz$ axis system; $D_t$ is torsional rigidity with respect to the x axis; w, v, and θ are the flapping, lead-lag, and torsional deflections, respectively; $\alpha$is the total twist angle of the helicopter rotor blade section at any location $x$; $l$ is the blade length; $e_1$ is the offset of the root of the blade from the rotating axis; $f_z$, $f_y$, and $m_x$ are the aerodynamic lift, drag, and pitching moment per unit length; ${\kappa }_m$ is the polar mass radius of gyration about the elastic axis; ${\kappa }_{m1}$, ${\kappa }_{m2}$ are the mass radii of gyration about the NA and the axis normal to chord through the shear center, respectively; $m$ is the mass per unit length; $T$ is the centrifugal tension due to the rotation of the blade; and $\mathsf{\Omega }$ is the angular velocity of the rotor blade. The rotor blades of the Bo 105 helicopter accommodate aerodynamic forces by flexing and reducing the drag forces at the root. The other properties of the Bo 105 rotor system are listed in Table 3 [24].

2.2.2. Boundary Conditions

Equations (21)–(23) represent an initial boundary value problem and involve ten boundary conditions (BCs) [24] for the most general case which are given as:

$$\mathrm{at} x=0: w\left(0,t\right)=0; w^{\prime }\left(0,t\right)=0; v\left(0,t\right)=0; { v}^{\prime }\left(0,t\right)=0; \theta \left(0,t\right)=0$$

(25)

$$\mathrm{at} x=l: S_z\left(l,t\right)=M_y\left(l,t\right)=S_y\left(l,t\right)=M_z\left(l,t\right)=Q\left(l,t\right)=0$$

(26)

where

$$S_z=-M_y^{\prime }+T{w}^{\prime }+{\mathsf{\Omega }}^{2}me(x+e_1)(sin\alpha +\theta cos\alpha )$$

(27)

$$S_y=-M_z^{\prime }+T{v}^{\prime }+{\mathsf{\Omega }}^{2}me\left(x+e_1\right)\left(cos\alpha -\theta sin\alpha \right)$$

(28)

$$M_y=D_{by}{w}^{″}+D_{bzy}{v}^{″}$$

(29)

$$M_z=D_{bz}v″+D_{bzy}{w}^{″}$$

(30)

$$Q_t=D_t{\theta }^{\prime }$$

(31)

with $S_y$ and $S_z$ are the shear forces along $y$ and $z$ axis, respectively, $M_y$ and $M_z$ are the bending moments about $y$ and $z$ axis, respectively, and $Q_t$ is the torque about the $x$ axis. Equations (21)–(23) and the BCs in Equations (25) and (26) contain the term $\alpha$, which does not affect the coupled vibration frequencies significantly. However, when calculating the forcing functions, $\alpha$ must be included appropriately since the aerodynamic loadings vary significantly even with a small change in $\alpha$.

2.2.3. Natural Frequencies of Free Vibration from the Modified Galerkin Method

Closed form solutions are not available for Equations (21)–(23) and approximate techniques are required for solving them in terms of natural frequencies. In this research, the MGM [23] is used to solve Equations (21)–(23) with some approximations in which the displacements $w\left(x\right)$, $v\left(x\right)$, and $\theta \left(x\right)$ are expressed in finite series of known functions $W_n\left(x\right)$, $V_n\left(x\right)$, and ${\mathsf{\Theta }}_n\left(x\right)$, respectively, which are the uncoupled out-of-plane bending, in-plane bending, and torsional vibration mode shapes of the rotor blade, respectively [23]. Involving these uncoupled mode shapes in relation to Equations (21)−(23), the following expressions are obtained:

$${\int }_0^l{W}_n{L}_w\left(w,v,\theta \right)dx+{\left[M_y{W}_n^{\prime }\right]}_0^l+{\left[{S_{z}W}_n\right]}_0^l=0$$

(32)

$${\int }_0^l{V}_n{L}_v\left(w,v,\theta \right)dx+{\left[M_z{V}_n^{\prime }\right]}_0^l+{\left[{S_{y}V}_n\right]}_0^l=0$$

(33)

$${\int }_0^l{\mathsf{\Theta }}_n{L}_{\theta }\left(w,v,\theta \right)dx+{\left[Q_t{\mathsf{\Theta }}_n\right]}_0^l=0$$

(34)

where $L_w(w,v,\theta )$, $L_v(w,v,\theta )$, and $L_{\theta }(w,v,\theta$) denote Equations (21), (22), and (23), respectively, with all the right hand terms transferred to the left hand side. For static or externally applied loads, the set of Equations (32)–(34) are nonhomogeneous; however, for a characteristic value problem, they will be homogeneous. Equations (32)–(34) are finally solved by introducing the frequency term $\omega$ into the functions $L_w(w,v,\theta )$, $L_v(w,v,\theta )$, and $L_{\theta }(w,v,\theta$) and the procedure is described in detail in the study conducted by Sarker [23] .

2.2.4. Natural Frequencies of Free Vibration from Finite Element Analysis

To verify the frequencies obtained from the MGM [23], a computational finite element model of free vibration for the composite rotor blade is developed in the commercial finite element package, Abaqus (Dassault Systèmes, 2024). The finite element model of the rotor blade incorporates properties of the blade section in Figure 2 obtained by the method described in Section 2.1.2, Section 2.1.3, Section 2.1.4, Section 2.1.5 and Section 2.1.6. Figure 6(a) and Figure 6(b) describe the finite element domain and magnified view of part of the meshed finite element domain, respectively, developed in Abaqus (Dassault Systèmes, 2024). For meshing the finite element domain, 8-node, liner brick, C3D8R elements are used, and 252,000 elements are considered for convergence of results. The finite element model of free vibration is solved for natural frequencies for both rotating and nonrotating cases.

2.3. Orthogonality Condition for Triply Couple Vibration and the State Space Model

2.3.1. Generalized Force and Moment for Coupled Vibration of the Rotor Blade

Since the rotor blade is subjected to triply coupled vibratory loads, the conventional, uncoupled orthogonality relationship will not be applicable, and the corresponding orthogonality relationship needs to be developed. Considering solutions of the governing equations (21)–(23) harmonic in nature, the displacements have the following forms:

$$w\left(x,t\right)=W_n\left(x\right)e^{i{\omega }_{n}t}$$

(35)

$$v\left(x,t\right)=V_n\left(x\right)e^{i{\omega }_{n}t}$$

(36)

$$\theta \left(x,t\right)={\mathsf{\Theta }}_n\left(x\right)e^{i{\omega }_{n}t}$$

(37)

where ${\omega }_n$ is the $n$^th natural frequency of vibration in rad/s. By plugging Equations (35)$-$(37) into Equations (21)$-$(23) for free vibration analysis with $f_z{=f}_y=m_x=$ 0, and for the nonrotating case with $\mathsf{\Omega }=$ 0 having no blade twist, the following expressions are obtained:

$$D_{by}{W}_n^{⁗}-m{\omega }_n^2{W}_n-me{\omega }_n^2{\mathsf{\Theta }}_n=0$$

(38)

$$D_{bz}{V}_n^{⁗}-m{\omega }_n^2{V}_n=0$$

(39)

$${-D}_t{\mathsf{\Theta }}_n^{″}-m{\omega }_n^2{\kappa }_m^2{\mathsf{\Theta }}_n-me{\omega }_n^2{W}_n=0$$

(40)

Multiplying Equation (38) with $W_m$ and integrating from $0$ to $l$,

$${\int }_0^l\left(D_{by}{W}_n^{⁗}{W}_m-m{\omega }_n^2{W}_n{W}_m-me{\omega }_n^2{\mathsf{\Theta }}_n{W}_m\right)dx=0$$

(41)

Interchanging $m$ and $n$ in Equation (41),

$${\int }_0^l\left(D_{by}{W}_m^{⁗}{W}_n-m{\omega }_m^2{W}_m{W}_n-me{\omega }_m^2{\mathsf{\Theta }}_m{W}_n\right)dx=0$$

(42)

Subtracting Equation (42) from Equation (41) and applying integration by parts using the BCs in Equation (25) and (26) with $\mathsf{\Omega }=\alpha =0$,

$$\left(m{\omega }_m^2-m{\omega }_n^2\right){\int }_0^l{W}_m{W}_{n}dx+m{\omega }_m^2{\int }_0^{l}e{W}_n{\mathsf{\Theta }}_{m}dx-m{\omega }_n^2{\int }_0^{l}e{W}_m{\mathsf{\Theta }}_{n}dx=0$$

(43)

Similarly, multiplying Equation (39) by $V_m$ and Equation (40) by ${\mathsf{\Theta }}_m$ and following the above procedure,

$$\left(m{\omega }_m^2-m{\omega }_n^2\right){\int }_0^l{V}_m{V}_{n}dx=0$$

(44)

$$\left(m{\omega }_m^2-m{\omega }_n^2\right){\int }_0^l{\kappa }_m^2{\mathsf{\Theta }}_m{\mathsf{\Theta }}_{n}dx+m{\omega }_m^2{\int }_0^{l}e{W}_m{\mathsf{\Theta }}_{n}dx-m{\omega }_n^2{\int }_0^{l}e{W}_n{\mathsf{\Theta }}_{m}dx=0$$

(45)

Adding Equations (43)–(45) together,

$$\left(m{\omega }_m^2-m{\omega }_n^2\right){\int }_0^l{[W}_m{W}_n+V_m{V}_n+e\left(W_m{\mathsf{\Theta }}_n+W_n{\mathsf{\Theta }}_m\right)+{\kappa }_m^2{\mathsf{\Theta }}_m{\mathsf{\Theta }}_n]dx=0$$

(46)

Since for two different modes $m$ and $n$, ${\omega }_m$ and ${\omega }_n$ will be different, therefore, the term $\left(m{\omega }_m^2-m{\omega }_n^2\right)\ne 0$, and the orthogonality condition for the triply coupled beam is written as,

$${\int }_0^l{[W}_m{W}_n+V_m{V}_n+e\left(W_m{\mathsf{\Theta }}_n+W_n{\mathsf{\Theta }}_m\right)+{\kappa }_m^2{\mathsf{\Theta }}_m{\mathsf{\Theta }}_n]dx={\delta }_{mn}$$

(47)

where ${\delta }_{mn}$ is the Kronecker delta defined as,

$${\delta }_{mn}=1;\hspace{1em}m=n$$

(48)

$${\delta }_{mn}=0;\hspace{1em}m\ne n$$

(49)

The orthogonality condition in Equation (47) will be used for estimation of the time-varying deflections of the rotor blade and can be expressed as the following form of series solutions:

$$w\left(x,t\right)=\sum _{n=1}^{\infty }{W}_n\left(x\right)q_n\left(t\right)$$

(50)

$$v\left(x,t\right)=\sum _{n=1}^{\infty }{V}_n\left(x\right)q_n\left(t\right)$$

(51)

$$\theta \left(x,t\right)=\sum _{n=1}^{\infty }{\mathsf{\Theta }}_n\left(x\right)q_n\left(t\right)$$

(52)

where $q_n$ is the generalized time coordinate for the $n$^th mode. Now, for the forced vibration of a nonrotating beam, Equations (21)–(23) with $\alpha =$ 0 reduce to,

$$D_{by}{w}^{⁗}+m\ddot{w}+me\ddot{\theta }=f_z$$

(53)

$$D_{bz}{v}^{⁗}+m\ddot{v}=f_y$$

(54)

$$-D_t{\theta }^{″}+m{\kappa }_m^2\ddot{\theta }+me\ddot{w}=m_x$$

(55)

Substituting Equations (50)–(52) into Equations (53)–(55),

$$\sum _{n=1}^{\infty }{D_{by}W}_n⁗q_n+\sum _{n=1}^{\infty }{mW}_n\ddot{q_n}+\sum _{n=1}^{\infty }me{\mathsf{\Theta }}_n\ddot{q_n}=f_z$$

(56)

$$\sum _{n=1}^{\infty }{D_{bz}V}_n⁗q_n+\sum _{n=1}^{\infty }m{V}_n\ddot{q_n}=f_y$$

(57)

$$\sum _{n=1}^{\infty }-D_t{\mathsf{\Theta }}_n″q_n+\sum _{n=1}^{\infty }m{\kappa }_m^2{\mathsf{\Theta }}_n\ddot{q_n}+\sum _{n=1}^{\infty }me{W}_n\ddot{q_n}=m_x$$

(58)

Multiplying Equations (56), (57), and (58) by $W_m$, $V_m$, and ${\mathsf{\Theta }}_m$, respectively, using the relationships from Equations (38)–(40), and adding them together,

$$\begin{gathered} \sum _{n=1}^{\infty }{\left[\right(m\omega }_n^2{W}_n+me{\omega }_n^2{\mathsf{\Theta }}_n)W_m{q}_n+{mW}_n{W}_m\ddot{q_n}+m{e\mathsf{\Theta }}_n{W}_m\ddot{q_n}+{{m\omega }_n^2{V}_n{V}_{m}q}_n+m{V}_n{V}_m\ddot{q_n} \\ +\left(m{\omega }_n^2{\kappa }_m^2{\mathsf{\Theta }}_n+me{\omega }_n^2{W}_n{){\mathsf{\Theta }}_{m}q}_n+m{\kappa }_m^2{\mathsf{\Theta }}_n{\mathsf{\Theta }}_m\ddot{q_n}+me{W}_n{\mathsf{\Theta }}_m\ddot{q_n}\right]=f_z{W}_m+f_y{V}_m+m_x{\mathsf{\Theta }}_m \end{gathered}$$

(59)

After rearranging the terms and integrating from $0$ to $l$,

$$\begin{gathered} \sum _{n=0}^{\infty }\{m\ddot{q_n}{\int }_0^l\left[W_n{W}_m+V_n{V}_m+e\left(W_m{\mathsf{\Theta }}_n+W_n{\mathsf{\Theta }}_m\right)+{\kappa }_m^2{\mathsf{\Theta }}_n{\mathsf{\Theta }}_m\right]dx \\ +m{\omega }_n^2{q}_n{\int }_0^l\left[W_n{W}_m+V_n{V}_m+e\left(W_m{\mathsf{\Theta }}_n+W_n{\mathsf{\Theta }}_m\right)+{\kappa }_m^2{\mathsf{\Theta }}_n{\mathsf{\Theta }}_m\right]dx\} \\ ={\int }_0^l(f_z{W}_m+f_y{V}_m+m_x{\mathsf{\Theta }}_m)dx \end{gathered}$$

(60)

Using the orthogonality condition from Equation (47) and Equation (48) for $n=m$,

$$\ddot{q_n}+{\omega }_n^2{q}_n={\displaystyle \frac{1}{m}}{\int }_0^l(f_z{W}_n+f_y{V}_n+m_x{\mathsf{\Theta }}_n)dx=F_z+F_y+M_x$$

(61)

with

$$F_z=F_z\left(x,t\right)={\displaystyle \frac{1}{m}}{\int }_0^l{f}_z{\left(x,t\right)W}_n\left(x\right)dx$$

(62)

$$F_y=F_y\left(x,t\right)={\displaystyle \frac{1}{m}}{\int }_0^l{f}_y{\left(x,t\right)V}_n\left(x\right)dx$$

(63)

$$M_x=M_x\left(x,t\right)={\displaystyle \frac{1}{m}}{\int }_0^l{m}_x\left(x,t\right){\mathsf{\Theta }}_n\left(x\right)dx$$

(64)

where $F_z\left(t\right)$ and $F_y\left(t\right)$ are the $n$^th generalized out-of-plane bending/flapping and in-plane bending/lead-lag forces, respectively; and $M_x\left(t\right)$ is the $n$^th generalized moment for torsional motion.

With damping considered, the form of Equation (61) takes the following form,

$$\ddot{q_n}+2 {\zeta }_n{\omega }_n\dot{q_n}+{\omega }_n^2{q}_n=F_z+F_y+M_x$$

(65)

where ${\zeta }_n$ represents the damping ratio corresponding to the n^th vibration mode which governs the rate of energy dissipation in the system. For the rotating blade, the natural frequency ${\omega }_n$ in Equation (65) needs to be replaced by that of the rotating case, ${\omega }_{nr}$, and the above equation becomes,

$$\ddot{q_n}+2 {\zeta }_n{\omega }_{nr}\dot{q_n}+{\omega }_{nr}^2{q}_n=F_z+F_y+M_x=F$$

(66)

The damping ratios for the flapping, lead–lag, and torsional modes are selected based on their respective aerodynamic and structural characteristics, with higher damping typically assigned to the flapping mode, e.g., ${\zeta }_f\approx 0.20$ owing to stronger aerodynamic coupling effects, whereas smaller values are used for the lead–lag and torsional modes, e.g., ${\zeta }_l\approx 0.02$ and ${\zeta }_t\approx 0.02$, respectively [31].

2.3.2. The State-Space Model

To design the controller for a multiple input-multiple output (MIMO) system, a state-space model is required to analyze the response in the time domain. In addition, the development of the state-space model requires coupled natural frequencies of vibration for the three DOFs of motion. From Equation (66), the state-space model of the rotor blade can be derived as follows by using the fundamental natural frequencies of the 3 DOFs of motion for $n=$1:

$$\ddot{q_{1f}}=-2 {\zeta }_{1f}{\omega }_{1rf}\dot{q_{1f}}-{\omega }_{1rf}^2{q}_{1f}+F$$

(67)

$$\ddot{q_{1l}}=-2 {\zeta }_{1l}{\omega }_{1rl}\dot{q_{2l}}-{\omega }_{1rl}^2{q}_{1l}+F$$

(68)

$$\ddot{q_{1t}}=-2 {\zeta }_{1t}{\omega }_{1rt}\dot{q_{3t}}-{\omega }_{1rt}^2{q}_{1t}+F$$

(69)

where ${\omega }_{1f}$, ${\omega }_{2l}$, and ${\omega }_{3t}$ are the first flapping, lead-lag, and torsion governed natural frequencies, respectively. Let,

$$x_1=q_{1f}\left(t\right)$$

(70)

$$x_2={\dot{q}}_{1f}\left(t\right)={\dot{x}}_1$$

(71)

$$x_3=q_{1l}\left(t\right)$$

(72)

$$x_4={\dot{q}}_{1l}\left(t\right)={\dot{x}}_3$$

(73)

$$x_5=q_{1t}\left(t\right)$$

(74)

$$x_6={\dot{q}}_{1t}\left(t\right)={\dot{x}}_5$$

(75)

where, the states of the model for the first mode of triply coupled vibration are expressed as,

$x_1$= tip displacement for flapping motion.

$x_2={\dot{x}}_1$= derivative of the tip displacement for flapping motion

$x_3$ = tip displacement for lead-lag motion

$x_4$= ${\dot{x}}_3$= derivative of tip displacement for lead-lag motion

$x_5$ = tip rotational displacement for torsional motion

$x_6$=${\dot{x}}_5$= derivative of the tip rotational displacement due to torsional motion

Following this, the state-space model is derived as,

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\\ {\dot{x}}_4\\ {\dot{x}}_5\\ {\dot{x}}_6\end{array}\right]=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ -{\omega }_{1rf}^2& -2{\zeta }_{1f}{\omega }_{1rf}& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& -{\omega }_{1rl}^2& -2{\zeta }_{1l}{\omega }_{1rl}& 0& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& -{\omega }_{1rt}^2& -2{\zeta }_{1t}{\omega }_{1rt}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\end{array}\right]+\left[\begin{array}{ccc}0& 0& 0\\ F& 0& 0\\ 0& 0& 0\\ 0& F& 0\\ 0& 0& 0\\ 0& 0& F\end{array}\right]\mathrm{u}$$

(76)

where $\mathrm{u}$ is the control input. In this analysis, the desired outputs are flapping, lead-lag, and torsional deflections at the blade tip. Therefore, the output matrix y assumes the following form,

$$\mathrm{y}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\end{array}\right]$$

(77)

2.4. Aerodynamic Force and Moment

Helicopter rotor blades experience different types of loading during hovering and forward flights. In this analysis, the steady-state lift, drag, and pitching moment are considered for hovering and forward flights as the primary loads on the rotor blade. The helicopter is treated as a lifting machine with low disk loading [32] for which the aerodynamic strip theory (AST) can be applied to estimate $f_z$, $f_y$, and $m_x$ to be used in Equation (61). The AST considers the individual strips of a rotor blade along with the blade geometry, elemental forces, and torque for the analysis. The overall lift $L$, drag $D_d$, and pitching moment $M_{ac}$ are obtained by integration along the entire blade over a rotor revolution. The formulas to calculate the aerodynamic force and moment for hovering and forward flights are approximated as [24]:

$$\mathrm{Lift} \mathrm{force}\text{:} L\left(x,t\right)={\displaystyle \frac{1}{2}}{\rho }_a{u}_t^2\left(x,t\right)c{c}_l\left(x\right)$$

(78)

$$\mathrm{Drag} \mathrm{force}\text{:} D_d\left(x,t\right)={\displaystyle \frac{1}{2}}{\rho }_a{u}_t^2\left(x,t\right)c{c}_d\left(x\right)$$

(79)

$$\mathrm{Pitching} \mathrm{moment}\text{:} M_{ac}\left(x,t\right)={\displaystyle \frac{1}{2}}{\rho }_a{u}_t^2\left(x,t\right)c^2{c}_m\left(x\right)$$

(80)

where ‘ac’ as the subscript of $M$ in Equation (80) stands for aerodynamic center meaning that the aerodynamic moment acts about the aerodynamic center of the rotor blade; $u_t$ is the tangential velocity of the rotor blade parallel to the rotor disk plane; ${\rho }_a$ is the density of air; $c_l$, $c_d$, and $c_m$ are the lift, drag, and pitching moment coefficients, respectively; and $c$ is the chord length of the blade section. It is to note that the blade twist angle $\alpha$ in Equations (21)–(23) is the summation of the inflow angle ${\alpha }_{in}$ and the blade angle of attack (AoA) ${\alpha }_a$. In addition to varying locations along the rotor blade length, the aerodynamic forces and moment in helicopter dynamics are also functions of the azimuth angle, $\psi$ where $\psi =\mathsf{\Omega }t$. For the forward flight, $u_t$ can be approximated as,

$$u_t=\mathsf{\Omega }x+V_{\infty }\mathrm{s}\mathrm{i}\mathrm{n}\psi$$

(81)

where $V_{\infty }$ is the forward speed of the helicopter. $c_l$, $c_d$, and $c_m$ mentioned in Equations (78)–(80) can be expressed as,

$$c_l=a{\alpha }_a$$

(82)

$$c_d=d_0+d_1{\alpha }_a+d_2{\alpha }_a^2$$

(83)

$$c_m=m_0+m_1{\alpha }_a$$

(84)

where $a$ is the lift curve slope, and $d_0$, $d_1$, $d_2$, $m_0$, and $m_1$ are empirically derived coefficients for the NACA 23012 airfoil [26]. For the small-angle assumption, $m_1$ can be neglected, and the expression for ${\alpha }_a$ can be found in the literature [24]. The values of the other parameters used in this analysis include $a=$ 5.7/rad, ${ d}_0=0$.0087, $d_1=-0.012$, $d_2=$ 0.4, and $c_m=-$0.008 [33]. For small ${\alpha }_{in}$, the aerodynamic loading $f_z$, $f_y$, and $m_x$ used in Equations (21)–(23) can be estimated approximately as [24,34]:

$$f_z\left(x,t\right)\cong L\left(x,t\right)$$

(85)

$$f_y\left(x,t\right)\cong D_d\left(x,t\right)$$

(86)

$$m_x(x,t)\cong M_{ac}\left(x,t\right)+L\left(x,t\right)e$$

(87)

For periodic vibration, the magnitudes of the vibratory loads on the rotor blade can generally be modeled as fractional components of $L(x,t)$, $D_d(x,t)$, and $m_x(x,t)$ with harmonic nature (typically 0.5–5% of the aerodynamic loadings), and are applied as the forcing functions along with the appropriate excitation frequency following the implementations of Sarker [23] and Uddin et al. [31].

2.5. Controller Design

2.5.1. Optimal Vibration Control Framework

The primary objective of the controller is to attenuate the vibratory loads transmitted through the rotor system by generating appropriate input signals that counteract structural responses. Optimal control strategies are widely employed in modern rotorcraft vibration control, wherein the control inputs are selected to minimize a quadratic cost function that balances vibration suppression against the control effort. Such controllers are particularly attractive because of their robustness to unpredictable aerodynamic excitations and adaptability to variations in system states. Unlike frequency-domain approaches, modern optimal control theory is formulated in the time domain, assuming linear time-invariant dynamics, which makes it well suited for real-time vibration attenuation problems.

In this study, the rotor blade is modeled as a multi-input multi-output (MIMO) system with three coupled DOFs of motion [37]. The governing equations are expressed in the standard state-space form as,

$$\dot{\mathrm{x}}\left(t\right)=\mathrm{A}\mathrm{x}\left(t\right)+\mathrm{B}\mathrm{u}\left(t\right)$$

(88)

$$\mathrm{y}\left(t\right)=\mathrm{C}\mathrm{x}\left(t\right)+\mathrm{D}\mathrm{u}\left(t\right)$$

(89)

where $\mathrm{x}\left(t\right)$ is the state vector containing the generalized displacements and velocities of the coupled modes, $\mathrm{u}$ is the control input vector, and $\mathrm{y}$ is the output vector representing the measurable responses. The time-invariant matrices $\mathrm{A}$, $\mathrm{B}$, $\mathrm{C}$, and $\mathrm{D}$ are obtained from coupled equations of motion. This formulation provides a systematic framework for implementing an optimal LQR to suppress vibratory loads across all three modes of the rotor blade.

2.5.2. Linear Quadratic Regulator

The LQR, formulated on state-feedback optimization principles, continues to serve as a fundamental and well-established benchmark in helicopter vibration control research owing to its analytical transparency and optimal performance formulation [35,36,37]. This is an alternative approach to the direct pole placement technique, in which the best pole locations are implicitly chosen based on the Linear Quadratic algorithm. Using the corresponding derived state-space model, the optimal control law is implemented to determine the best feedback controller gain matrix $\mathrm{K}$ that reduces the vibration of the considered rotor blade. Figure 7 shows a block diagram of the LQR algorithm, where $\mathrm{r}$ is the reference signal.

For the regulator, the value of reference signal $\mathrm{r}$ is zero. Therefore, the state feedback control law is defined as:

$$\mathrm{u}\left(t\right)=-\mathrm{K}\mathrm{x}\left(t\right)$$

(90)

A quadratic cost function $J$ is defined to be minimized by optimizing the control effort and errors. Control errors are described by the squared values of the state variables, and control efforts are described by the squared values of the control input. To minimize the control errors, more control effort is required, and a reduction in the control effort increases the control error. By introducing the relative cost with two parameter matrices, $\mathrm{Q}$ and $\mathrm{R}$, these two contradictory objectives are quantified as follows:

$$J={\displaystyle \frac{1}{2}}{\int }_0^{\infty }\left[{\mathrm{x}}^{\mathrm{T}}\left(t\right)\mathrm{Q}\mathrm{x}\left(t\right)+{\mathrm{u}}^{\mathrm{T}}\left(t\right)\mathrm{R}\mathrm{u}\left(t\right)\right]dt$$

(91)

where the superscript $\mathrm{T}$ indicates the transpose of x. After determining the value of $\mathrm{R},$ the controller gain matrix $\mathrm{K}$ is expressed as follows:

$$\mathrm{K}={\mathrm{R}}^{-1}{\mathrm{B}}^{\mathrm{T}}\mathrm{P}$$

(92)

where $\mathrm{P}$ is the positive semi-definite unique matrix defined by the algebraic Riccati equation,

$$\mathrm{P}\mathrm{A}+{\mathrm{A}}^{\mathrm{T}}\mathrm{P}+{\mathrm{C}}^{\mathrm{T}}\mathrm{C}-\mathrm{P}\mathrm{B}{\mathrm{R}}^{-1}{\mathrm{B}}^{\mathrm{T}}\mathrm{P}=0$$

(93)

with $\mathrm{C}$ as the time-invariant matrix. For realizing the LQR controller, several assumptions are necessary, such as:

Matrix $\mathrm{A}$and $\mathrm{B}$are stable and detectable, and

(ii) The solution for the matrix $\mathrm{P}$ is always symmetric.

In the context of LQR design, the weighting matrices $\mathrm{Q}$ and $\mathrm{R}$ in the quadratic cost function directly influence the trade-off between vibration suppression and control effort. A higher value of $\mathrm{Q}$ at a constant R place more emphasis on minimizing state deviations, such as flapping, lead-lag, and torsional displacements, which results in more aggressive vibration reduction at the expense of increased control input demand. Conversely, a lower $\mathrm{Q}$ reduces the priority of state accuracy, potentially allowing larger residual vibrations but conserving the control effort.

3. Model Validation

3.1. Experimental Validation of Natural Frequencies

The state-space model regarding the coupled vibration is validated by comparing the natural frequency of the rotor blade obtained from this study with that obtained from experimental results reported by Jacklin et al. [16], who investigated the pitch-harmonic response of a full-scale IBC system in two major wind-tunnel campaigns at NASA Ames Research Center. Table 4 presents the fundamental natural frequencies of vibration of the triply coupled rotor blade for the rotating case, obtained from the MGM [23] and from the experimental study. In Table 4, $f_{1f}$, $f_{1l}$, and $f_{1t}$ indicate the fundamental frequencies In Hz where subscript 1 indicates the first mode and subscripts f, l, and t denote flapping, lead-lag, and torsional vibration modes, respectively. The predicted frequencies show a reasonable agreement with that of the experimental data demonstrating the reliability of the modified Galerkin approach used in this study for estimating the vibration. It is to note that the torsional case exhibits the highest deviation in Table 4. This is because the experimental $f_{1t}$ [16] corresponds to the uncoupled torsional mode, where the MGM [23] presented in this study yields the coupled torsion-dominated frequency leading to a larger mismatch than the flapping and lead–lag cases, for which the experimental results [16] are recorded as coupled frequencies.

3.2. Validation for Vibration Control

The vibration attenuation levels obtained in the theoretical LQR simulations closely matched those reported by Camino and Santos [38] who applied a periodic LQR to a fully actuated two-dimensional, four-blade rotor system. Table 5 illustrates the comparison of vibrations where the current findings are in good agreement with the results obtained by Camino and Santos [38].

In both studies, the flapping mode shows the highest relative reduction followed by the lead−lag mode (torsional deflection was not studied by them). The comparable percentages across the flapping and lead-lag modes suggest that the mathematical formulation presented in the current study captures the essential dynamics of the rotor blade vibration and its suppression. While Camino and Santos [38] implement their control law on a model with physical actuators and validate the stability via Floquet multiplier analysis, the approach in this research isolates the controller–system interaction from actuator limitations, allowing for a clean, parametric exploration of $\mathrm{Q}$−$\mathrm{R}$ trade-offs. The close agreement In attenuation trends reinforces confidence that the proposed control framework would exhibit similar performance in an actuator-integrated experimental setup. This consistency in trends and reduction percentages supports the validity of our theoretical approach and its relevance to future experimental implementation.

4. Results and Discussion

4.1. Vibration Analysis

4.1.1. Convergence Study

The FEA results for the dynamic response of the rotor blade are obtained after conducting a mesh convergence study focusing on fundamental natural frequencies. Figure 8 shows the variation of fundamental natural frequencies of the triply coupled vibration of the rotor blade with the total number of DOFs for nonrotating and rotating cases, where $f$ stands for the natural frequency in Hz with other subscripts as explained in Section 3.1. From Figure 8, as the number of DOFs increases, the trends of the natural frequencies become flattened after showing some initial changes indicating that the solution is already converged.

4.1.2. Natural Frequencies and Mode Shapes-Nonrotating Case

Table 6 lists the natural frequencies of free vibration of the triply coupled, nonrotating rotor blade for the first three modes obtained by the MGM [23] and the FEA using Abaqus (Dassault Systèmes, 2024). Here, $f_{nf}$, $f_{nl}$, and $f_{nt}$ indicate the $n$^th natural frequencies in Hz with f, l, and t denoting flapping, lead-lag, and torsional vibration modes, respectively. Table 6 shows a good agreement between the frequencies computed from the MGM and the FEA which is justified from the error analysis shown in Table 7.

From Table 7, the error level of $f_{nf}$ is slightly lower than that of $f_{nt}$. $f_{nl}$ has the least error level because of relatively higher in-plane bending stiffness than out-of-plane bending or torsional stiffness. It can be concluded that $f_{nt}$ would be less sensitive to error than $f_{nf}$ for higher modes. One reason for this is, for a particular mode, $f_{nt}$ is much higher in magnitude than $f_{nf}$ or $f_{nl}$. Another explanation is that the blade can be considered as a long, thin plate due to its geometry which is more prone to bending than torsion.

Figure 9 depicts the mode shapes of the triply coupled vibration of the nonrotating rotor blade obtained from the finite element simulation in Abaqus (Dassault Systèmes, 2024). The fundamental flapping and lead-lag mode shapes are given in Figure 9(a) and (b), respectively, with no node along the length of the blade while Figure 9(c) shows the fundamental torsion governed mode having a blue axis running parallel to the elastic axis of the blade. This axis goes through the shear center suggesting that the blade is about to rotate with respect to this axis. Figure 9(d)–(f) explain the similar phenomena for the second mode governed by the flapping, lead-lag, and torsion, respectively, each having one node. Following this, Figure 9(g)–(i) describe the third mode of the triply coupled vibration each having two nodes.

4.1.3. Natural Frequencies and Mode Shapes-Rotating Case

The natural frequencies for the first three modes of free vibration of the triply coupled, rotating helicopter blade obtained by the MGM [23] and FEA in Abaqus are shown in Table 8. From Table 8, flapping frequencies for the rotating case are much higher than that for the nonrotating case as listed in Table 6. This is due to the addition of the rotational stiffness to the structural stiffness developed from the angular velocity of the rotor blade. In contrast, minor differences are observed in the lead–lag and torsional frequencies between the nonrotating and rotating cases in Table 6 and Table 8, respectively, indicating only slight effect of centrifugal stiffening for those modes. Therefore, smaller frequency shifts for the lead-lag and torsional vibration modes between the rotating and nonrotating cases are observed and for the torsional vibration mode, it is almost negligible. The percentage error between the MGM [23] and FEA results for the rotating case remains within acceptable limits as presented in Table 9. The trends in error magnitude for the three different modes for the rotating blade are consistent with those observed for the nonrotating case.

Figure 10 describes the mode shapes of the triply coupled vibration of the rotating helicopter blade. These mode shapes are similar to that of Figure 9; however, in Figure 10, the modes are affected by the additional stiffness raised from the angular velocity of the rotor blade.

4.2. Rotor-Blade Vibration Control

After analyzing the vibration characteristics of the triply coupled rotor-blade system, the derived state-space models for both hovering and forward-flight conditions are numerically simulated under sinusoidal and step excitation forces, without control and subsequently with the LQR controller applied. For each case, the controlled and uncontrolled vibratory responses are presented on the same plot, enabling direct comparison of the controller’s effectiveness in attenuating flapping, lead–lag, and torsional motions. This study tuned the $\mathrm{Q}$ parameter, which is critical for achieving a balanced response, especially in systems where actuator capability or energy constraints are limiting factors. In theoretical studies, the systematic variation in $\mathrm{Q}$ provides insight into the achievable bounds of vibration attenuation without being restricted by physical actuator limits.

4.2.1. Controlled Response in the Hovering Flight

For hovering flight, two different types of forces, namely, periodic and impulse forces, are considered to visualize the controlled displacements at the tip of the helicopter rotor blade. Figure 11, Figure 12 and Figure 13 represent the actual steady-state flapping, lead-lag, and torsional deflections, respectively, and the corresponding controlled deflections at the tip of the helicopter rotor blade during the hovering flight. The displacements are estimated based on the contributions of the vibratory loadings generated from the aerodynamic lift, drag, and pitching moment with different values of $\mathrm{Q}$. The LQR controller is tuned to obtain the desired controlled deflection by varying the weighting matrix $\mathrm{Q}$ while keeping $\mathrm{R}=\mathrm{I}$. Figure 11 shows that for the flapping motion, the deflection of the blade tip fluctuates between -0.02 m and 0.02 m, and the nature is periodic due to the sinusoidal nature of the input force. The corresponding controlled flapping displacement shows peaks at 0.0046 m and $-$0.002 m for $\mathrm{Q}$ = 800$\mathrm{I}$.

Compared to the flapping response, the lead-lag motion shows a much smaller magnitude owing to the higher bending rigidity in the in-plane bending direction. Therefore, Figure 12 shows the amplitude of the original lead-lag deflection at the tip of the blade ranging from -0.0035 m to 0.004 m. After applying the LQR control, the corresponding peak amplitudes fluctuate between -0.0016 m and 0.0016 m. It should be noted that the magnitude of the drag force is much less than that of the lift force, and the corresponding flapping and lead-lag displacements are reflected in Figure 11 and Figure 12, respectively. Figure 13 shows the amplitudes of the original and controlled torsional deflections of the rotor blade at the tip owing to the torsional loading produced by the aerodynamic pitching moment and for $\mathrm{Q}=$ 30000$\mathrm{I}$, the trend of which is similar to Figure 11 and Figure 12.

Unlike the periodic input, Figure 14, Figure 15 and Figure 16 show the response of the state-space model for the step excitation force in terms of the controlled flapping, lead-lag, and torsional displacements for the hovering flight. It is evident that the controlled peak deflection for step excitation is higher than that of sinusoidal excitation for the flapping case (Figure 11). Similar phenomena are also observed for the lead-lag and torsional deflection cases in Figure 12 and Figure 13, respectively.

4.2.2. Controlled Response in the Forward Flight

In forward flight, the helicopter rotor blade experiences asymmetric velocity distributions and time-varying periodic loads based on the speed of the helicopter, which, in another form, is represented by the advanced ratio. Therefore, the AoA in forward flight varies with time or with the azimuth angle. In this study, the analysis is demonstrated for an AoA of 4$°$ and an advanced ratio of 0.3 and can be similarly performed for any other value of the AoA to calculate the controlled response. Figure 17 shows the maximum flapping deflection of 0.9 m at the blade tip, which is measured from the hub plane of the rotor. The corresponding controlled deflection gets level at 0.4 m for $\mathrm{Q}=$ 100$\mathrm{I}$. In this case, both the original and controlled flapping deflections show higher magnitudes than the hovering case because of the inclusion of the displacement generated by the coning angle measured from the hub plane.

Figure 18 shows the uncontrolled and controlled deflections for the lead-lag case for $\mathrm{Q}=$ 500$\mathrm{I}$ which reveals that the magnitudes of the displacements are much lower than those of the flapping case because of the smaller magnitude of the drag force compared to the lift force. Figure 19 shows the controlled blade tip deflection owing to the torsional moment ranging from 0.026 rad to 0.037 rad, which is observed as harmonic in nature.

Similarly, for the impulse load during the hovering flight, the LQR controller is finally used to simulate the controlled displacements for all three DOFs of motion for the forward flight with a step excitation force. Figure 20, Figure 21 and Figure 22 depict the original flapping, lead-lag, and torsional displacements of the blade tip, respectively, owing to the step input loading along with the corresponding controlled displacements. From the simulation results, it can be inferred that the behavior trend for each DOF of motion is similar to that for hovering flight subjected to step excitation. The control of the flapping displacement is always faster than the lead-lag and torsional displacements because of the larger magnitude of the aerodynamic damping for the flapping motion than that of the lead-lag or torsional motions.

Figure 23 shows the variation in the helicopter rotor vibration amplitudes for flapping, lead–lag, and torsional motions with the LQR weighting parameter $\mathrm{Q}$. The amplitudes are normalized by the maximum value of each case across the $\mathrm{Q}$ sweep, enabling direct comparison of attenuation trends among the three vibrational DOFs of motion. Figure 23 shows that increasing the Q parameter leads to a monotonic reduction in normalized vibration amplitude for all three modes, confirming the effectiveness of the LQR weighting in penalizing excessive state excursions. Among the three, the torsional response decreases most rapidly, suggesting that torsional motion is more sensitive to the state penalty, while the lead–lag response attenuates more gradually, reflecting its greater inertia-dominated behavior. The flapping mode lies in between, with a clear but lesser steep attenuation trend than that of the torsional response. Together, these results demonstrate that careful tuning of the parameter $\mathrm{Q}$ allows targeted suppression of specific vibration modes and validates the robustness of the controller across coupled, multiple DOFs of motion.

At extremely high values of $\mathrm{Q}$, the controller penalizes state excursions so strongly that the system response saturates, with further increases in $\mathrm{Q}$ producing negligible additional vibration reductions. Such over-penalization not only yields diminishing returns, but also renders the controller overly conservative, diminishing its ability to respond effectively to dynamic disturbances. By contrast, the optimal range of $\mathrm{Q}=$ ${10}^2$–${10}^4$ provides a balanced trade-off between attenuation and responsiveness, ensuring optimum vibration suppression without unnecessary conservatism.

5. Conclusions

In this study, a comprehensive analysis is conducted to explore the vibration characteristics of a flexible, composite helicopter rotor blade and develop the corresponding state-space model to attenuate the vibration that considers coupled, three degrees-of-freedom motion of the rotor blade. The sectional properties of the rotor blade are estimated analytically and used in the modified Galerkin and finite element approaches to find the natural frequencies of vibration. The orthogonality condition for the triply coupled vibration is formulated to derive the state-space equations of the optimal controller. Controlled vibration responses in both hovering and forward-flight conditions under step and sinusoidal excitations confirm the effectiveness of the Linear Quadratic Regulator (LQR), yielding 60–90% and 50–70% reductions in vibration for periodic and impulsive loads, respectively, and showing substantial controllability for attenuating flapping, lead-lag, and torsional vibrations.

The resulting state-space model provides a general and adaptable framework capable of predicting controlled displacements at any spanwise location of the rotor blade for rotating and nonrotating configurations, provided modal parameters are known. Although the present work focuses on the first vibration mode, the methodology can readily be extended to higher-order modes by including relevant natural frequencies and damping ratios. Consistent with prior studies, the results confirm the advantages of the optimal control over classical approaches, offering enhanced stability and reduced tuning effort. Systematic tuning of the LQR weighting parameter revealed an optimal range of 10²−10⁴, providing a practical design guideline for future actuator-based implementations.

Overall, the proposed analytical framework explains the coupled vibration behavior of the composite rotor-blade system, demonstrates substantial vibration suppression across all modes, and establishes a robust theoretical foundation for validating more advanced observer-based, adaptive, or actuator-integrated control strategies. Future work will incorporate unsteady aerodynamic effects to improve the forward-flight fidelity and will focus on the experimental validation, actuator integration, and development of hybrid control architectures such as fuzzy-augmented or PID-assisted schemes to further enhance vibration suppression performance.