# Cyclic Symmetry Finite Element Forced Response Analysis of a

## Transcript Of Cyclic Symmetry Finite Element Forced Response Analysis of a

Cyclic Symmetry Finite Element Forced Response Analysis of a Distortion-Tolerant Fan with Boundary Layer Ingestion

J.B. Min, T.S.R. Reddy†, M.A. Bakhle

and

R.M. Coroneos, G.L. Stefko A.J. Provenza, K.P. Duffy†

NASA Glenn Research Center, Cleveland, Ohio †University of Toledo, Toledo, Ohio

Accurate prediction of the blade vibration stress is required to determine overall durability of fan blade design under Boundary Layer Ingestion (BLI) distorted flow environments. Traditional single blade modeling technique is incapable of representing accurate modeling for the entire rotor blade system subject to complex dynamic loading behaviors and vibrations in distorted flow conditions. A particular objective of our work was to develop a high-fidelity full-rotor aeromechanics analysis capability for a system subjected to a distorted inlet flow by applying cyclic symmetry finite element modeling methodology. This reduction modeling method allows computationally very efficient analysis using a small periodic section of the full rotor blade system. Experimental testing by the use of the 8-foot by 6-foot Supersonic Wind Tunnel Test facility at NASA Glenn Research Center was also carried out for the system designated as the Boundary Layer Ingesting Inlet/Distortion-Tolerant Fan (BLI2DTF) technology development. The results obtained from the present numerical modeling technique were evaluated with those of the wind tunnel experimental test, toward establishing a computationally efficient aeromechanics analysis modeling tool facilitating for analyses of the full rotor blade systems subjected to a distorted inlet flow conditions. Fairly good correlations were achieved hence our computational modeling techniques were fully demonstrated. The analysis result showed that the safety margin requirement set in the BLI2DTF fan blade design provided a sufficient margin with respect to the operating speed range.

I. Introduction

THE concept of relocating the engines from a podded position on the wing or fuselage and embedding them within the airframe [1] (Figure 1) can improve vehicle fuel burn, weight, noise, and emissions (through reduced fuel burn and engine cycle improvements). This Hybrid-Wing-Body (HWB) architecture is particularly suited to embedded engines, as balance requirements already place them near the aft of the airframe.[2] The embedded engines on a HWB can be integrated into the fuselage and ducted to exhaust through the rear of the vehicle. A main feature of this proposed boundary layer ingested inlet distortion-tolerant fan aircraft is the embedded aft-mounted engines which ingest part of the fuselage boundary layer air flow to reduce fuel burn. This embedded aft engines may also provide reduced susceptibility to bird strike since engines are partly hidden head-on especially at take-off, and possibly enabling to reduce the engine noise to the ground from positioning the engines above the airframe. The National Aeronautics and Space Administration (NASA) and United Technologies Research Center (UTRC) with contributions from Virginia Polytechnic and State University (Virginia Tech) have developed a BLI propulsor to significantly reduce the fuel consumption. Accordingly, a high-efficiency Boundary Layer Ingestion (BLI) flow Distortion-Tolerant Fan (DTF) blade rotor system has been studied as part of the Boundary Layer Ingesting Inlet/Distortion-Tolerant Fan (BLI2DTF) Project.

1

Inlet airflow distortion ingestion into a fan gives rise to additional aero-mechanical vibration and a decrease in stability margin. An embedded BLI propulsion system inherently involves a persistent and severe inlet distortion reaching the fan at most operating conditions, resulting in high dynamic stresses due to aeroelastic forced responses and the possibility of flutter. If not adequately addressed by the propulsor design, these aeromechanics phenomena could cause fan structural failure due to high cycle fatigue (HCF) or unstable self-excited failure, respectively. Numerical simulation modelling for this dynamically interconnected problem is a challenging problem.

Figure 1: Conceptual aircraft design and boundary layer ingestion inlet distortion-tolerant rotor-fan system (inset at the bottom: BLI2DTF fan in the 8’x6’ Supersonic Wind Tunnel test at NASA GRC)

Aero-mechanical vibration stresses in fan blades due to inlet airflow distortion is of particular interest in this study since the increase of stress in a blade can lead to its premature failure. Blade failure is usually occurring within short time by unstable self-excited (flutter) overload and high cycle fatigue (HCF) dynamic load. Turbomachinery blade vibration stress problem has been reported as the single cause to termination in new engine development effort. “While over 90 percent of the potential HCF problems are uncovered during development testing of a new engine, the remaining few account for nearly 30 percent of the total development cost and are responsible for over 25 percent of all engine distress events.”[3]

One of the main research objectives of the BLI2DTF Project was to develop and apply aeromechanics analysis tools to predict the effects of unique BLI distortion on “designed to be distortion-tolerant” fan system.[4,5] Traditional single blade modeling techniques are incapable of accurately modeling the entire rotor blade system due to complex dynamic loading behaviors and vibrations in distorted flow conditions. Thus, the present aeromechanics analysis on BLI2DTF fan rotor system had four specific objectives: (1) establish a model based approach that enables a highfidelity full-rotor aeromechanics analysis capability for performance, safety, and computational efficiency in design analysis of a fan rotor system subjected to a distorted inlet flow, (2) complete forced response analyses to determine dynamic stresses in DTF fan blade due to continual operation in a distorted flow resulting from BLI using cyclic symmetry modeling techniques, (3) assess the concerns of fatigue problem using Goodman Diagram approach by examining predicted stress values regarding the HCF behavior due to blade static and vibration dynamic stresses, and (4) provide guidance on the strain measurement by determining optimal strain gauge locations on the blades towards a cost-effective experimental test measurement from overall and/or directional strain characteristics predicted due to a loading condition aforementioned. This latest high-fidelity modeling capability enables the analysis of arbitrary inlet distortion patterns that can be specified through variations of flow properties, together with prescribed characteristics of DTF fan blade vibrations. The present analysis approach will ensure operability of DTF fan blades in the embedded propulsion system of HWB aircraft.

II. Cyclic symmetry finite element modeling technique

Understanding complex mechanical and aerodynamic interaction of turbomachine blading in this distorted flow environment has become topmost to the success of BLI2DTF fan design and durability. As mentioned above, typical single blade analysis models widely employed in the turbomachinery research community are incapable of representing the design for both understanding the harmonic modal blade dynamic behaviors and accurately predicting

2

corresponding forced response information of the blades subjected to complex distorted flow dynamic loading and vibrations. Nevertheless, analyzing this type of the problems by modeling a full 360 degree rotor blades structural system is not computationally efficient and practical, especially in large size systems. With addressing this issue in the present DTF fan analysis, reduction modeling techniques were essential and a novel reduction modeling capability was developed for fan distortion analysis to model non-axisymmetric large turbomachinery BLI2DTF rotor blades system. To achieve this capability, cyclic symmetry analysis modeling formulations were incorporated with the Ansys finite element (FE) computer program[6] to analyze the BLI2DTF forced responses.

Cyclic symmetry is a type of symmetry in which identical sectors are attached in a circular pattern about a rotation axis. It is assumed that in each individual sector’s local coordinate system the stiffness and mass matrices for all of the sectors are identical. The reduction of degrees-of-freedom is accomplished at the sector level. The compatibility of sector interfaces generates a reduced set of equations of the complete structure. A key element of cyclic symmetry method is to determine the cyclic symmetry transformation. That is, the solution of the complete structure can be reduced to the solution of the cyclic symmetric substructures. Once the solution for the substructure is determined, the solution of the complete structure is given by the cyclic symmetry transformation. The general procedure for the implementation of this method is given in reference 7.

The sector geometry of a BLI2DTF fan rotor system was meshed[8] with higher order 20-node hexahedral elements for the airfoil region. Higher order 10-node tetrahedral elements for the airfoil platform and the disk part were more densely meshed for better prediction of stress results. This sector model represents a cyclic symmetric sector of a whole DTF fan rotor system, shown in Figure 2.

(a)

(b)

(c)

Figure 2: a) FE mesh of a cyclically symmetric “base sector” model of BLI2DTF (close-up), b) description of cyclic symmetry interface boundaries, c) fully expanded mesh by cyclic symmetry formulation

Eighteen blades were placed symmetrically around the disk, so a 20 degree sector was used to produce a symmetric sector. The airfoil itself was meshed with 4 hexahedral elements through the thickness in order to properly capture the bending stiffness. The resulting mesh of a sector geometry having 557,554 elements and 992,548 nodes is shown in Figure 2 which requires to solve 5,912,580 equations for a particular engine operating condition. It was anticipated to be able to analyze this extremely large model by solving the equations only for a smaller sized sector model without losing the computational accuracy rather than analyzing the entire 360 degree rotor blade model requiring unrealistic computing time and hardware capacity.

