A harmonic balance solution for the intrinsic 1D nonlinear equations of the beams

Abstract

In this paper, a harmonic balance method is introduced to solve the intrinsic 1D nonlinear equations presented by Hodges for analyzing initially curved and twisted anisotropic rotating beams. First, the nonlinear first-order partial differential equations are discretized in space domain using the Galerkin approach. Then, the time response of the system of equations is approximated by using first and the second harmonic terms under the harmonic excitations applied to the structure. The harmonic balance algorithm has been developed for the linearized system of equations about the steady-state solution and for fully nonlinear equations. The differences between these two systems have been addressed under consideration of different levels of external loads. Fully analytical formulae have been presented for the HB solution of the linearized equations in this work. In the proposed approach for the nonlinear system, the Jacobian matrix utilized in the numerical solution part is obtained analytically to increase the time efficiency of the solution method. The results of the proposed approach for an initially twisted and curved blade are compared with a different algorithm based on a finite difference approach to discretize the equations and a time marching method to solve the system of equations in time domain. The comparison between the outputs of the harmonic balance method and the time marching method demonstrate high correlations between these distinct methods. It should be noted that the accuracy of the time marching method solution has been evaluated in another publication by the authors using the output of FLIGHTLAB for a blade with different types of boundary conditions. The proposed harmonic balance-based method delivers a very time-efficient approach for time analysis of twisted/curved beams and blades under harmonic loads in comparison with other employed methods. The developed approach can be a powerful computational tool for design optimization of composite blades.

1. Introduction

Recent demand and interest in aerial mobility enhance the importance of analytical and computation tools for designing light and effective components like composite blades and wings. One-dimensional (1D) models play a crucial role in structural analysis considering their economical computational efforts and they could provide the designer with simple tools to analyze numerous problems. One-dimensional beam models can be investigated using geometrically nonlinear regimes that include different levels of accuracy to evaluate the static and dynamic characteristics. These equations have been developed based on a fully intrinsic formulation and they do not need the displacement and rotation variables. Therefore, they are free of infinite-degree nonlinearities found in other types of formulations.

In recent years, the variational-asymptotic method (VAM) has received considerable interest from researchers in the field of dynamics. Variational-asymptotic method is computationally more efficient than 3D FEA.

To provide the cross-sectional properties for 1D models, the VAM has been used. The process starts by considering the elastic energy functional based on small parameters associated with slender beams. Then, the problem is solved by minimizing this function. Berdichevsky (1981) and Hodges (1990) appear to be the first to decompose a 3D geometrically nonlinear elasticity analysis of slender structures into a nonlinear 1D analysis in the spanwise direction and a cross-sectional 2D analysis. There are several publications and more details about this analysis (Cesnik, 1994; Cesnik and Hodges, 1997; Hodges 2003; Yu et al., 2002). According to these references, the developed solution for the cross-sectional analysis is based on a perturbation method. This approximation solution method for the set of nonlinear equations has been explained in detail by Cesnik (1994) and Hodges (2006).

Many papers have been published on the general subject of aeroelasticity using the 1D model of beams and blades (Cesnik et al., 1996; Ghorashi and Nitzsche, 2007; Hodges et al., 1991; Ho et al., 2010; Shang et al., 1999; Yu et al., 2012). A comprehensive review of the state-of-the art by focusing on the vibrations, dynamics, and control aspects of rotating blade was provided by Ghorashi and Nitzsche (2008). The nonlinear dynamic response of an accelerating composite rotor blade using perturbations was studied (Ghorashi and Nitzsche, 2009). The energy-consistent Galerkin approach was developed by Patil and Althoff (2011) for the nonlinear dynamics of beams using intrinsic equations. The fully intrinsic first-order equations were derived considering different configurations for the connection between the hub and the blade by Sotoudeh and Hodges (2013). Verification of the presented formulations in one specific configuration was done by checking the beam natural frequencies. A computational model was also developed to study the nonlinear steady-state static response and free vibration of fiber multiscale composite beams using the intrinsic 1D nonlinear equations (Rafiee et al., 2016). The structural dynamics analysis of rotating helicopter blades with non-uniform properties subjected to realistic boundary conditions using 1D nonlinear equations was presented by Siami et al. (2021). Flexible joints were employed between the hub and blade, and the results were evaluated by comparing the time responses and the Campbell diagrams against the commercial package FLIGHTLAB.

Other beam theories have been used for structural analysis of composite beams. The post-buckling of simply supported functionally graded material (FGM) beams have been studied using different theories by Amara et al. (2016). The results showed the advantages of higher-order theories like parabolic and exponential shear deformation beam theories compared to the classical beam theory. A refined shear deformation theory was employed for thermo-mechanical buckling analysis of laminated composite beams (Derbale et al., 2021). The thermal effects on the behavior of composite beams were investigated by Bouazza et al. (2019a). An analytical approach based on theory of hyperbolic refined shear deformation was developed for buckling of composite beams (Bouazza et al., 2019b; Bouazza and Zenkour, 2020). The effect of different parameters was investigated on thermal buckling of laminated composite beams. The thermal buckling behavior of the functionality graded beam attached with piezoelectric layers was studied by Ellali et al. (2022) using the shear deformation beam theory.

The harmonic balance method (HBM) is extensively used in different vibration simulations. A completely numerical harmonic balance approach was used for solving the equation of bladed rotor (Peters and Morillo, 1997). The generalized HBM was used by Luo and Huang (2011) for analytical approximate solutions of periodic motions in nonlinear dynamical systems, and the nonlinear damping, periodically forced, Duffing oscillator was investigated as an example. Multi-mode coupled differential equations have been developed for simply supported and clamped-hinged flexible beams on an elastic foundation subject to an external harmonic excitation by Hasan et al. (2014). The result from the proposed HBM agrees well with the result from Runge–Kutta’s numerical integration method. Hasan and Rahman (2022) used the multi-level residue HBM for solving the geometrically cubic nonlinear vibration of the simply supported beam. The results were evaluated using a numerical integration method. Fitzgibbon et al. (2019) employed the HBM for rotor blade performance predictions. This study was presented for two types of blades in hover and forward flight. The outputs of HBM were compared with results from time marching simulations. The authors showed that the HBM is a promising alternative to computationally expensive time marching simulations. The predictions of the rotor loading at significantly reduced computational costs in both hover and forward flight were very good.

In this paper, the HBM is used to solve the 1D nonlinear equations under the harmonic loads. It should be noted that based on the knowledge of the authors, the HBM has not yet been applied for solving the intrinsic 1D Hodges equation. In other words, the main contribution of this work is to develop a harmonic balance solution method for the 1D nonlinear intrinsic equations of beams. The proposed methods include an analytical solution for a linearized form of the equations and a numerical solution scheme with an analytical form of the Jacobian matrix for the numerical approach. The equations are discretized in space domain by using a Galerkin approach. Then, the first and second harmonic terms are considered to approximate the time response of the system under the harmonic loads. The Jacobian matrix of the numerical solution part is derived analytically to increase the time efficiency of the introduced solution approach. The accuracy of the solution is evaluated by comparing the results with the outputs of a time marching method (TMM) evaluated previously by the authors (Siami et al., 2021). The harmonic balance response for the linearized system of equations about the steady-state solution is presented using analytical formulations and the outputs are compared with the solution of the nonlinear system. The main advantage of this method is the time efficiency of the proposed solutions compared to the other existing solution methods for the 1D nonlinear equations. The proposed solutions can be very useful in design process of the blades and also in checking the effect of control strategies on the structural dynamics of the beams and blades. The proposed methods are provided time responses of the structure up to the second harmonic terms of external loads.

This paper is organized as follows: In Section 2, the 1D nonlinear equations and the discretization process in space domain are explained. The HBM used for time analysis of the equations is presented in this section. The results of the proposed method are discussed in Section 3. In addition, in this section, the accuracy of the HB solutions is evaluated by comparing some calculated variables with those obtained from the TMM. Section 4 collects the conclusion of this paper.

2. Theory

2.1. Nonlinear fully intrinsic beam equations

The intrinsic equations governing the dynamics of a general, non-uniform, twisted, curved, anisotropic beam undergoing large deformation are given as follows (Hodges, 1990):

$F^{\prime }+\left(\stackrel{\sim }{k}+\stackrel{\sim }{K}\right)F+f_{ext}=\dot{P}+\stackrel{\sim }{\mathrm{\Omega }}P$

$M^{\prime }+\left(\stackrel{\sim }{k}+\stackrel{\sim }{K}\right)M+\left({\stackrel{\sim }{e}}_1+\stackrel{\sim }{\Gamma }\right)F+m_{ext}=\dot{H}+\stackrel{\sim }{\mathrm{\Omega }}H+\stackrel{\sim }{V}P$

$V^{\prime }+\left(\stackrel{\sim }{k}+\stackrel{\sim }{K}\right)V+\left({\stackrel{\sim }{e}}_1+\stackrel{\sim }{\Gamma }\right)\mathrm{\Omega }=\dot{\Gamma }$

${\mathrm{\Omega }}^{\prime }+\left(\stackrel{\sim }{k}+\stackrel{\sim }{K}\right)\mathrm{\Omega }=\dot{K}$

(1)