The cyclic sector model has two edges that align along the surfaces of cyclic symmetry (Fig. 2, (a)) in angles. The edge having the algebraically lower angle, θ, in the cylindrical coordinate system is called the low edge and the one having the higher angle, θ + α, is called the high edge. The angle, α, between the two successive surfaces of cyclic symmetry is called the sector angle. The architecture of the cyclic symmetry solution process depends upon how the compatibility and equilibrium conditions of the cyclic sector are enforced in the matrix-solution process using the duplicate sector method. [7]

The principal features of this cyclic symmetry modeling method are briefly described below. Along with the notations illustrated in Figure 2(b), all loading, boundary conditions, and coupling and constraint equations being

3

present on the basic sector (A) are applied to the duplicate sector (B). Constraint relationships (Eqn. (1)) can be defined to relate the lower (θ = 0) and higher (θ = α) angle edges of the basic sector to allow calculations related to a given number of harmonic indices (k) defined below. Essentially, the basic sector is duplicated in the analysis steps to satisfy the required constraint relationships and to obtain nodal displacements.

{𝑢𝐴′ } = [ 𝑐𝑜𝑠𝑘𝛼 𝑠𝑖𝑛𝑘𝛼 ] {𝑢𝐴}

(1)

𝑢𝐵′ −𝑠𝑖𝑛𝑘𝛼 𝑐𝑜𝑠𝑘𝛼 𝑢𝐵

where:

𝑢𝐴, 𝑢𝐵 = calculated displacement on lower angle side of basic and duplicated sectors (A and B, respectively) 𝑢𝐴′ , 𝑢𝐵′ = displacements on higher angle side of basic and duplicated sectors (A and B, respectively) determined from constraint relationships