where ${\left(\right)}^{\prime }$ denotes the derivative with respect to the beam reference line and $\dot{\left(\right)}$ is the absolute time derivative. $F\left(x,t\right)$ and $M\left(x,t\right)$ denote the inertial force and moment vector, $P\left(x,t\right)$ and $H\left(x,t\right)$ are the linear and angular momentum vectors, $\Gamma \left(x,t\right)$ and $K\left(x,t\right)$ are the beam strain and local curvatures. In addition, $V\left(x,t\right)$ and $\mathrm{\Omega }\left(x,t\right)$ are the linear and angular velocities. The external force and moment caused by aerodynamic effects, for example, are denoted by $f_{ext}\left(x,t\right)$ and $m_{ext}\left(x,t\right)$ , respectively. The initial curvature and twist of the beam are grouped in the vector ${k=\left[\begin{array}{ccc}{k}_1& k_2& k_3\end{array}\right]}^T$ and $e_1={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T$ is the unit vector in the axial direction, $x$ according to Figure 1. It should be noted that the tilde transforms the vector cross operation $y\times$ into its matrix multiplication equivalent, $\stackrel{\sim }{y}$ .

Figure 1.

A twisted blade and its coordinate system.

The intrinsic beam equations provide 4 vector equations for 8 vector unknowns, $\left(F,M,P,H,\Gamma ,K,V,\mathrm{\Omega }\right)$ . Therefore, four more vector equations are needed. Two equations relate the generalized forces ( $F,M$ ) and the generalized strains ( $\Gamma$ , K) via the beam cross-section stiffness matrix. The relationships between the generalized momenta ( $P,H$ ) and the generalized velocities ( $V,\mathrm{\Omega }$ ) can be built using the beam cross-section inertia matrix. These additional equations are (Hodges, 1990):

$\left\{\begin{array}{c}F\\ M\end{array}\right\}={\left[\begin{array}{cc}\mathbb{R}& \mathbb{S}\\ {\mathbb{S}}^T& \mathbb{T}\end{array}\right]}^{-1}\left\{\begin{array}{c}\mathrm{\Gamma }\\ K\end{array}\right\}$

(2)

$\left\{\begin{array}{c}P\\ H\end{array}\right\}=\left[\begin{array}{cc}\mathbb{G}& \mathbb{K}\\ {\mathbb{K}}^T& \mathbb{I}\end{array}\right]\left\{\begin{array}{c}V\\ \mathrm{\Omega }\end{array}\right\}$

where $\mathbb{R}$ , $\mathbb{S}$ , and $\mathbb{T}$ are the cross-sectional flexibilities. The inertia-related matrices have the following form:

$\mathbb{G}=\mu I=\left[\begin{array}{ccc}\mu & 0& 0\\ 0& \mu & 0\\ 0& 0& \mu \end{array}\right],\mathbb{K}=\left[\begin{array}{ccc}0& \mu {\overline{\xi }}_3& -\mu {\overline{\xi }}_2\\ -\mu {\overline{\xi }}_3& 0& 0\\ \mu {\overline{\xi }}_2& 0& 0\end{array}\right],\mathbb{I}=\left[\begin{array}{ccc}{\mathbb{i}}_2+{\mathbb{i}}_3& 0& 0\\ 0& {\mathbb{i}}_2& {\mathbb{i}}_{23}\\ 0& {\mathbb{i}}_{23}& {\mathbb{i}}_3\end{array}\right]$

(3)

where $\mu$ , ${\overline{\xi }}_2$ , ${\overline{\xi }}_3$ , and ${\mathbb{i}}_2$ , ${\mathbb{i}}_3$ and ${\mathbb{i}}_{23}$ are the mass per unit length, offsets from the reference line of the cross-sectional mass centroid, and the two cross-sectional mass momenta and the product of inertia per unit length, respectively.

The shifted Legendre polynomials are used here as the independent trial functions. These functions are presented as follows:

$P^0\left(\overline{x}\right)=1,P^1\left(\overline{x}\right)=2\overline{x}-1$

$P^{k+1}\left(\overline{x}\right)=\frac{\left(2k+1\right)\left(2\overline{x}-1\right)P^k\left(\overline{x}\right)-k{P}^{k-1}\left(\overline{x}\right)}{k+1}$

(4)

within the normalized (shifted) interval $0\le \overline{x}\le 1$ .

These functions are orthogonal over the shifted $\left(\mathrm{0,1}\right)$ interval:

${\int }_0^1{P}^i\left(\overline{x}\right)P^j\left(\overline{x}\right)d\overline{x}={\delta }_{ij}\frac{1}{2i+1}$

(5)

In addition, ${\mathrm{\Phi }}_k\left(\overline{x}\right)$ are the mode shapes along the 3 directions, organized as follows:

${\mathrm{\Phi }}^k\left(\overline{x}\right)=\left[\begin{array}{ccc}{\varphi }^k\left(\overline{x}\right)& 0& 0\\ 0& {\varphi }^k\left(\overline{x}\right)& 0\\ 0& 0& {\varphi }^k\left(\overline{x}\right)\end{array}\right]=\left[\begin{array}{ccc}{P}^{k-1}\left(\overline{x}\right)& 0& 0\\ 0& P^{k-1}\left(\overline{x}\right)& 0\\ 0& 0& P^{k-1}\left(\overline{x}\right)\end{array}\right]$

(6)

where $\overline{x}=x/L$ and $L$ is the length of the beam along the span and $k$ presents the order of the mode shape. The parameters that will be the solution of discretized system can be found approximately by separating variables in the space and time domains. So, the unknown variables at each element are written as

$V\left(x,t\right)=\left[\begin{array}{ccc}\begin{array}{cc}{\mathrm{\Phi }}^1\left(x\right)& {\mathrm{\Phi }}^2\left(x\right)\end{array}& \dots & {\mathrm{\Phi }}^m\left(x\right)\end{array}\right]{\left[\begin{array}{cccc}{v}^1\left(t\right)& v^2\left(t\right)& \dots & v^m\left(t\right)\end{array}\right]}^T\hspace{0.17em}\Omega \left(x,t\right)=\left[\begin{array}{ccc}\begin{array}{cc}{\mathrm{\Phi }}^1\left(x\right)& {\mathrm{\Phi }}^2\left(x\right)\end{array}& \dots & {\mathrm{\Phi }}^m\left(x\right)\end{array}\right]{\left[\begin{array}{cccc}{\omega }^1\left(t\right)& {\omega }^2\left(t\right)& \dots & {\omega }^m\left(t\right)\end{array}\right]}^T\hspace{0.17em}\Gamma \left(x,t\right)=\left[\begin{array}{ccc}\begin{array}{cc}{\mathrm{\Phi }}^1\left(x\right)& {\mathrm{\Phi }}^2\left(x\right)\end{array}& \dots & {\mathrm{\Phi }}^m\left(x\right)\end{array}\right]{\left[\begin{array}{cccc}{\gamma }^1\left(t\right)& {\gamma }^2\left(t\right)& \dots & {\gamma }^m\left(t\right)\end{array}\right]}^T\hspace{0.17em}{\rm K}\left(x,t\right)=\left[\begin{array}{ccc}\begin{array}{cc}{\mathrm{\Phi }}^1\left(x\right)& {\mathrm{\Phi }}^2\left(x\right)\end{array}& \dots & {\mathrm{\Phi }}^m\left(x\right)\end{array}\right]{\left[\begin{array}{cccc}{\kappa }^1\left(t\right)& {\kappa }^2\left(t\right)& \dots & {\kappa }^m\left(t\right)\end{array}\right]}^T$

(7)

where in the above equations, the superscripts indicate the order of the trial function. It should be noted in this work the assumed modes for each variable and each direction are chosen to be identical. Based on equation (7), the time-dependent part of the unknown variables of each element is defined as

$q\left(t\right)={\left[v^1\left(t\right),{\omega }^1\left(t\right),{\gamma }^1\left(t\right),{\kappa }^1\left(t\right),\dots ,v^m\left(t\right),{\omega }^m\left(t\right),{\gamma }^m\left(t\right),{\kappa }^m\left(t\right)\right]}^T$

(8)

Considering a cantilever beam, the boundary conditions of the blade are chosen to constrain the generalized velocity at the root ( $x=0$ ) and the generalized forces at the tip ( $x=L$ ) of the beam. The boundary conditions can be written as follows:

$V\left(0\right)=V_0,\mathrm{\Omega }\left(0\right)={\mathrm{\Omega }}_0,F\left(L\right)=F_L,M\left(L\right)=M_L$

(9)

Assuming uniform structural properties along each element and using the separation of variables given in equation (7), the intrinsic 1D nonlinear equations presented in equation (1) can be written as (Siami et al., 2021)

$A^{ki}{\dot{q}}^i\left(t\right)+B^{ki}{q}^i\left(t\right)+C^{kij}{q}^i\left(t\right)q^j\left(t\right)+D^k+F_{ext}^k\left(t\right)+M_{ext}^k\left(t\right)=0$

(10)

where $q^i\left(t\right)$ has been defined in equation (8). It should be noted that $F_{ext}^k\left(t\right)$ and $M_{ext}^k\left(t\right)$ are the terms related to external loads. They will appear in the equations if external distributed loads are applied to the elements. The matrix coefficients presented in equation (10) can be found in another study by Siami et al. (2021).

For modal analysis, the equation of the system needs to be linearized around the steady-state solution. The steady-state solution of equation (10) can be determined by solving the following equations:

$B^{ki}{q}_s^i+C^{kij}{q}_s^i{q}_s^j+{\overline{D}}^k=0$

(11)

where the vector $q^0$ collects the steady-state solution due to steady-state forcing vector $\overline{D}$ . This solution is obtained using the Newton–Raphson iterative method. The Jacobian required for the iteration process can be easily verified to be

${\mathcal{J}}_s^{ki}\left(q_s\right)=B^{ki}+C^{kij}{q}_s^j+C^{kji}{q}_s^j$

(12)

The one-dimensional beam equations are linearized about the equilibrium state. For this purpose, the variable vector can be expressed in terms of a perturbation about the steady-state condition as follows:

$q\left(t\right)=q_s+q_p\left(t\right)$

(13)

where $q_p\left(t\right)$ is a small perturbation about the steady-state solution. By introducing equation (13) into equation (10), the linearized 1D nonlinear equation about the steady-state solution can be obtained:

${\widehat{A}}^{ki}{\dot{q}}_p^i+{\widehat{B}}^{ki}{q}_p^i=0$

(14)

where

${\widehat{A}}^{ki}=A^{ki},{\widehat{B}}^{ki}=B^{ki}+\left(C^{kij}+C^{kji}\right)q_s^j$

(15)

The natural frequencies and corresponding mode shapes are identified using the unforced blade model given by equation (14).

2.2. Harmonic balance method

In this section, the formulations of the deployed HB approach for the intrinsic 1D nonlinear equations given by equation (10) are presented. The periodic external forces and moments applied to the blade are considered as follows:

${\mathcal{F}}_{ext}^k\left(t\right)={\mathcal{F}}_{C1}^k\cos \mathrm{\Omega }t+{\mathcal{F}}_{C2}^k\cos 2\mathrm{\Omega }t+{\mathcal{F}}_{S1}^k\sin \mathrm{\Omega }t+{\mathcal{F}}_{S2}^k\sin 2\mathrm{\Omega }t\hspace{0.17em}{\mathcal{M}}_{ext}^k\left(t\right)={\mathcal{M}}_{C1}^k\cos \mathrm{\Omega }t+{\mathcal{M}}_{C2}^k\cos 2\mathrm{\Omega }t+{\mathcal{M}}_{S1}^k\sin \mathrm{\Omega }t+{\mathcal{M}}_{S2}^k\sin 2\mathrm{\Omega }t$

(16)

The approximated periodic solution for the nonlinear intrinsic equation given in equation (10) under the harmonic loads given by equation (16) can be found by using the first and the second harmonic solutions:

$q^i\left(t\right)=q_h^i\left(t\right)=Q_0^i+Q_1^i\cos \mathrm{\Omega }t+S_1^i\hspace{0.17em}\sin \mathrm{\Omega }t+Q_2^i\cos 2\mathrm{\Omega }t+S_2^i\hspace{0.17em}\sin 2\mathrm{\Omega }t$

(17)

Substituting equations (16) and (17) into equation (10) yields the following relation:

$A^{ki}\left(-\mathrm{\Omega }{Q}_1^i\hspace{0.17em}\sin \mathrm{\Omega }t+\mathrm{\Omega }{S}_1^i\cos \mathrm{\Omega }t-2\mathrm{\Omega }{Q}_2^i\hspace{0.17em}\sin 2\mathrm{\Omega }t+2\mathrm{\Omega }{S}_2^i\cos 2\mathrm{\Omega }t\right)+B^{ki}\left(Q_0^i+Q_1^i\cos \mathrm{\Omega }t+S_1^i\hspace{0.17em}\sin \mathrm{\Omega }t+Q_2^i\cos 2\mathrm{\Omega }t+S_2^i\hspace{0.17em}\sin 2\mathrm{\Omega }t\right)+D^k+C^{kij}{\mathcal{T}}^{ij}+{\mathcal{F}}_{C1}^k\cos \mathrm{\Omega }t+{\mathcal{F}}_{C2}^k\cos 2\mathrm{\Omega }t+{\mathcal{F}}_{S1}^k\sin \mathrm{\Omega }t+{\mathcal{F}}_{S2}^k\sin 2\mathrm{\Omega }t+{\mathcal{M}}_{C1}^k\cos \mathrm{\Omega }t+{\mathcal{M}}_{C2}^k\cos 2\mathrm{\Omega }t+{\mathcal{M}}_{S1}^k\sin \mathrm{\Omega }t+{\mathcal{M}}_{S2}^k\sin 2\mathrm{\Omega }t=0$

(18)

where

${\mathcal{T}}^{ij}=Q_0^i{Q}_0^j+Q_0^i{Q}_1^j\cos \mathrm{\Omega }t+Q_0^i{S}_1^j\sin \mathrm{\Omega }t+Q_0^i{Q}_2^j\cos 2\mathrm{\Omega }t+Q_0^i{S}_2^j\sin 2\mathrm{\Omega }t+Q_1^i{Q}_0^j\cos \mathrm{\Omega }t+Q_1^i{Q}_1^j{\cos }^2\mathrm{\Omega }t+Q_1^i{S}_1^j\cos \mathrm{\Omega }t\sin \mathrm{\Omega }t+Q_1^i{Q}_2^j\cos \mathrm{\Omega }t\cos 2\mathrm{\Omega }t+Q_1^i{S}_2^j\cos \mathrm{\Omega }t\sin 2\mathrm{\Omega }t+{S_1^{i}Q}_0^j\sin \mathrm{\Omega }t+{S_1^{i}Q}_1^j\sin \mathrm{\Omega }t\cos \mathrm{\Omega }t+S_1^i{S}_1^j{\sin }^2\mathrm{\Omega }t+S_1^i{Q}_2^j\sin \mathrm{\Omega }t\cos 2\mathrm{\Omega }t+{S_1^{i}S}_2^j\sin \mathrm{\Omega }t\sin 2\mathrm{\Omega }t+{Q_2^{i}Q}_0^j\cos 2\mathrm{\Omega }t+{Q_2^{i}Q}_1^j\cos 2\mathrm{\Omega }t\cos \mathrm{\Omega }t+{Q_2^{i}S}_1^j\cos 2\mathrm{\Omega }t\sin \mathrm{\Omega }t+{Q_2^{i}Q}_2^j{\cos }^{2}2\mathrm{\Omega }t+{Q_2^{i}S}_2^j\cos 2\mathrm{\Omega }t\sin 2\mathrm{\Omega }t+{S_2^{i}Q}_0^j\sin 2\mathrm{\Omega }t+{S_2^{i}S}_1^j\sin 2\mathrm{\Omega }t\sin \mathrm{\Omega }t+{S_2^{i}Q}_2^j\sin 2\mathrm{\Omega }t\cos 2\mathrm{\Omega }t+{S_2^{i}S}_2^j{\sin }^{2}2\mathrm{\Omega }t$

(19)

The above equation can be rewritten as

${\mathcal{T}}^{ij}=Q_0^i{Q}_0^j+\left(Q_0^i{Q}_1^j+Q_1^i{Q}_0^j\right)\cos \mathrm{\Omega }t+\left(Q_0^i{S}_1^j+S_1^i{Q}_0^j\right)\sin \mathrm{\Omega }t+\left(Q_0^i{Q}_2^j+Q_2^i{Q}_0^j\right)\cos 2\mathrm{\Omega }t+\left(Q_0^i{S}_2^j+S_2^i{Q}_0^j\right)\sin 2\mathrm{\Omega }t+Q_1^i{Q}_1^j\left(\frac{\cos 2\mathrm{\Omega }t+1}{2}\right)+Q_1^i{S}_1^j\left(\frac{\sin 2\mathrm{\Omega }t}{2}\right)+Q_1^i{Q}_2^j\left(\frac{\cos 3\mathrm{\Omega }t+\cos \mathrm{\Omega }t}{2}\right)+Q_1^i{S}_2^j\left(\frac{\sin \mathrm{\Omega }t+\sin 3\mathrm{\Omega }t}{2}\right)+S_1^i{Q}_1^j\left(\frac{\sin 2\mathrm{\Omega }t}{2}\right)+S_1^i{S}_1^j\left(\frac{1-\cos 2\mathrm{\Omega }t}{2}\right)+S_1^i{Q}_2^j\left(\frac{\sin 3\mathrm{\Omega }t-\sin \mathrm{\Omega }t}{2}\right)+S_1^i{S}_2^j\left(\frac{\cos \mathrm{\Omega }t-\cos 3\mathrm{\Omega }t}{2}\right)+Q_2^i{Q}_1^j\left(\frac{\cos 3\mathrm{\Omega }t+\cos \mathrm{\Omega }t}{2}\right)+Q_2^i{S}_1^j\left(\frac{\sin 3\mathrm{\Omega }t-\sin \mathrm{\Omega }t}{2}\right)+Q_2^i{Q}_2^j\left(\frac{\cos 4\mathrm{\Omega }t+1}{2}\right)+Q_2^i{S}_2^j\left(\frac{\sin 4\mathrm{\Omega }t}{2}\right)+S_2^i{Q}_1^j\left(\frac{\sin \mathrm{\Omega }t+\sin 3\mathrm{\Omega }t}{2}\right)+S_2^i{S}_1^j\left(\frac{\cos \mathrm{\Omega }t-\cos 3\mathrm{\Omega }t}{2}\right)+S_2^i{Q}_2^j\left(\frac{\sin 4\mathrm{\Omega }t}{2}\right)+S_2^i{S}_2^j\left(\frac{1-\cos 4\mathrm{\Omega }t}{2}\right)$

(20)

Considering the first and the second harmonic terms and ignoring the higher-order terms, equation (20) can be represented as

${\mathcal{T}}^{ij}=Q_0^i{Q}_0^j+\left(Q_0^i{Q}_1^j+Q_1^i{Q}_0^j\right)\cos \mathrm{\Omega }t+\left(Q_0^i{S}_1^j+S_1^i{Q}_0^j\right)\sin \mathrm{\Omega }t+\left(Q_0^i{Q}_2^j+Q_2^i{Q}_0^j\right)\cos 2\mathrm{\Omega }t+\left(Q_0^i{S}_2^j+S_2^i{Q}_0^j\right)\sin 2\mathrm{\Omega }t+Q_1^i{Q}_1^j\left(\frac{\cos 2\mathrm{\Omega }t+1}{2}\right)+Q_1^i{S}_1^j\left(\frac{\sin 2\mathrm{\Omega }t}{2}\right)+Q_1^i{Q}_2^j\left(\frac{\cos \mathrm{\Omega }t}{2}\right)+Q_1^i{S}_2^j\left(\frac{\sin \mathrm{\Omega }t}{2}\right)+S_1^i{Q}_1^j\left(\frac{\sin 2\mathrm{\Omega }t}{2}\right)+S_1^i{S}_1^j\left(\frac{1-\cos 2\mathrm{\Omega }t}{2}\right)+S_1^i{Q}_2^j\left(-\frac{\sin 2\mathrm{\Omega }t}{2}\right)+S_1^i{S}_2^j\left(\frac{\cos \mathrm{\Omega }t}{2}\right)+Q_2^i{Q}_1^j\left(\frac{\cos \mathrm{\Omega }t}{2}\right)+Q_2^i{S}_1^j\left(-\frac{\sin \mathrm{\Omega }t}{2}\right)+Q_2^i{Q}_2^j\left(\frac{1}{2}\right)+S_2^i{Q}_1^j\left(\frac{\sin \mathrm{\Omega }t}{2}\right)+S_2^i{S}_1^j\left(\frac{\cos \mathrm{\Omega }t}{2}\right)+S_2^i{S}_2^j\left(\frac{1}{2}\right)$

(21)

Substituting equation (21) into equation (18) yields the following equations:

$B^{ki}{Q}_0^i+D^k+C^{kij}\left(Q_0^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_1^j+\frac{1}{2}{S}_1^i{S}_1^j+\frac{1}{2}{Q}_2^i{Q}_2^j+\frac{1}{2}{S}_2^i{S}_2^j\right)=0$

(22)

$\left[{\mathrm{\Omega }A}^{ki}{S}_1^i+B^{ki}{Q}_1^i+C^{kij}\left(Q_0^i{Q}_1^j+Q_1^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_2^j+\frac{1}{2}{S}_1^i{S}_2^j+\frac{1}{2}{Q}_2^i{Q}_1^j+\frac{1}{2}{S}_2^i{S}_1^j\right)+{\mathcal{F}}_{c1}^k+{\mathcal{M}}_{c1}^k\right]\cos \mathrm{\Omega }t=0$

(23)

$\left[{-\mathrm{\Omega }A}^{ki}{Q}_1^i+B^{ki}{S}_1^i+C^{kij}\left(Q_0^i{S}_1^j+S_1^i{Q}_0^j+\frac{1}{2}{Q}_1^i{S}_2^j-\frac{1}{2}{S}_1^i{Q}_2^j-\frac{1}{2}{Q}_2^i{S}_1^j+\frac{1}{2}{S}_2^i{Q}_1^j\right)+{\mathcal{F}}_{s1}^k+{\mathcal{M}}_{s1}^k\right]\sin \mathrm{\Omega }t=0$

(24)

$\left[2{\mathrm{\Omega }A}^{ki}{S}_2^i+B^{ki}{Q}_2^i+C^{kij}\left(Q_0^i{Q}_2^j+Q_2^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_1^j-\frac{1}{2}{S}_1^i{S}_1^j\right)+{\mathcal{F}}_{c2}^k+{\mathcal{M}}_{c2}^k\right]\cos 2\mathrm{\Omega }t=0$

(25)

$\left[-2{\mathrm{\Omega }A}^{ki}{Q}_2^i+B^{ki}{S}_2^i+C^{kij}\left(Q_0^i{S}_2^j+S_2^i{Q}_0^j+\frac{1}{2}{Q}_1^i{S}_1^j+\frac{1}{2}{S}_1^i{Q}_1^j\right)+{\mathcal{F}}_{s2}^k+{\mathcal{M}}_{s2}^k\right]\sin 2\mathrm{\Omega }t=0$

(26)

According to the above equations, the following relations can be written:

$E_1\left(X\right)=B^{ki}{Q}_0^i+D^k+C^{kij}\left(Q_0^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_1^j+\frac{1}{2}{S}_1^i{S}_1^j+\frac{1}{2}{Q}_2^i{Q}_2^j+\frac{1}{2}{S}_2^i{S}_2^j\right)=0$

(27)

$E_2\left(X\right)={\mathrm{\Omega }A}^{ki}{S}_1^i+B^{ki}{Q}_1^i+C^{kij}\left(Q_0^i{Q}_1^j+Q_1^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_2^j+\frac{1}{2}{S}_1^i{S}_2^j+\frac{1}{2}{Q}_2^i{Q}_1^j+\frac{1}{2}{S}_2^i{S}_1^j\right)+{\mathcal{F}}_{c1}^k+{\mathcal{M}}_{c1}^k=0$

(28)

$E_3\left(X\right)={-\mathrm{\Omega }A}^{ki}{Q}_1^i+B^{ki}{S}_1^i+C^{kij}\left(Q_0^i{S}_1^j+S_1^i{Q}_0^j+\frac{1}{2}{Q}_1^i{S}_2^j-\frac{1}{2}{S}_1^i{Q}_2^j-\frac{1}{2}{Q}_2^i{S}_1^j+\frac{1}{2}{S}_2^i{Q}_1^j\right)+{\mathcal{F}}_{s1}^k+{\mathcal{M}}_{s1}^k=0$

(29)

$E_4\left(X\right)=2{\mathrm{\Omega }A}^{ki}{S}_2^i+B^{ki}{Q}_2^i+C^{kij}\left(Q_0^i{Q}_2^j+Q_2^i{Q}_0^j+\frac{1}{2}{Q}_1^i{Q}_1^j-\frac{1}{2}{S}_1^i{S}_1^j\right)+{\mathcal{F}}_{c2}^k+{\mathcal{M}}_{c2}^k=0$

(30)

$E_5\left(X\right)=-2{\mathrm{\Omega }A}^{ki}{Q}_2^i+B^{ki}{S}_2^i+C^{kij}\left(Q_0^i{S}_2^j+S_2^i{Q}_0^j+\frac{1}{2}{Q}_1^i{S}_1^j+\frac{1}{2}{S}_1^i{Q}_1^j\right)+{\mathcal{F}}_{s2}^k+{\mathcal{M}}_{s2}^k=0$

(31)

where

$X={\left[\begin{array}{ccccc}{Q}_0& Q_1& S_1& Q_2& S_2\end{array}\right]}^T$

(32)

Now, one can represent equations (27) to (31) in a compact form as

$\mathcal{E}\left(X\right)={\left[\begin{array}{ccccc}{E}_1\left(X\right)& E_2\left(X\right)& E_3\left(X\right)& E_4\left(X\right)& E_5\left(X\right)\end{array}\right]}^T=0$

(33)

The solution of the system of equations presented in equation (33) can be obtained by using the Newton–Raphson method as

$X^{n+1}=X^n-{\mathcal{J}}^{-1}\left(X^n\right)\mathcal{E}\left(X^n\right)$

(34)

where $\mathcal{J}$ is the Jacobian matrix, calculated as

$\mathcal{J}\left(X^n\right)={\frac{\partial \mathcal{E}\left(X\right)}{\partial X}|}_{X=X^n}$

(35)

Using equations (33) and (35), the Jacobian matrix can be represented as

$\mathcal{J}\left(X^n\right)={\left[\begin{array}{c}\frac{\partial E_1}{\partial X}\frac{\partial E_2}{\partial X}\frac{\partial E_3}{\partial X}\frac{\partial E_4}{\partial X}\frac{\partial E_5}{\partial X}\end{array}\right]}_{X=X^n}^T$

(36)

The above relation for the Jacobian matrix produces:

$\mathcal{J}\left(X^n\right)={\left[\begin{array}{c}\frac{\partial E_1}{\partial Q_0}\frac{\partial E_1}{\partial Q_1}\frac{\partial E_1}{\partial S_1}\frac{\partial E_1}{\partial Q_2}\frac{\partial E_1}{\partial S_2}\\ \frac{\partial E_2}{\partial Q_0}\frac{\partial E_2}{\partial Q_1}\frac{\partial E_2}{\partial S_1}\frac{\partial E_2}{\partial Q_2}\frac{\partial E_2}{\partial S_2}\\ \frac{\partial E_3}{\partial Q_0}\frac{\partial E_3}{\partial Q_1}\frac{\partial E_3}{\partial S_1}\frac{\partial E_3}{\partial Q_2}\frac{\partial E_3}{\partial S_2}\\ \frac{\partial E_4}{\partial Q_0}\frac{\partial E_4}{\partial Q_1}\frac{\partial E_4}{\partial S_1}\frac{\partial E_4}{\partial Q_2}\frac{\partial E_4}{\partial S_2}\\ \frac{\partial E_5}{\partial Q_0}\frac{\partial E_5}{\partial Q_1}\frac{\partial E_5}{\partial S_1}\frac{\partial E_5}{\partial Q_2}\frac{\partial E_5}{\partial S_2}\end{array}\right]}_{X=X^n}$

(37)

The different partitions of the Jacobian matrix are derived using equations (27), (28), (29), (30), and (31) and the definition of this matrix as shown in equation (37). These matrices have been presented in Appendix.

2.3. Harmonic balance solution of the linearized system

The proposed HBM can be applied to obtain the linearized form of the intrinsic 1D nonlinear equations. In this case, an analytical, closed-form solution is used to calculate the coefficient vectors of the harmonic terms. The closed-form formulation obtained for the linearized system is then employed as the initial guess for the Newton–Raphson method used for solving the nonlinear equations. Using the linearized equation given in equation (14) and considering the harmonic external loads presented in equation (16), the approximate periodic solution of the linearized system up to the second harmonic terms can be found:

$q_{L_H}^i\left(t\right)=Q_{L1}^i\cos \mathrm{\Omega }t+S_{L1}^i\hspace{0.17em}\sin \mathrm{\Omega }t+Q_{L2}^i\cos 2\mathrm{\Omega }t+S_{L2}^i\hspace{0.17em}\sin 2\mathrm{\Omega }t$

(38)

Substituting equations (16) and (38) into equation (14) yields the following equations:

$A^{ki}\left(-\mathrm{\Omega }{Q}_{L1}^i\hspace{0.17em}\sin \mathrm{\Omega }t+\mathrm{\Omega }{S}_{L1}^i\cos \mathrm{\Omega }t-2\mathrm{\Omega }{Q}_{L2}^i\hspace{0.17em}\sin 2\mathrm{\Omega }t+2\mathrm{\Omega }{S}_{L2}^i\cos 2\mathrm{\Omega }t\right)+{\widehat{B}}^{ki}\left(Q_{L1}^i\cos \mathrm{\Omega }t+S_{L1}^i\hspace{0.17em}\sin \mathrm{\Omega }t+Q_{L2}^i\cos 2\mathrm{\Omega }t+S_{L2}^i\hspace{0.17em}\sin 2\mathrm{\Omega }t\right)+{\mathcal{F}}_{c1}^k\cos \mathrm{\Omega }t+{\mathcal{F}}_{c2}^k\cos 2\mathrm{\Omega }t+{\mathcal{F}}_{s1}^k\sin \mathrm{\Omega }t+{\mathcal{F}}_{s2}^k\sin 2\mathrm{\Omega }t+{\mathcal{M}}_{c1}^k\cos \mathrm{\Omega }t+{\mathcal{M}}_{c2}^k\cos 2\mathrm{\Omega }t+{\mathcal{M}}_{s1}^k\sin \mathrm{\Omega }t+{\mathcal{M}}_{s2}^k\sin 2\mathrm{\Omega }t=0$

(39)

From the above equation, the following relation can be extracted:

$\mathrm{\Omega }{A}^{ki}{S}_{L1}^i+{\widehat{B}}^{ki}{Q}_{L1}^i+{\mathcal{F}}_{c1}^k+{\mathcal{M}}_{c1}^k=0$

(40)

$-\mathrm{\Omega }{A}^{ki}{Q}_{L1}^i+{\widehat{B}}^{ki}{S}_{L1}^i+{\mathcal{F}}_{s1}^k+{\mathcal{M}}_{s1}^k=0$

(41)

$2\mathrm{\Omega }{A}^{ki}{S}_{L2}^i+{\widehat{B}}^{ki}{Q}_{L2}^i+{\mathcal{F}}_{c2}^k+{\mathcal{M}}_{c2}^k=0$

(42)

$-2\mathrm{\Omega }{A}^{ki}{Q}_{L2}^i+{\widehat{B}}^{ki}{S}_{L2}^i+{\mathcal{F}}_{s2}^k+{\mathcal{M}}_{s2}^k=0$

(43)

From equations (40) and (41), one can write the following relation:

$Q_{L1}^i={\left({A^{ki}}^{-1}{\widehat{B}}^{ki}+{\mathrm{\Omega }}^2{{\widehat{B}}^{ki}}^{-1}{A}^{ki}\right)}^{-1}\left(\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{F}}_{s1}^k+\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{M}}_{s1}^k-{A^{ki}}^{-1}{\mathcal{F}}_{c1}^k-{A^{ki}}^{-1}{\mathcal{M}}_{c1}^k\right)$

(44)

$S_{L1}^i={-\left({A^{ki}}^{-1}{\widehat{B}}^{ki}+{\mathrm{\Omega }}^2{{\widehat{B}}^{ki}}^{-1}{A}^{ki}\right)}^{-1}\left(\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{F}}_{c1}^k+\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{M}}_{c1}^k+{A^{ki}}^{-1}{\mathcal{F}}_{s1}^k+{A^{ki}}^{-1}{\mathcal{M}}_{s1}^k\right)$

(45)

In addition, using equations (42) and (43), the following relation for $Q_{L2}^i$ and $S_{L2}^i$ can be obtained:

$Q_{L2}^i={\left({A^{ki}}^{-1}{\widehat{B}}^{ki}+{4\mathrm{\Omega }}^2{{\widehat{B}}^{ki}}^{-1}{A}^{ki}\right)}^{-1}\left(2\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{F}}_{s2}^k+2\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{M}}_{s2}^k-{A^{ki}}^{-1}{\mathcal{F}}_{c2}^k-{A^{ki}}^{-1}{\mathcal{M}}_{c2}^k\right)$

(46)

$S_{L2}^i={-\left({A^{ki}}^{-1}{\widehat{B}}^{ki}+4{\mathrm{\Omega }}^2{{\widehat{B}}^{ki}}^{-1}{A}^{ki}\right)}^{-1}\left(2\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{F}}_{c2}^k+2\mathrm{\Omega }{{\widehat{B}}^{ki}}^{-1}{\mathcal{M}}_{c2}^k+{A^{ki}}^{-1}{\mathcal{F}}_{s2}^k+{A^{ki}}^{-1}{\mathcal{M}}_{s2}^k\right)$