𝑁 𝑖𝑓 𝑁 𝑖𝑠 𝑒𝑣𝑒𝑛 k = harmonic index = 0, 1, 2…. {𝑁2−1

2 𝑖𝑓 𝑁 𝑖𝑠 𝑜𝑑𝑑 𝛼 = 2𝜋⁄𝑁 = sector angle N = total number of sectors in 360 degree (N = 18)

The cyclic symmetric solution sequences consist of three basic steps. The first step transforms applied loads to cyclic symmetric components using finite Fourier theory and enforces cyclic symmetry constraint equations (Eqn. (1)) for each harmonic index (i.e. nodal diameter) (k = 0, 1, . . ., N/2). At the solution stage of a cyclic symmetry analysis, the cyclic symmetry compatibility conditions are enforced for each harmonic index solution via coupling and/or constraint equations (CEs) connecting the nodes on the low- and high-edge components on the basic and duplicate sectors. As stated earlier, the cyclic symmetry solution process depends upon how the compatibility and equilibrium conditions of the cyclic sector are enforced in the matrix-solution process. The solution stage performs each harmonic index solution as a separate load step where the real (basic sector) and imaginary (duplicate sector) parts of the solution are calculated, as expressed below in Eqn. (2) – Eqn. (5). The Fourier transformation from physical components, u, to the different harmonic index components, 𝑢̅, is given by,

For harmonic index, k = 0 (symmetric mode), 𝑢̅𝑘=0 = 𝑁1 ∑𝑁𝑛=1 𝑢𝑛 (2)

For harmonic index, 0 < k < N/2 (mix mode),

basic sector: (𝑢̅𝑘)𝐴 = 𝑁2 ∑𝑁𝑛=1 𝑢𝑛 cos(𝑛 − 1) 𝑘𝛼 (3)

duplicate sector: (𝑢̅𝑘)𝐵 = 𝑁2 ∑𝑁𝑛=1 𝑢𝑛 sin(𝑛 − 1) 𝑘𝛼 (4)

For N even only, k = N/2 (anti-symmetric mode): 𝑢̅𝑘=𝑁/2 = 𝑁1 ∑𝑁𝑛=1(−1)(𝑛−1)𝑢𝑛 (5)

where: u = displacements (or any physical component, such as, forces, pressure loads, temperatures, and inertial loads 𝑢̅ = cyclic displacements (or any cyclic symmetric component) n = sector number, varies from 1 to N

In the second step, any applied load on the full 360 degree fan model is treated through a Fourier transformation process above and applied on to the cyclic basic sector. For each value of harmonic index, k, the procedure solves the corresponding equation. In the third step, the responses are then combined to get the complete response of the full fan structure via the Fourier expansion given by Eqn. (6). That is, using the cyclic symmetry solution of the basic and duplicate sectors connecting low and high edges of basic and duplicate sectors combines the solutions from the two sectors by performing computations on the selected load step to combine the results of the two sectors. Note that only a harmonic index zero solution (k=0) is valid for a static solution with cyclic loading.

4

The transformation to physical components, u, from the cyclic symmetry, 𝑢̅ , components is recovered by this equation,

𝑢𝑛= ̅𝑢𝑘=0 + ∑𝑁𝑘=1[ (𝑢̅𝑘)𝐴 cos(𝑛 − 1)𝑘𝛼 + (𝑢̅𝑘)𝐵sin(𝑛 − 1)𝑘𝛼] + (−1)(𝑛−1) ̅𝑢𝑘=𝑁/2

(6)

The last term (−1)(𝑛−1) ̅𝑢𝑘=𝑁/2 exists only for a case of N even. In cyclic symmetry modal analysis, the mode shape in each sector is obtained from the eigenvalue solution. The

cyclic symmetry modal solution occurs after the static cyclic symmetry solution, so that the modal analysis yields the mode shapes of the fan blade in the pre-stressed state. The displacement components (x, y, or z) at any node in sector n for harmonic index k, in the full fan structure are yielded by Eqn. (6).

If the mode shapes are normalized to the mass matrix in the modal analysis, the normalized displacement components in the full fan blade system are given by,

𝑢 √𝑁

𝑖𝑓 𝑘 = 0 𝑜𝑟 𝑘 = 𝑁⁄2

𝑢𝑛𝑜𝑟𝑚𝑎𝑙𝑖𝑧𝑒𝑑 = { 𝑢

𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

(7)

√𝑁⁄2

In a cyclic symmetry modal analysis, it is necessary to understand the concepts of harmonic indices and nodal diameters. Examining at a corresponding point in each sector, the mode shapes are sinusoidal in the circumferential direction, which indicates nodal lines on the bladed disk system called as nodal diameters (ND). Most mode shapes contain lines of zero out-of-plane displacement which cross the entire disk. This modal behavior is described in more detail in Section III.

In addition to cyclic symmetry constraint equations defined above, surface-to-surface contact pairs were used to define contact behaviors at the blade disk dovetail assembly as illustrated in Figure 3. In the contact between two bodies, the surface of one body was taken as a contact surface and the surface of the other body as a target surface. The “contact-target” pair concept has been widely used in finite element simulations. The contact surface was meshed using contact elements which are a 3D, 8-node, higher order quadrilateral element located on the surfaces of blade dovetail area. The target surface was meshed using target elements which are a 3D, 8-node, higher order quadrilateral to represent 3D target surface for the associated contact elements on the surfaces of disk dovetail socket area. For these surface-to-surface contact elements, the augmented Lagrangian penalty method,[6,9] incorporated with a contact "spring" stiffness called the contact stiffness and target stiffness was used to establish a relationship between the contact-target surfaces for bonded contact in conjunction with cyclic symmetry constraint conditions. Contact problems are highly nonlinear and require significant computer resources. In the present analysis, it was assumed that both contacting bodies are deformable with similar stiffness so flexible-to-flexible contact type was employed. Generally the regions of contact are not known until an actual solution process is evolved with the loads, materials, boundary conditions, and other factors in a largely unpredictable manner.

contact elements on surface in contact with

disk

Figure 3: Surface-to-surface contact pairs between disk hub and blade in dovetail attachment

The augmented Lagrangian method used in this study is an iterative series of penalty updates to find the Lagrange multipliers (i.e., contact tractions) and was investigated with the magnitude of the contact stiffness coefficient from a common range of 0.01 - 1.0. It was observed from this investigation that the augmented Lagrangian penalty method was insensitive to the magnitude of the contact stiffness coefficient. Hence, this method was employed in the present

5

DTF fan forced response analysis, presumably well-suited for a wide range of frictional contact situations governed by a Coulomb law[9] at variable rotor speeds of DTF fan blade. Further boundary conditions are applied at the inner radius of the disk sector to restrain motion in the radial and circumferential directions to prevent the rigid body motions. Also, surface elements manipulated by thickness set to zero were overlaid onto an area face of the blade elements as a modeling technique to apply steady and unsteady aerodynamic pressure loads on the blade surfaces.

In materials used for this analysis, the blade part was made from a titanium alloy and the disk part was made from a steel.

III. Methods for forced response computation

A finite element forced response analysis method utilizing the cyclic symmetry modeling approach was employed to provide a comprehensive study of a variety of structural characteristics of a rotating DTF fan blade system. The first task (section A) was to examine the steady stresses including nonlinear rotational effects of spin and stress stiffening due to the centrifugal loading caused by rotational velocity and steady state aerodynamic blade pressure load. The second task (section B) was to perform a prestressed vibration modal analysis, allowing for the classification of every natural mode within the frequency range of interest. This computation also provided information about areas of high stress and how the relative magnitude of these stresses varies with frequency. Finally, section C, a cyclic symmetry finite element analysis model incorporated with Turbo aerodynamic analysis computer program[10] and Ansys structural dynamics computer program was developed to compute the forced response stress distribution in fan blade and as a result the blade structural safety margins for the targeted DTF blade design requirements are presented in section D. While all the equations in the following sections are represented in the cyclic symmetry coordinates, added symbols to the equations are avoided for a brevity in the descriptions.

A. CYCLIC SYMMETRY CENTRIFUGAL STATIC STRESS ANALYSIS MODEL SUBJECTED TO STEADY-STATE AERODYNAMIC PRESSURE AT FAN SPEEDS

The first task was to develop a cyclic symmetry finite element static analysis model incorporated with Turbo steady-state aerodynamic blade pressure load as a function of the fan blade rotation. To capture the true blade behavior in this rotational condition, geometric nonlinearity was to be accounted for providing the blade deflections potentially to cause significant changes in the blade structural geometry. In such events, while solving the equations of equilibrium by utilizing an iterative process would be employed to obtain the accurate prediction for the largely deformed configuration, excessively high computing time is required for the size of the present problem. Even if a small deflection analysis may not directly account for change in geometry, the effect for large deflection would be possible to be accounted for by an adjustment of the structural stiffness matrix in a small deflection solution approach. This approach allows adjusting the stiffness of a rotating blade to account for dynamic mass effects. Equilibrium of the system and centrifugal forces on the mass using small deflection logic can be expressed as,

[𝐾]{𝑢} = {𝛺2}[𝑀]{𝑟}

(8)

where:

[K] = global stiffness matrix {𝑢} = displacement of the mass from the rest position {𝑟} = rest position of the mass with respect to the axis of rotation 𝛺 = angular velocity of rotation [𝑀] = global mass matrix = ∑𝑒𝑒==𝑛1 𝑀𝑒 e = element 1, 2, , n [𝑀𝑒] = element mass matrix = ∫𝑉[𝑁]𝑇[𝑁]𝜌𝑑𝑉 [𝑁] = shape function matrix

ρ= element density

n= number of elements

However, to account for deflection effects, Eqn. (8) must be expanded to:

[𝐾]{𝑢} = {𝛺2}[𝑀]{𝑟 + 𝑢}

(9)

Rearranging terms,

([𝐾] − {𝛺2}[𝑀]){𝑢} = {𝛺2}[𝑀]{𝑟}

(10)

6

Defining:

[𝐾̅] = [𝐾] − {𝛺2}[𝑀]

(11)

[𝐹̅] = {𝛺2}[𝑀]{𝑟}

Eqn. (10) becomes simply,

[𝐾̅]{𝑢} = [𝐹̅]

(12)

[𝐾̅] is the stiffness needed in a small deflection solution to account for large deflection effect, and [𝐹̅] is the same as that derived from small deflection logic. Rotational velocities are combined with the element mass matrices to form a body force load vector. Thus, the large deflection effects are included in a small deflection solution. This decrease in the effective stiffness matrix is called spin-softening.[11] Additionally, with assuming that the latest stress state in the computation may be large enough to change the blade shape under high centrifugal rotation, inclusion of such stiffness effects may need to be examined for the potential changing geometry due to this deformation by large stress which may no longer be neglected. This effect called stress stiffening effect was also included in the computation. This stress stiffening is the stiffening (or weakening) of a structure due to its stress state. This stiffening effect normally needs to be considered for thin structures with bending stiffness very small compared to axial stiffness, such as thin DTF blades which couples the in-plane and out-of-plane displacements. Accordingly, this stress stiffness matrix was added to the stiffness matrix in Eqn. (12) in order to update the total stiffness for accurate predictions. Hence, the equilibrium equation for static analysis is given by,

([𝐾̅] + [𝑆̅]){𝑢} = [𝐹̅]

(13)

where:

{𝑢} = nodal displacement vector [𝐾̅] = global stiffness including spin softening stiffness due to centrifugal force [𝑆̅] = stress stiffening (or weakening) due to its stress state [𝑆̅] = ∫ [𝐺]𝑇[𝜏] [𝐺]𝑑𝑉 when necessary

𝑉

[𝐺] = matrix shape function derivatives

[𝜏]= matrix of current Cauchy (true) stresses in global coordinate [𝐹̅] = global load vector due to rotating body force and steady state aerodynamic pressure load

As for the mass term, the mass summary is calculated from accumulating each element contribution which does not reflect the boundary conditions, coupled degrees of freedom, constraint equations, and multipoint constraint in contact elements.

B. PRESTRESSED CYCLIC SYMMETRY VIBRATION MODAL ANALYSIS MODEL SUBJECTED TO CENTRIFUGAL STRESS WITH STEADY-STATE AERODYNAMIC PRESSURE AT FAN SPEEDS

Given that the cyclic symmetry harmonic forced responses are to be computed in terms of the unsteady state aerodynamic engine order (EO) excitation vibration loads, a perturbed prestressed cyclic symmetry modal analysis model was developed after first performing a cyclic symmetry static analysis for identifying frequencies and mode shapes of blades vibration and critical speed ranges for a given set of nodal diameters (as illustrated in Figure 4). This allows for the classification of every natural mode within the frequency range of interest.

This vibration modal analysis assumed that there is no damping and the structure has no time varying forces applied (i.e., free vibration) to determine the cyclic symmetry natural frequencies and mode shapes. The equation of motion for an undamped system can be expressed in matrix notation,

[𝑀]{𝑢̈ (𝑡)} + [𝐾̃]{𝑢(𝑡)} = {0}

(14)

where: [𝐾̃] = [𝐾̅] + [𝑆̅]

For a linear system, free vibrations are harmonic of the form,

{𝑢(𝑡)} = {𝑢} = {𝜙}𝑟𝑐𝑜𝑠𝜔𝑟𝑡

(15)

7

where: {𝜙}𝑟 = eigenvector representing the mode shape of the rth natural frequency 𝜔𝑟 = rth natural frequency (cycles/second, Hz); 𝜔𝑟 = 2𝑓𝑟𝜋 where 𝑓𝑟 = natural circular frequencies (radians/second)

0 ND

1 ND

2 ND

3 ND

4 ND

…

9 ND

Note: ND = Nodal Diameter; Mode refers to Mode Family; Magnitude of displacement is plotted

Figure 4: Illustration of nodal diameter (ND) patterns from BLI2DTF model shown with higher deformation in red color

As indicated in Fig. 4, the cyclic symmetry prestressed modal analysis solution provides analytical results for the interaction between the disk vibration and the blade vibration of DTF fan blade. The blade frequencies were solved corresponding to a particular disk mode (or nodal diameter, ND). Different disk nodal diameter information is to capture the blade frequencies and mode shapes for each nodal diameter.

The cyclic symmetry modal solution uses the same low- and high-edge components defined in the static cyclic symmetry analysis. To calculate the frequencies and mode shapes of a deformed structure, the linear perturbation solution method based on the Newton-Raphson procedure[12] after performing a static analysis is needed for a prestressed modal analysis of cyclic symmetry structures. The linear perturbation modal analysis procedure is designed to solve a problem from this preloaded condition. The linear perturbation modal analysis can be viewed as an iteration solution on top of a prior analysis. During the linear perturbation process, all of the effects from the static analysis are taken into account and are “snapshot solution” so that the perturbed loads can be generated linearly by using the "snapshot solution” matrices and material properties. With a restart procedure from the static analysis, the iterative process is used for a restart solution by allowing modifications of perturbed loads to generate new matrices of perturbation loads, mass density, and damping. The load history vector from the modal analysis is updated for the next iteration solution and subsequently the modal coordinates are updated by solving eigenvalue solution. During this phase of the analysis, the frequency dependent perturbation load vector is calculated and stored in the mode shape file and the element load information files for a subsequent mode-superposition modal-based forced response analysis. This procedure allows for saving significant time in computation and reusing the mode shape information from an earlier modal analysis to improve the accuracy of a mode-superposition harmonic forced response computation. Modal extracting computation process typically requires more time than element loads generation, perturbation residual load vector calculation, and enforced snapshot solution mode calculation.

Hence, the equation of motion for an undamped system can be expressed for the linear perturbation vibration modal analysis to compute the natural frequencies and mode shapes,

8

(−𝜔𝑟2[𝑀] + [𝐾̃𝑙]){𝜙}𝑟 = {0}

(16)

where: [𝐾̃𝑙] = stiffness update at the iteration of perturbation load step “l” for prestressed modal analysis where spin-

softening effects and stress stiffening effects included

This equality is satisfied if either {𝜙}𝑟 = {0} or if the determinant of (−𝜔𝑟2[𝑀] + [𝐾̃𝑙]) is zero. The first option is the trivial one and, therefore, is not of interest. Thus, the second one gives the solution,

|[𝐾̃𝑙] − 𝜔𝑟2[𝑀]| = 0

(17)

This is an eigenvalue problem which may be solved for up to ndof values of 𝜔𝑟2 and ndof eigenvectors {𝜙}𝑟 which satisfy Eqn. (16) where ndof is the number of degrees-of-freedom (DOFs) of the model.

While several eigenvalue and eigenvector extraction procedures are available for solving Eqn. (17), the present study employed the Block Lanczos method [13] to extract the requested eigenvalues using the sparse direct solving technique that performs forward and backward substitution using the factors. During this solution process, the frequency dependent perturbation load is calculated and stored as nodal forces for a subsequent mode-superposition modal-based harmonic forced response analysis.

Repeat analysis was performed to calculate blade vibration frequencies at different rotational speeds. Also, as will be described in detail in Section IV, an additional objective of this prestressed cyclic symmetry modal analysis was to guide the optimal strain gauge placement on the fan blades by providing the information regarding how the magnitude of modal stresses varies with frequencies so that the optimal locations of the strain gauges are determined on a blade for the effective strain gauge measurement in the experimental wind tunnel test of fan blades.

C. CYCLIC SYMMETRY HARMONIC FORCED RESPONSE ANALYSIS MODEL SUBJECTED TO A COMPLETE CYCLE OF ENGINE ORDER FORCING VIBRATION BY UNSTEADY-STATE AERODYNAMIC PRESSURE LOAD AT FAN SPEEDS

Next, a cyclic symmetry harmonic forced response analysis model was developed utilizing the mode superposition method. The results of a cyclic symmetry modal analysis described in Section III.B was retrieved for a harmonic response analysis described in this section. This cyclic symmetry harmonic forced response analysis model allows for a complete solution of DTF fan blade vibration dynamic stresses throughout a complete cycle of engine order forcing vibration.

The finite element semi-discrete equation of motion for rotation including prestressed structural stiffening effects can be expressed in matrix form,

[𝑀]{𝑢̈ (𝑡)} + [𝐶]{𝑢̇ (𝑡)} + [𝐾̃𝑙]{𝑢(𝑡)} = {𝐹𝑎(𝑡)}

(18)

where: [𝐶] = structural damping matrix {𝑢̈ (𝑡)} = nodal acceleration vector {𝑢̇ (𝑡)} = nodal velocity vector {𝑢(𝑡)} = nodal displacement vector {𝐹𝑎(𝑡)} = unsteady aerodynamic EO excitation pressure load

t = time representing the current equilibrium iteration

l = load step in perturbation solutions

The cyclic symmetry harmonic forced response analysis solves the time-dependent equations of motion (Eqn. (18)) for DTF fan blade structures undergoing harmonic forced vibration. This analysis assumed that the entire DTF blade structure has frequency-dependent stiffness, damping, and mass effects and all loads and displacements vary sinusoidally at the same known frequency (although not necessarily in phase). Mathematically, Eqn. (18) represents a system of linear differential equations of second order. Two procedures can be considered for the solution of Eqn. (18): direct integration method and mode superposition method. Given that a traveling wave EO excitation occurs due to circumferential disturbances in the BLI flow field from upstream of the fan and a mode superposition method shows the numerical effectiveness over a direct integration method, engine order loading (circumferentially traveling wave excitation) is modeled with a mode-superposition analysis method.[6,14] Engine order excitation is periodic

9

forcing arising from circumferential disturbance in a set of blades on the bladed rotor. The mode-superposition method sums factored mode shapes using the natural frequencies and mode shapes obtained from a linear perturbed modal analysis to characterize the dynamic response of a harmonically excited structure. In addition, since a mode superposition method may need to consider only a few modes, the mode superposition method can be much effective than a direct integration which would be prohibitively expensive in computations if the order of the matrices is large.[14]

Eqn. (18) was solved by the mode-superposition method to the interest of specific EO exciting unsteady aerodynamic pressure load on the blades. {𝐹𝑎} is the time-varying load vector given by,

{𝐹𝑎} = {𝐹𝑛𝑑} + 𝑠𝑓. {𝐹𝑠}

(19)

where: {𝐹𝑛𝑑} = time varying unsteady aerodynamic nodal forces

sf = load vector scale factor {𝐹𝑠} = load vector read in from the linear perturbation modal analysis in Section III.B

Note that all loads from the linear perturbation modal analysis are scaled, including forces and pressures. All loads from the modal solution load vector during a given load step in the linear perturbation modal analysis is applied as engine order loads. Damping in Eqn. (18) should be specified; otherwise, the response will be infinite at the resonant frequencies. Additional information on the transformation of Eqn. (18) from cyclic symmetry coordinates to physical coordinates and Fourier transformation in cyclic symmetry forced response analysis are described in reference 6.

Hence, similarly given in reference 14, a general treatment of dynamic systems can be formulated by Lagrange’s equation from the scalar quantities of kinetic energy, potential energy, and work where the equations of motion of a system can be formulated in a number of independent coordinate systems. Such independent coordinates are called generalized coordinates and can be denoted by the character 𝑞𝑟 in the present study. Eqn. (18) can be transferred into modal space and hence, a set of modal coordinates 𝑞𝑟was defined such that,

{𝑢} = ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟

(20)

where: {𝜙𝑟} = rth mode shape nmod = number of modes to be used from the preceding modal analysis in Section III.B

Substituting Eq. (20) into Eqn. (18)

[𝑀] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟}𝑞𝑟̈ + [𝐶] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟̇ + [𝐾̃𝑙] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟 = {𝐹𝑎(𝑡)}

(21)

After few steps of the matrix transformation and substitution with applying its orthogonality and normality conditions, and only the r = j terms remain, Eqn. (21) can be expressed in the modal coordinates as,

𝑞𝑗̈ + 2𝜔𝑗𝜉𝑗𝑞𝑗̇ + 𝜔𝑗2𝑞𝑗 = 𝑓𝑗

(22)

where: 𝑞𝑗 = displacement in modal coordinates in terms of harmonic indices 𝑓𝑗 = force in modal coordinates in terms of harmonic indices 𝜉𝑗 = fraction of critical damping for mode j 𝜔𝑗 = natural frequency of mode j

Since j represents any mode, Eqn. (22) represents nmod uncoupled equations in the nmod unknowns 𝑞𝑗 . The advantage of this uncoupled system with the mode-superposition method is that all the computationally expensive

matrix algebra has already been done in the eigenvalue solver, so extensive time-dependent equations of motion (Eqn. (18)) can be analyzed inexpensively in modal coordinates along with Eqn. (20). The modal coordinates 𝑞𝑗 are converted back into the system response displacements {𝑢} to the applied loading using Eqn. (20). In other words, the individual modal responses 𝑞𝑗 are superimposed to obtain the actual system response (i.e., mode-superposition). This set of uncoupled cyclic sector equations is solved while enforcing the compatibility boundary conditions between the sectors. In a harmonic forced vibration, the force term can be expressed as,

𝑓𝑗 = 𝑓𝑗𝑐𝑒𝑖Ω𝑡

(23)

10

J.B. Min, T.S.R. Reddy†, M.A. Bakhle

and

R.M. Coroneos, G.L. Stefko A.J. Provenza, K.P. Duffy†

NASA Glenn Research Center, Cleveland, Ohio †University of Toledo, Toledo, Ohio

Accurate prediction of the blade vibration stress is required to determine overall durability of fan blade design under Boundary Layer Ingestion (BLI) distorted flow environments. Traditional single blade modeling technique is incapable of representing accurate modeling for the entire rotor blade system subject to complex dynamic loading behaviors and vibrations in distorted flow conditions. A particular objective of our work was to develop a high-fidelity full-rotor aeromechanics analysis capability for a system subjected to a distorted inlet flow by applying cyclic symmetry finite element modeling methodology. This reduction modeling method allows computationally very efficient analysis using a small periodic section of the full rotor blade system. Experimental testing by the use of the 8-foot by 6-foot Supersonic Wind Tunnel Test facility at NASA Glenn Research Center was also carried out for the system designated as the Boundary Layer Ingesting Inlet/Distortion-Tolerant Fan (BLI2DTF) technology development. The results obtained from the present numerical modeling technique were evaluated with those of the wind tunnel experimental test, toward establishing a computationally efficient aeromechanics analysis modeling tool facilitating for analyses of the full rotor blade systems subjected to a distorted inlet flow conditions. Fairly good correlations were achieved hence our computational modeling techniques were fully demonstrated. The analysis result showed that the safety margin requirement set in the BLI2DTF fan blade design provided a sufficient margin with respect to the operating speed range.

I. Introduction

THE concept of relocating the engines from a podded position on the wing or fuselage and embedding them within the airframe [1] (Figure 1) can improve vehicle fuel burn, weight, noise, and emissions (through reduced fuel burn and engine cycle improvements). This Hybrid-Wing-Body (HWB) architecture is particularly suited to embedded engines, as balance requirements already place them near the aft of the airframe.[2] The embedded engines on a HWB can be integrated into the fuselage and ducted to exhaust through the rear of the vehicle. A main feature of this proposed boundary layer ingested inlet distortion-tolerant fan aircraft is the embedded aft-mounted engines which ingest part of the fuselage boundary layer air flow to reduce fuel burn. This embedded aft engines may also provide reduced susceptibility to bird strike since engines are partly hidden head-on especially at take-off, and possibly enabling to reduce the engine noise to the ground from positioning the engines above the airframe. The National Aeronautics and Space Administration (NASA) and United Technologies Research Center (UTRC) with contributions from Virginia Polytechnic and State University (Virginia Tech) have developed a BLI propulsor to significantly reduce the fuel consumption. Accordingly, a high-efficiency Boundary Layer Ingestion (BLI) flow Distortion-Tolerant Fan (DTF) blade rotor system has been studied as part of the Boundary Layer Ingesting Inlet/Distortion-Tolerant Fan (BLI2DTF) Project.

1

Inlet airflow distortion ingestion into a fan gives rise to additional aero-mechanical vibration and a decrease in stability margin. An embedded BLI propulsion system inherently involves a persistent and severe inlet distortion reaching the fan at most operating conditions, resulting in high dynamic stresses due to aeroelastic forced responses and the possibility of flutter. If not adequately addressed by the propulsor design, these aeromechanics phenomena could cause fan structural failure due to high cycle fatigue (HCF) or unstable self-excited failure, respectively. Numerical simulation modelling for this dynamically interconnected problem is a challenging problem.

Figure 1: Conceptual aircraft design and boundary layer ingestion inlet distortion-tolerant rotor-fan system (inset at the bottom: BLI2DTF fan in the 8’x6’ Supersonic Wind Tunnel test at NASA GRC)

Aero-mechanical vibration stresses in fan blades due to inlet airflow distortion is of particular interest in this study since the increase of stress in a blade can lead to its premature failure. Blade failure is usually occurring within short time by unstable self-excited (flutter) overload and high cycle fatigue (HCF) dynamic load. Turbomachinery blade vibration stress problem has been reported as the single cause to termination in new engine development effort. “While over 90 percent of the potential HCF problems are uncovered during development testing of a new engine, the remaining few account for nearly 30 percent of the total development cost and are responsible for over 25 percent of all engine distress events.”[3]

One of the main research objectives of the BLI2DTF Project was to develop and apply aeromechanics analysis tools to predict the effects of unique BLI distortion on “designed to be distortion-tolerant” fan system.[4,5] Traditional single blade modeling techniques are incapable of accurately modeling the entire rotor blade system due to complex dynamic loading behaviors and vibrations in distorted flow conditions. Thus, the present aeromechanics analysis on BLI2DTF fan rotor system had four specific objectives: (1) establish a model based approach that enables a highfidelity full-rotor aeromechanics analysis capability for performance, safety, and computational efficiency in design analysis of a fan rotor system subjected to a distorted inlet flow, (2) complete forced response analyses to determine dynamic stresses in DTF fan blade due to continual operation in a distorted flow resulting from BLI using cyclic symmetry modeling techniques, (3) assess the concerns of fatigue problem using Goodman Diagram approach by examining predicted stress values regarding the HCF behavior due to blade static and vibration dynamic stresses, and (4) provide guidance on the strain measurement by determining optimal strain gauge locations on the blades towards a cost-effective experimental test measurement from overall and/or directional strain characteristics predicted due to a loading condition aforementioned. This latest high-fidelity modeling capability enables the analysis of arbitrary inlet distortion patterns that can be specified through variations of flow properties, together with prescribed characteristics of DTF fan blade vibrations. The present analysis approach will ensure operability of DTF fan blades in the embedded propulsion system of HWB aircraft.

II. Cyclic symmetry finite element modeling technique

Understanding complex mechanical and aerodynamic interaction of turbomachine blading in this distorted flow environment has become topmost to the success of BLI2DTF fan design and durability. As mentioned above, typical single blade analysis models widely employed in the turbomachinery research community are incapable of representing the design for both understanding the harmonic modal blade dynamic behaviors and accurately predicting

2

corresponding forced response information of the blades subjected to complex distorted flow dynamic loading and vibrations. Nevertheless, analyzing this type of the problems by modeling a full 360 degree rotor blades structural system is not computationally efficient and practical, especially in large size systems. With addressing this issue in the present DTF fan analysis, reduction modeling techniques were essential and a novel reduction modeling capability was developed for fan distortion analysis to model non-axisymmetric large turbomachinery BLI2DTF rotor blades system. To achieve this capability, cyclic symmetry analysis modeling formulations were incorporated with the Ansys finite element (FE) computer program[6] to analyze the BLI2DTF forced responses.

Cyclic symmetry is a type of symmetry in which identical sectors are attached in a circular pattern about a rotation axis. It is assumed that in each individual sector’s local coordinate system the stiffness and mass matrices for all of the sectors are identical. The reduction of degrees-of-freedom is accomplished at the sector level. The compatibility of sector interfaces generates a reduced set of equations of the complete structure. A key element of cyclic symmetry method is to determine the cyclic symmetry transformation. That is, the solution of the complete structure can be reduced to the solution of the cyclic symmetric substructures. Once the solution for the substructure is determined, the solution of the complete structure is given by the cyclic symmetry transformation. The general procedure for the implementation of this method is given in reference 7.

The sector geometry of a BLI2DTF fan rotor system was meshed[8] with higher order 20-node hexahedral elements for the airfoil region. Higher order 10-node tetrahedral elements for the airfoil platform and the disk part were more densely meshed for better prediction of stress results. This sector model represents a cyclic symmetric sector of a whole DTF fan rotor system, shown in Figure 2.

(a)

(b)

(c)

Figure 2: a) FE mesh of a cyclically symmetric “base sector” model of BLI2DTF (close-up), b) description of cyclic symmetry interface boundaries, c) fully expanded mesh by cyclic symmetry formulation