(47)

Using equations (14) and (38), the total response of the linearized system can be written as follows:

$q^i\left(t\right)={q_s+q}_{L_H}^i\left(t\right)$

(48)

3. Results

3.1. Evaluation of the HB approach

In this section, the results of the proposed HBM obtained for solving the nonlinear system given by equation (10) are discussed, and the outputs of the algorithm are compared with a different time solution method for the 1D nonlinear equations. This distinct method is based on the time marching approach, and it uses the central finite difference method to discretize the equations in time domain. It should be noted that the accuracy of the TMM has been evaluated by the author in Siami et al. (2021) by comparing the outputs of the algorithm with the results obtained with the commercial program FLIGHTLAB for a blade.

For the evaluation part, a twisted and curved slender blade with NACA0012 airfoil is used. Figure 2 presents the 3D model of the blade created in ANSYS Workbench™. The geometrical and the material properties of this blade are listed in Table 1. The cross-section of the blade with the generated 1032 linear quadrilateral elements is plotted in Figure 3. It should be stressed that the assigned values for both the initial twist, $k_1$ and the edgewise curvature, $k_3$ are relatively high in comparison with other references (Ho et al., 2010; Yu et al., 2012). This grid has been used in the cross-sectional analysis package developed by the authors to calculate the structural properties required for the 1D analysis. It should be noted that the results of the cross-sectional tool developed by the authors has been evaluated using the results of VABS and the 3D model of a blade in ANSYS. More details about the results of the cross-sectional analysis and the 3D model analysis in ANSYS can be found in another publication by the authors (Siami and Nitzsche, 2023).

Figure 2.

Blade with high initial twist and curvature generated in ANSYS.

Table 1.

Properties of an initially twisted and curved blade.

Geometrical properties
Airfoil	NACA0012
Chord length	0.135 $\left[m\right]$
Length of blade	1 $\left[m\right]$
Initial twist ( $k_1$ )	0.7854 $\left[rad/m\right]$
Initial curvature ( $k_3$ )	0.20 $\left[rad/m\right]$
Material properties
Young’s modulus ( $E$ )	$7.1\times {10}^{10}\hspace{0.17em}\left[Pa\right]$
Shear modulus ( $G$ )	$2.6692\times {10}^{10}\hspace{0.17em}\left[Pa\right]$
Poisson’s ratio ( $\nu$ )	0.33
Density ( $\rho$ )	$2770\hspace{0.17em}\left[kg/m^3\right]$

Figure 3.

Cross-section of the beam (NACA0012 airfoil) with linear quadrilateral elements used for the cross-sectional analysis.

The external harmonic forces applied to the blade are presented in Figure 4. Here, it is assumed that the harmonic forces are applied in two orthogonal directions: flapwise and edgewise directions. The directions of the external loads per unit of length are shown in Figure 4. These forces can be identified as the lift force (flapwise component) and the drag force (edgewise component) applied to the blade. Considering Figure 4 and the length of blade, the maximum values of the harmonic forces are $1058N$ in the flapwise direction and $390N$ in the edgewise direction. The forces are uniformly distributed along the span of the blade; however, different types of load distributions can be considered in the developed HB solution approach. Here, clamped-free boundary conditions are used, as seen in equation (9).