Eighteen blades were placed symmetrically around the disk, so a 20 degree sector was used to produce a symmetric sector. The airfoil itself was meshed with 4 hexahedral elements through the thickness in order to properly capture the bending stiffness. The resulting mesh of a sector geometry having 557,554 elements and 992,548 nodes is shown in Figure 2 which requires to solve 5,912,580 equations for a particular engine operating condition. It was anticipated to be able to analyze this extremely large model by solving the equations only for a smaller sized sector model without losing the computational accuracy rather than analyzing the entire 360 degree rotor blade model requiring unrealistic computing time and hardware capacity.

The cyclic sector model has two edges that align along the surfaces of cyclic symmetry (Fig. 2, (a)) in angles. The edge having the algebraically lower angle, θ, in the cylindrical coordinate system is called the low edge and the one having the higher angle, θ + α, is called the high edge. The angle, α, between the two successive surfaces of cyclic symmetry is called the sector angle. The architecture of the cyclic symmetry solution process depends upon how the compatibility and equilibrium conditions of the cyclic sector are enforced in the matrix-solution process using the duplicate sector method. [7]

The principal features of this cyclic symmetry modeling method are briefly described below. Along with the notations illustrated in Figure 2(b), all loading, boundary conditions, and coupling and constraint equations being

3