Figure 4.

Applied harmonic forces along the span of the blade in two orthogonal directions (the blue line is in the flapwise direction, and the red line is in the edgewise direction).

The displacement obtained from the HB solution at the tip of blade in z-direction is plotted in Figure 5(a). The calculated values at the tip of blade under the harmonic loads from TMM are presented in this figure to compare the results. It can be observed the calculated displacements from two distinct methods are almost identical. In Figure 5(b), the displacement along the span of blade in z-direction at 10 consecutive time steps have been presented for the HB and TM methods. This figure illustrates a very good correlation between the results of the mentioned approaches in space and time domains. In addition, Figure 5(a) shows that the maximum displacement at the tip of blade is $0.04\hspace{0.17em}m$ , and it illustrate the level of applied force to the blade with $1\hspace{0.17em}m$ length along the span does not cause a high value deflection at the free end of the blade.

Figure 5.

(a) Calculated displacement under the harmonic loads at the tip of blade in z-direction. (b) Calculated displacement under the harmonic loads along the span of blade at 10 consecutive time steps (blue lines are used for the harmonic balance method (HBM) results and black lines for the time marching method (TMM) results).

In Figure 6, the obtained linear velocities from the HBM and TMM are plotted in the y- and z-directions. In addition, the comparisons among the results obtained for forces and moments at the root of blade from the two approaches are presented in Figure 7 and Figure 8. According to these figures, it can be seen that the results of proposed method agree very well with those obtained with the time marching method.

Figure 6.

(a) Calculated velocity in y-direction under the harmonic loads at the tip of blade. (b) Calculated velocity in z-direction under the harmonic loads at the tip of blade (blue lines are used for the HBM results and black lines for the TMM results).

Figure 7.

(a) Calculated force in y-direction under the harmonic loads at the root of blade. (b) Calculated force in y-direction under the harmonic loads at the root of blade (blue lines are used for the HBM results and black lines for the TMM results).

Figure 8.

(a) Calculated moment in y-direction under the harmonic loads at the root of blade. (b) Calculated moment in z-direction under the harmonic loads at the root of blade (blue lines are used for the HBM results and black lines for the TMM results).

3.2. Effect of linearization of 1D nonlinear equations

In this section, the responses of the intrinsic 1D nonlinear equations are compared with the linearized form given by equation (14) to evaluate the capability of using this analytical formulation for time domain analysis of blades. As it has been developed in Section 2, for the linearized system, the analytical HB solution is found using equations (44), (45), (46), (47), and (48). As such, the solution of this system is completely independent from numerical algorithms and produces a very low-cost solution for the 1D intrinsic equations. In the proposed scheme for the HBM solution of the fully nonlinear system, the Newton method is employed to find the coefficient vectors of the HBM based on the developed analytical Jacobian matrix. The linearized solution is applied as the initial guess of the numerical solution.

The results are calculated for a blade with the same cross-sectional properties given in Table 1 and same boundary conditions depicted in equation (9). The length of blade for the 1D analysis is 2 $m$ . In addition, to highlight the effect of external loads, two different levels of harmonic loads are applied to the blade as shown in Figure 9. In these plots, the second case presents a condition in which the external harmonic forces in both orthogonal directions (i.e. the flapwise and edgewise directions) are 10 times higher than those of the first case. The frequency of the harmonic load is 148 $rad/s$ and includes a combination of sine and cosine terms.

Figure 9.

Applied harmonic forces along the span of the blade in two distinct cases: (a) applied forces in flapwise direction, (b) applied forces in edgewise direction (blue lines are the load components for Case 1 and the black lines are those associated to Case 2).

The comparison between the HBM solution of the linearized equations with that obtained for the corresponding nonlinear system is presented in Figure 10 for the lower level of excitation loads. This figure shows the velocity versus displacement of the blade at the free end in the flapwise and edgewise directions. Based on the presented graphs, it can be observed that both the analytical solution of the linearized system and the corresponding solution of the nonlinear system can well predict the dynamic behavior of the structure. However, under the higher level of excitation loads, the differences between the two solution approaches appear as indicated in Figure 11. Based on the velocity-displacement plots at the blade free end, the HB solution of the linearized equations fails to follow the HB solution of the nonlinear system. In other words, Figure 11 shows by increasing the values of applied forces, the nonlinear behavior of the blade is appeared, and it demonstrates the importance of using the solution of nonlinear system of equations compared to the linearized form of equations to predict the dynamic behavior of the blade. It also should be noted that the TMM method when the higher level of load was present could not provide acceptable results. In addition, it required for the fully nonlinear equations a very fine time domain resolution, therefore demanding much higher simulation times when compared to the proposed HBM.

Figure 10.

Velocity versus displacement under lower-level loads (Case 1) for the linearized and fully nonlinear systems: (a) y-direction, (b) z-direction.

Figure 11.

Velocity versus displacement under higher-level loads (Case 2) for the linearized and fully nonlinear systems: (a) y-direction, (b) z-direction.

4. Conclusion

In this work, the HBM combined with the Galerkin approach has been applied to time domain analysis of the intrinsic 1D nonlinear equations. First, the system of equations is discretized spatially by the Galerkin approach. Then, the first and second harmonic components are used to obtain the time response of the structure under the applied harmonic excitations. The harmonic balance method (HBM) solutions have been separately presented for both the fully nonlinear and the linearized form of the equations. For the numerical solution used in the nonlinear case, the Jacobian matrix has been derived analytically to increase the time efficiency of the method and the linearized solution was used as the initial guess. The accuracy of the developed approach was evaluated by comparing its results with the results of a time marching solution for an initially twisted and curved blade. In addition, the outputs of the HBM linearized form of equations are compared with the fully nonlinear counterpart under different level of external loads. The developed analytical HBM solution for the linearized system provides very similar results to the HBM solution of the nonlinear system under the lower level of harmonic external loads. However, it diverged from the results of the nonlinear solution under the higher level of loads as demonstrated by the velocity-displacement graphs of two load cases. By increasing the level of applied loads on the blade, the effect of nonlinearity described by the 1D equations becomes more significant, and the analytical solution provided for its linearized form loses accuracy. However, the analytical solution can be used as an initial guess for the fully nonlinear problem.

The proposed method using the Galerkin approach and the HBM approach can provide a powerful tool for the structural dynamic analysis of twisted/curved anisotropic slender beams such as blades. The algorithm is very time efficient in comparison to other common tools for analyzing these types of structures. The accuracy of the proposed formulation has been addressed and evaluated for different levels of external loads.