present on the basic sector (A) are applied to the duplicate sector (B). Constraint relationships (Eqn. (1)) can be defined to relate the lower (θ = 0) and higher (θ = α) angle edges of the basic sector to allow calculations related to a given number of harmonic indices (k) defined below. Essentially, the basic sector is duplicated in the analysis steps to satisfy the required constraint relationships and to obtain nodal displacements.

{𝑢𝐴′ } = [ 𝑐𝑜𝑠𝑘𝛼 𝑠𝑖𝑛𝑘𝛼 ] {𝑢𝐴}

(1)

𝑢𝐵′ −𝑠𝑖𝑛𝑘𝛼 𝑐𝑜𝑠𝑘𝛼 𝑢𝐵

where:

𝑢𝐴, 𝑢𝐵 = calculated displacement on lower angle side of basic and duplicated sectors (A and B, respectively) 𝑢𝐴′ , 𝑢𝐵′ = displacements on higher angle side of basic and duplicated sectors (A and B, respectively) determined from constraint relationships

𝑁 𝑖𝑓 𝑁 𝑖𝑠 𝑒𝑣𝑒𝑛 k = harmonic index = 0, 1, 2…. {𝑁2−1

2 𝑖𝑓 𝑁 𝑖𝑠 𝑜𝑑𝑑 𝛼 = 2𝜋⁄𝑁 = sector angle N = total number of sectors in 360 degree (N = 18)

The cyclic symmetric solution sequences consist of three basic steps. The first step transforms applied loads to cyclic symmetric components using finite Fourier theory and enforces cyclic symmetry constraint equations (Eqn. (1)) for each harmonic index (i.e. nodal diameter) (k = 0, 1, . . ., N/2). At the solution stage of a cyclic symmetry analysis, the cyclic symmetry compatibility conditions are enforced for each harmonic index solution via coupling and/or constraint equations (CEs) connecting the nodes on the low- and high-edge components on the basic and duplicate sectors. As stated earlier, the cyclic symmetry solution process depends upon how the compatibility and equilibrium conditions of the cyclic sector are enforced in the matrix-solution process. The solution stage performs each harmonic index solution as a separate load step where the real (basic sector) and imaginary (duplicate sector) parts of the solution are calculated, as expressed below in Eqn. (2) – Eqn. (5). The Fourier transformation from physical components, u, to the different harmonic index components, 𝑢̅, is given by,

For harmonic index, k = 0 (symmetric mode), 𝑢̅𝑘=0 = 𝑁1 ∑𝑁𝑛=1 𝑢𝑛 (2)

For harmonic index, 0 < k < N/2 (mix mode),

basic sector: (𝑢̅𝑘)𝐴 = 𝑁2 ∑𝑁𝑛=1 𝑢𝑛 cos(𝑛 − 1) 𝑘𝛼 (3)

duplicate sector: (𝑢̅𝑘)𝐵 = 𝑁2 ∑𝑁𝑛=1 𝑢𝑛 sin(𝑛 − 1) 𝑘𝛼 (4)

For N even only, k = N/2 (anti-symmetric mode): 𝑢̅𝑘=𝑁/2 = 𝑁1 ∑𝑁𝑛=1(−1)(𝑛−1)𝑢𝑛 (5)

where: u = displacements (or any physical component, such as, forces, pressure loads, temperatures, and inertial loads 𝑢̅ = cyclic displacements (or any cyclic symmetric component) n = sector number, varies from 1 to N

In the second step, any applied load on the full 360 degree fan model is treated through a Fourier transformation process above and applied on to the cyclic basic sector. For each value of harmonic index, k, the procedure solves the corresponding equation. In the third step, the responses are then combined to get the complete response of the full fan structure via the Fourier expansion given by Eqn. (6). That is, using the cyclic symmetry solution of the basic and duplicate sectors connecting low and high edges of basic and duplicate sectors combines the solutions from the two sectors by performing computations on the selected load step to combine the results of the two sectors. Note that only a harmonic index zero solution (k=0) is valid for a static solution with cyclic loading.

4

The transformation to physical components, u, from the cyclic symmetry, 𝑢̅ , components is recovered by this equation,

𝑢𝑛= ̅𝑢𝑘=0 + ∑𝑁𝑘=1[ (𝑢̅𝑘)𝐴 cos(𝑛 − 1)𝑘𝛼 + (𝑢̅𝑘)𝐵sin(𝑛 − 1)𝑘𝛼] + (−1)(𝑛−1) ̅𝑢𝑘=𝑁/2

(6)

The last term (−1)(𝑛−1) ̅𝑢𝑘=𝑁/2 exists only for a case of N even. In cyclic symmetry modal analysis, the mode shape in each sector is obtained from the eigenvalue solution. The

cyclic symmetry modal solution occurs after the static cyclic symmetry solution, so that the modal analysis yields the mode shapes of the fan blade in the pre-stressed state. The displacement components (x, y, or z) at any node in sector n for harmonic index k, in the full fan structure are yielded by Eqn. (6).

If the mode shapes are normalized to the mass matrix in the modal analysis, the normalized displacement components in the full fan blade system are given by,

𝑢 √𝑁

𝑖𝑓 𝑘 = 0 𝑜𝑟 𝑘 = 𝑁⁄2

𝑢𝑛𝑜𝑟𝑚𝑎𝑙𝑖𝑧𝑒𝑑 = { 𝑢

𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

(7)

√𝑁⁄2

In a cyclic symmetry modal analysis, it is necessary to understand the concepts of harmonic indices and nodal diameters. Examining at a corresponding point in each sector, the mode shapes are sinusoidal in the circumferential direction, which indicates nodal lines on the bladed disk system called as nodal diameters (ND). Most mode shapes contain lines of zero out-of-plane displacement which cross the entire disk. This modal behavior is described in more detail in Section III.

In addition to cyclic symmetry constraint equations defined above, surface-to-surface contact pairs were used to define contact behaviors at the blade disk dovetail assembly as illustrated in Figure 3. In the contact between two bodies, the surface of one body was taken as a contact surface and the surface of the other body as a target surface. The “contact-target” pair concept has been widely used in finite element simulations. The contact surface was meshed using contact elements which are a 3D, 8-node, higher order quadrilateral element located on the surfaces of blade dovetail area. The target surface was meshed using target elements which are a 3D, 8-node, higher order quadrilateral to represent 3D target surface for the associated contact elements on the surfaces of disk dovetail socket area. For these surface-to-surface contact elements, the augmented Lagrangian penalty method,[6,9] incorporated with a contact "spring" stiffness called the contact stiffness and target stiffness was used to establish a relationship between the contact-target surfaces for bonded contact in conjunction with cyclic symmetry constraint conditions. Contact problems are highly nonlinear and require significant computer resources. In the present analysis, it was assumed that both contacting bodies are deformable with similar stiffness so flexible-to-flexible contact type was employed. Generally the regions of contact are not known until an actual solution process is evolved with the loads, materials, boundary conditions, and other factors in a largely unpredictable manner.

contact elements on surface in contact with

disk

Figure 3: Surface-to-surface contact pairs between disk hub and blade in dovetail attachment

The augmented Lagrangian method used in this study is an iterative series of penalty updates to find the Lagrange multipliers (i.e., contact tractions) and was investigated with the magnitude of the contact stiffness coefficient from a common range of 0.01 - 1.0. It was observed from this investigation that the augmented Lagrangian penalty method was insensitive to the magnitude of the contact stiffness coefficient. Hence, this method was employed in the present

5

DTF fan forced response analysis, presumably well-suited for a wide range of frictional contact situations governed by a Coulomb law[9] at variable rotor speeds of DTF fan blade. Further boundary conditions are applied at the inner radius of the disk sector to restrain motion in the radial and circumferential directions to prevent the rigid body motions. Also, surface elements manipulated by thickness set to zero were overlaid onto an area face of the blade elements as a modeling technique to apply steady and unsteady aerodynamic pressure loads on the blade surfaces.

In materials used for this analysis, the blade part was made from a titanium alloy and the disk part was made from a steel.

III. Methods for forced response computation

A finite element forced response analysis method utilizing the cyclic symmetry modeling approach was employed to provide a comprehensive study of a variety of structural characteristics of a rotating DTF fan blade system. The first task (section A) was to examine the steady stresses including nonlinear rotational effects of spin and stress stiffening due to the centrifugal loading caused by rotational velocity and steady state aerodynamic blade pressure load. The second task (section B) was to perform a prestressed vibration modal analysis, allowing for the classification of every natural mode within the frequency range of interest. This computation also provided information about areas of high stress and how the relative magnitude of these stresses varies with frequency. Finally, section C, a cyclic symmetry finite element analysis model incorporated with Turbo aerodynamic analysis computer program[10] and Ansys structural dynamics computer program was developed to compute the forced response stress distribution in fan blade and as a result the blade structural safety margins for the targeted DTF blade design requirements are presented in section D. While all the equations in the following sections are represented in the cyclic symmetry coordinates, added symbols to the equations are avoided for a brevity in the descriptions.

A. CYCLIC SYMMETRY CENTRIFUGAL STATIC STRESS ANALYSIS MODEL SUBJECTED TO STEADY-STATE AERODYNAMIC PRESSURE AT FAN SPEEDS

The first task was to develop a cyclic symmetry finite element static analysis model incorporated with Turbo steady-state aerodynamic blade pressure load as a function of the fan blade rotation. To capture the true blade behavior in this rotational condition, geometric nonlinearity was to be accounted for providing the blade deflections potentially to cause significant changes in the blade structural geometry. In such events, while solving the equations of equilibrium by utilizing an iterative process would be employed to obtain the accurate prediction for the largely deformed configuration, excessively high computing time is required for the size of the present problem. Even if a small deflection analysis may not directly account for change in geometry, the effect for large deflection would be possible to be accounted for by an adjustment of the structural stiffness matrix in a small deflection solution approach. This approach allows adjusting the stiffness of a rotating blade to account for dynamic mass effects. Equilibrium of the system and centrifugal forces on the mass using small deflection logic can be expressed as,

[𝐾]{𝑢} = {𝛺2}[𝑀]{𝑟}

(8)

where:

[K] = global stiffness matrix {𝑢} = displacement of the mass from the rest position {𝑟} = rest position of the mass with respect to the axis of rotation 𝛺 = angular velocity of rotation [𝑀] = global mass matrix = ∑𝑒𝑒==𝑛1 𝑀𝑒 e = element 1, 2, , n [𝑀𝑒] = element mass matrix = ∫𝑉[𝑁]𝑇[𝑁]𝜌𝑑𝑉 [𝑁] = shape function matrix

ρ= element density

n= number of elements

However, to account for deflection effects, Eqn. (8) must be expanded to:

[𝐾]{𝑢} = {𝛺2}[𝑀]{𝑟 + 𝑢}

(9)

Rearranging terms,

([𝐾] − {𝛺2}[𝑀]){𝑢} = {𝛺2}[𝑀]{𝑟}

(10)

6

Defining:

[𝐾̅] = [𝐾] − {𝛺2}[𝑀]

(11)

[𝐹̅] = {𝛺2}[𝑀]{𝑟}

Eqn. (10) becomes simply,

[𝐾̅]{𝑢} = [𝐹̅]

(12)

[𝐾̅] is the stiffness needed in a small deflection solution to account for large deflection effect, and [𝐹̅] is the same as that derived from small deflection logic. Rotational velocities are combined with the element mass matrices to form a body force load vector. Thus, the large deflection effects are included in a small deflection solution. This decrease in the effective stiffness matrix is called spin-softening.[11] Additionally, with assuming that the latest stress state in the computation may be large enough to change the blade shape under high centrifugal rotation, inclusion of such stiffness effects may need to be examined for the potential changing geometry due to this deformation by large stress which may no longer be neglected. This effect called stress stiffening effect was also included in the computation. This stress stiffening is the stiffening (or weakening) of a structure due to its stress state. This stiffening effect normally needs to be considered for thin structures with bending stiffness very small compared to axial stiffness, such as thin DTF blades which couples the in-plane and out-of-plane displacements. Accordingly, this stress stiffness matrix was added to the stiffness matrix in Eqn. (12) in order to update the total stiffness for accurate predictions. Hence, the equilibrium equation for static analysis is given by,

([𝐾̅] + [𝑆̅]){𝑢} = [𝐹̅]

(13)

where:

{𝑢} = nodal displacement vector [𝐾̅] = global stiffness including spin softening stiffness due to centrifugal force [𝑆̅] = stress stiffening (or weakening) due to its stress state [𝑆̅] = ∫ [𝐺]𝑇[𝜏] [𝐺]𝑑𝑉 when necessary

𝑉

[𝐺] = matrix shape function derivatives

[𝜏]= matrix of current Cauchy (true) stresses in global coordinate [𝐹̅] = global load vector due to rotating body force and steady state aerodynamic pressure load

As for the mass term, the mass summary is calculated from accumulating each element contribution which does not reflect the boundary conditions, coupled degrees of freedom, constraint equations, and multipoint constraint in contact elements.

B. PRESTRESSED CYCLIC SYMMETRY VIBRATION MODAL ANALYSIS MODEL SUBJECTED TO CENTRIFUGAL STRESS WITH STEADY-STATE AERODYNAMIC PRESSURE AT FAN SPEEDS

Given that the cyclic symmetry harmonic forced responses are to be computed in terms of the unsteady state aerodynamic engine order (EO) excitation vibration loads, a perturbed prestressed cyclic symmetry modal analysis model was developed after first performing a cyclic symmetry static analysis for identifying frequencies and mode shapes of blades vibration and critical speed ranges for a given set of nodal diameters (as illustrated in Figure 4). This allows for the classification of every natural mode within the frequency range of interest.

This vibration modal analysis assumed that there is no damping and the structure has no time varying forces applied (i.e., free vibration) to determine the cyclic symmetry natural frequencies and mode shapes. The equation of motion for an undamped system can be expressed in matrix notation,

[𝑀]{𝑢̈ (𝑡)} + [𝐾̃]{𝑢(𝑡)} = {0}

(14)

where: [𝐾̃] = [𝐾̅] + [𝑆̅]

For a linear system, free vibrations are harmonic of the form,

{𝑢(𝑡)} = {𝑢} = {𝜙}𝑟𝑐𝑜𝑠𝜔𝑟𝑡

(15)

7

where: {𝜙}𝑟 = eigenvector representing the mode shape of the rth natural frequency 𝜔𝑟 = rth natural frequency (cycles/second, Hz); 𝜔𝑟 = 2𝑓𝑟𝜋 where 𝑓𝑟 = natural circular frequencies (radians/second)

0 ND

1 ND

2 ND

3 ND

4 ND

…

9 ND

Note: ND = Nodal Diameter; Mode refers to Mode Family; Magnitude of displacement is plotted

Figure 4: Illustration of nodal diameter (ND) patterns from BLI2DTF model shown with higher deformation in red color

As indicated in Fig. 4, the cyclic symmetry prestressed modal analysis solution provides analytical results for the interaction between the disk vibration and the blade vibration of DTF fan blade. The blade frequencies were solved corresponding to a particular disk mode (or nodal diameter, ND). Different disk nodal diameter information is to capture the blade frequencies and mode shapes for each nodal diameter.

The cyclic symmetry modal solution uses the same low- and high-edge components defined in the static cyclic symmetry analysis. To calculate the frequencies and mode shapes of a deformed structure, the linear perturbation solution method based on the Newton-Raphson procedure[12] after performing a static analysis is needed for a prestressed modal analysis of cyclic symmetry structures. The linear perturbation modal analysis procedure is designed to solve a problem from this preloaded condition. The linear perturbation modal analysis can be viewed as an iteration solution on top of a prior analysis. During the linear perturbation process, all of the effects from the static analysis are taken into account and are “snapshot solution” so that the perturbed loads can be generated linearly by using the "snapshot solution” matrices and material properties. With a restart procedure from the static analysis, the iterative process is used for a restart solution by allowing modifications of perturbed loads to generate new matrices of perturbation loads, mass density, and damping. The load history vector from the modal analysis is updated for the next iteration solution and subsequently the modal coordinates are updated by solving eigenvalue solution. During this phase of the analysis, the frequency dependent perturbation load vector is calculated and stored in the mode shape file and the element load information files for a subsequent mode-superposition modal-based forced response analysis. This procedure allows for saving significant time in computation and reusing the mode shape information from an earlier modal analysis to improve the accuracy of a mode-superposition harmonic forced response computation. Modal extracting computation process typically requires more time than element loads generation, perturbation residual load vector calculation, and enforced snapshot solution mode calculation.

Hence, the equation of motion for an undamped system can be expressed for the linear perturbation vibration modal analysis to compute the natural frequencies and mode shapes,

8

(−𝜔𝑟2[𝑀] + [𝐾̃𝑙]){𝜙}𝑟 = {0}

(16)

where: [𝐾̃𝑙] = stiffness update at the iteration of perturbation load step “l” for prestressed modal analysis where spin-

softening effects and stress stiffening effects included

This equality is satisfied if either {𝜙}𝑟 = {0} or if the determinant of (−𝜔𝑟2[𝑀] + [𝐾̃𝑙]) is zero. The first option is the trivial one and, therefore, is not of interest. Thus, the second one gives the solution,

|[𝐾̃𝑙] − 𝜔𝑟2[𝑀]| = 0

(17)

This is an eigenvalue problem which may be solved for up to ndof values of 𝜔𝑟2 and ndof eigenvectors {𝜙}𝑟 which satisfy Eqn. (16) where ndof is the number of degrees-of-freedom (DOFs) of the model.

While several eigenvalue and eigenvector extraction procedures are available for solving Eqn. (17), the present study employed the Block Lanczos method [13] to extract the requested eigenvalues using the sparse direct solving technique that performs forward and backward substitution using the factors. During this solution process, the frequency dependent perturbation load is calculated and stored as nodal forces for a subsequent mode-superposition modal-based harmonic forced response analysis.

Repeat analysis was performed to calculate blade vibration frequencies at different rotational speeds. Also, as will be described in detail in Section IV, an additional objective of this prestressed cyclic symmetry modal analysis was to guide the optimal strain gauge placement on the fan blades by providing the information regarding how the magnitude of modal stresses varies with frequencies so that the optimal locations of the strain gauges are determined on a blade for the effective strain gauge measurement in the experimental wind tunnel test of fan blades.

C. CYCLIC SYMMETRY HARMONIC FORCED RESPONSE ANALYSIS MODEL SUBJECTED TO A COMPLETE CYCLE OF ENGINE ORDER FORCING VIBRATION BY UNSTEADY-STATE AERODYNAMIC PRESSURE LOAD AT FAN SPEEDS

Next, a cyclic symmetry harmonic forced response analysis model was developed utilizing the mode superposition method. The results of a cyclic symmetry modal analysis described in Section III.B was retrieved for a harmonic response analysis described in this section. This cyclic symmetry harmonic forced response analysis model allows for a complete solution of DTF fan blade vibration dynamic stresses throughout a complete cycle of engine order forcing vibration.

The finite element semi-discrete equation of motion for rotation including prestressed structural stiffening effects can be expressed in matrix form,

[𝑀]{𝑢̈ (𝑡)} + [𝐶]{𝑢̇ (𝑡)} + [𝐾̃𝑙]{𝑢(𝑡)} = {𝐹𝑎(𝑡)}

(18)

where: [𝐶] = structural damping matrix {𝑢̈ (𝑡)} = nodal acceleration vector {𝑢̇ (𝑡)} = nodal velocity vector {𝑢(𝑡)} = nodal displacement vector {𝐹𝑎(𝑡)} = unsteady aerodynamic EO excitation pressure load

t = time representing the current equilibrium iteration

l = load step in perturbation solutions

The cyclic symmetry harmonic forced response analysis solves the time-dependent equations of motion (Eqn. (18)) for DTF fan blade structures undergoing harmonic forced vibration. This analysis assumed that the entire DTF blade structure has frequency-dependent stiffness, damping, and mass effects and all loads and displacements vary sinusoidally at the same known frequency (although not necessarily in phase). Mathematically, Eqn. (18) represents a system of linear differential equations of second order. Two procedures can be considered for the solution of Eqn. (18): direct integration method and mode superposition method. Given that a traveling wave EO excitation occurs due to circumferential disturbances in the BLI flow field from upstream of the fan and a mode superposition method shows the numerical effectiveness over a direct integration method, engine order loading (circumferentially traveling wave excitation) is modeled with a mode-superposition analysis method.[6,14] Engine order excitation is periodic

9

forcing arising from circumferential disturbance in a set of blades on the bladed rotor. The mode-superposition method sums factored mode shapes using the natural frequencies and mode shapes obtained from a linear perturbed modal analysis to characterize the dynamic response of a harmonically excited structure. In addition, since a mode superposition method may need to consider only a few modes, the mode superposition method can be much effective than a direct integration which would be prohibitively expensive in computations if the order of the matrices is large.[14]

Eqn. (18) was solved by the mode-superposition method to the interest of specific EO exciting unsteady aerodynamic pressure load on the blades. {𝐹𝑎} is the time-varying load vector given by,

{𝐹𝑎} = {𝐹𝑛𝑑} + 𝑠𝑓. {𝐹𝑠}

(19)

where: {𝐹𝑛𝑑} = time varying unsteady aerodynamic nodal forces

sf = load vector scale factor {𝐹𝑠} = load vector read in from the linear perturbation modal analysis in Section III.B

Note that all loads from the linear perturbation modal analysis are scaled, including forces and pressures. All loads from the modal solution load vector during a given load step in the linear perturbation modal analysis is applied as engine order loads. Damping in Eqn. (18) should be specified; otherwise, the response will be infinite at the resonant frequencies. Additional information on the transformation of Eqn. (18) from cyclic symmetry coordinates to physical coordinates and Fourier transformation in cyclic symmetry forced response analysis are described in reference 6.

Hence, similarly given in reference 14, a general treatment of dynamic systems can be formulated by Lagrange’s equation from the scalar quantities of kinetic energy, potential energy, and work where the equations of motion of a system can be formulated in a number of independent coordinate systems. Such independent coordinates are called generalized coordinates and can be denoted by the character 𝑞𝑟 in the present study. Eqn. (18) can be transferred into modal space and hence, a set of modal coordinates 𝑞𝑟was defined such that,

{𝑢} = ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟

(20)

where: {𝜙𝑟} = rth mode shape nmod = number of modes to be used from the preceding modal analysis in Section III.B

Substituting Eq. (20) into Eqn. (18)

[𝑀] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟}𝑞𝑟̈ + [𝐶] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟̇ + [𝐾̃𝑙] ∑𝑛𝑟=𝑚1𝑜𝑑{𝜙𝑟} 𝑞𝑟 = {𝐹𝑎(𝑡)}

(21)

After few steps of the matrix transformation and substitution with applying its orthogonality and normality conditions, and only the r = j terms remain, Eqn. (21) can be expressed in the modal coordinates as,

𝑞𝑗̈ + 2𝜔𝑗𝜉𝑗𝑞𝑗̇ + 𝜔𝑗2𝑞𝑗 = 𝑓𝑗

(22)

where: 𝑞𝑗 = displacement in modal coordinates in terms of harmonic indices 𝑓𝑗 = force in modal coordinates in terms of harmonic indices 𝜉𝑗 = fraction of critical damping for mode j 𝜔𝑗 = natural frequency of mode j

Since j represents any mode, Eqn. (22) represents nmod uncoupled equations in the nmod unknowns 𝑞𝑗 . The advantage of this uncoupled system with the mode-superposition method is that all the computationally expensive

matrix algebra has already been done in the eigenvalue solver, so extensive time-dependent equations of motion (Eqn. (18)) can be analyzed inexpensively in modal coordinates along with Eqn. (20). The modal coordinates 𝑞𝑗 are converted back into the system response displacements {𝑢} to the applied loading using Eqn. (20). In other words, the individual modal responses 𝑞𝑗 are superimposed to obtain the actual system response (i.e., mode-superposition). This set of uncoupled cyclic sector equations is solved while enforcing the compatibility boundary conditions between the sectors. In a harmonic forced vibration, the force term can be expressed as,

𝑓𝑗 = 𝑓𝑗𝑐𝑒𝑖Ω𝑡

(23)

10