# Fluid/Structure Coupled Aeroelastic Computations for Transonic Flows in Turbomachinery

Hirofumi Doi, Juan Alonso

► To cite this version:

Hirofumi Doi, Juan Alonso. Fluid/Structure Coupled Aeroelastic Computations for Transonic Flows in Turbomachinery. Turbo Expo, International Gas Turbine Institute, Jun 2002, Amsterdam, Netherlands. pp.787-794, 10.1115/GT2002-30313. hal-01337402

HAL Id: hal-01337402

https://hal.science/hal-01337402

Submitted on 25 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FLUID/STRUCTURE COUPLED AEROELASTIC COMPUTATIONS FOR TRANSONIC FLOWS IN TURBOMACHINERY

**Hirofumi Doi**

Department of Aeronautics and Astronautics
Stanford University
Stanford, California 94305
Email: hdoi@stanfordalumni.org

**Juan J. Alonso**

Department of Aeronautics and Astronautics
Stanford University
Stanford, California 94305
Email: jjalonso@stanford.edu

## ABSTRACT

The present study demonstrates the capabilities of a fluid/structure coupled computational approach which consists of an unsteady three-dimensional Navier-Stokes flow solver, *TFLO*, and a finite element structural analysis package, *MSC/NASTRAN*. The parallelized flow solver relies on a multi-block cell-centered finite volume discretization and the dual time stepping time integration scheme with multigrid for convergence acceleration. High accuracy is pursued with respect to load transfer, deformation tracking and synchronization between the two disciplines. As a result, the program successfully predicts the aeroelastic responses of a high performance fan, *NASA Rotor 67*, over a range of operational conditions. The results show that the unsteady pressure generated at the shock may act to damp or excite the blade motion mainly depending on the inter-blade phase angle. It is concluded that the level of sophistication in the individually sophisticated disciplines together with an accurate coupling interface will allow for accurate prediction of flutter boundaries of turbomachinery components.

## INTRODUCTION

The unstable, self-excited or forced vibrations of rotor blades must be avoided in the design of high performance turbomachinery components because they may induce structural fail-

ures. In order to predict these instabilities, the presence of strong shocks in the flow needs to be accounted for, especially for transonic flows. It is, therefore, necessary to use the Euler equations or the Reynolds-averaged Navier-Stokes (RANS) equations to represent the unsteady flow fields. Unsteady aerodynamics around oscillating cascades has been studied using these nonlinear equations by many researchers. For these approaches, time-marching methods (Gerolymos, 1993; He, 1994; Bakhle, 1997; Ji, 1999) or time-linearized methods (Hall, 1993; Ning, 1998) are used to calculate aerodynamic work per cycle over a period of oscilation prescribing frequencies and mode shapes.

On the other hand, the coupling of a structural model and a fully nonlinear aerodynamic model requires a time-marching method and determines the frequency of the problem rather than specifying it as an input parameter. The most noticeable work on the coupled computations is done by Vahdati, Imregun and their colleagues (Vahdati, 1995; Chew, 1998). In their work, a mode superposition of the structure is incorporated into a finite element RANS solver. However, since the methodologies of each individual discipline have matured independently, each methodology has usually evolved to use a different type of grid generation, a different discretization method and a different time integration scheme so that high accuracy and efficiency individually can be achieved. In order to take advantage of the maturity of both, a reasonable alternative would be to construct an interface procedure between a flow solver and a structural solver for direct integrations of structural equations, in which the two solvers exchange the interface information given by updating the fluid and
structural variables alternatively.

The present study explores this possibility by integrating an unsteady RANS flow solver and a finite element structural solver for aeroelastic problems in turbomachinery, using advanced fluid/structure coupling techniques. The flow solver used here is an unsteady three-dimensional RANS solver called TFLO (Yao, 2000; Yao, 2001), originally developed to simulate unsteady flows due to blade row interactions in turbomachinery. The structural solver used here is one of the industrial standards, MSC/NASTRAN. Furthermore, the interface procedure is based on the approach proposed by Brown (Brown, 1997). Using this coupling approach, the aeroelastic responses of a compressor rotor are predicted, and the influence of the flow fields on the stabilities is observed.

## DESCRIPTION OF THE METHOD

In order to predict the dynamic response of a flexible structure in a fluid flow, the equations of motion of the structure and the fluid equations must be interacted. One difficulty in handling the fluid/structure coupling numerically comes from the fact that the structural equations are usually formulated with material (Lagrangian) coordinates while the fluid equations are typically written using spatial (Eulerian) coordinates. In such an approach, the procedure is advanced in time by solving the flow field and the structural deformation alternatively by independent flow and structural solvers which exchange information on the structural body surface as illustrated in Figure 1. The flow solver provides the aerodynamic loads to the structural solver in order for the structural solver to calculate the displacement field of the structure. In return, the structural solver provides the surface deflections to the flow solver which changes the flow fields through the boundary conditions on the structural body surface.

### Flow Solver: TFLO

TFLO is an unsteady, three-dimensional, turbomachinery flow solver in which the unsteady RANS equations are solved using a cell-centered discretization on multiblock meshes. The solution procedure is based on explicit Runge-Kutta methods with several convergence acceleration techniques such as multigrid, residual smoothing and local time-stepping. A parallelization strategy is adapted to TFLO based on the multiblock decomposition (Reuther, 1997) using the Message Passage Interface (MPI) standard (www-unix.mcs.anl.gov/mpi/index.html). For turbulence models, the *k*-ω model (Yao, 2000) is found to be the most suitable choice among the implemented models for complicated viscous transonic flows in turbomachinery with timely changing grids in this particular study.

For an aeroelastic analysis in a time-marching fashion, it is necessary to deform the fluid mesh at each physical time step so that it continuously conforms to the instantaneous shape of

<figure><img src="image_2.png" /><figcaption>Figure 1. FLUID/STRUCTURE COUPLING</figcaption></figure>

the aeroelastically deforming body under consideration. In the present study, a high quality mesh of the whole computational domain is divided into an appropriate multiblock system for parallelization purposes. Because the flow solver assumes a point-to-point match between adjacent blocks in the mesh, each block may be independently perturbed, provided that modified surfaces are treated continuously across block boundaries. For edge, face and block perturbation in these procedures, the arc-length attenuation algorithm (Reuther, 1996) is used.

### Fluid-Structure Interface

While the flow solver employs relatively finer multiblock structured meshes, the structural solver employs triangular elements. This mismatch of mesh points on the interface gives rise to the necessity of a numerical approximation to transfer the pressure distribution on the fluid surface mesh to a corresponding nodal load distribution on the structure and also for the transfer of the displacements of the structural nodes to corresponding perturbations of the fluid mesh.

In addition, the fluid system usually requires a finer temporal resolution than the structural vibration for most aeroelastic problems. In computations, while the flow solver employs the Runge-Kutta time integration which has a stability limit, the structural solver uses the Newmark time integration (Sitton, 1997) which is unconditionally stable. This difference in the stability limits need to be taken care of by a synchronization approach between the two solvers alternatively advancing in time.

**Conservation of Loads and Energy** One of the principles that is needed to support the maintenance of coupling accuracy is the conservation of loads and energy. In general while
the fluid system addresses the pressure field on the cell surfaces on the interface, the structural system is solved based on a set of concentrated forces at the nodes on the interface. A distributed pressure load, therefore, must be first transfered into equivalent nodal forces. Such a transformation must satisfy two requirements. One is that the nodal forces must yield the same net forces as the original distributed pressure loads do. Thus,

$$ \sum_m \mathbf{f}^{(m)} = \int_{\partial\Omega} p d\mathbf{S}, \qquad (1) $$

where $\mathbf{f}^{(m)}$ is the nodal force vector at the node $m$ in the structure, and $p$ is the pressure distribution on the surface $\partial\Omega$ whose projected area vector is $\mathbf{S}$. The second requirement states the maintenance of a proper energy balance. Equating the virtual work performed by $\mathbf{f}^{(m)}$ working on a virtual nodal displacement $\delta\mathbf{q}^{(m)}$ and that by $p$ moving through the equivalent distributed virtual displacement $\delta\boldsymbol{\mu}$, the second requirement is given as follows,

$$ \sum_m \mathbf{f}^{(m)} \delta\mathbf{q}^{(m)} = \int_{\partial\Omega} p \delta\boldsymbol{\mu} d\mathbf{S}. \qquad (2) $$

In this study, the approach by Brown (Brown, 1997) is chosen to ensure that the transfer of the pressure fields to the nodal forces is both consistent and conservative. Brown's approach can also be used in the extrapolation of the nodal displacement in the structural system to the mesh deformation on the surface of the fluid system as explained in the following section.

**Deformation Tracking System** Consider one of a set of finite elements describing a structural model and a fluid mesh point on the blade surface shown in Figure 2. Let $\mathbf{x}$ denote the closest point on the closest element from a fluid mesh point $\mathbf{X}$, it can be assumed that the vector connecting both points remain perpendicular to the element after deformation. Then, the displacements $\boldsymbol{\mu}(\mathbf{X})$ and rotations $\boldsymbol{\mu}_{\theta}(\mathbf{X})$ at $\mathbf{X}$ can be expressed as follows,

$$ \boldsymbol{\mu}(\mathbf{X}) = \mathbf{u}(\mathbf{x}) - (\mathbf{X} - \mathbf{x}) \times \mathbf{u}_{\theta}(\mathbf{x}), \quad \boldsymbol{\mu}_{\theta}(\mathbf{X}) = \mathbf{u}_{\theta}(\mathbf{x}), \qquad (3) $$

where $\mathbf{u}(\mathbf{x})$ and $\mathbf{u}_{\theta}(\mathbf{x})$ are the displacement and rotation at $\mathbf{x}$. In finite element analysis, the displacement at any point within the domain of an element in the model can be determined by the assumed interpolation functions in terms of the nodal displacements on the element. However, since MSC/NASTRAN is used as a structural solver and the interpolation function used in it is

<figure><img src="image_3.png" /><figcaption>Figure 2. DEFORMATION TRACKING SYSTEM</figcaption></figure>

unknown, all of the finite elements were assumed to use the simple standard iso-parametric interpolation function based on the area coordinate as the weighting (Reuther, 1999). When $\mathbf{X}$ is associated with $\mathbf{x}$ on the $m$-th element whose displacements are denoted as $\mathbf{q}^{(m)}$, this interpolation function $\eta^{(m)}$ and $\eta_{\theta}^{(m)}$ can be expressed as,

$$ \mathbf{u}(\mathbf{x}) = [\eta^{(m)}(\mathbf{x})] \cdot \mathbf{q}^{(m)}, \quad \mathbf{u}_{\theta}(\mathbf{x}) = [\eta_{\theta}^{(m)}(\mathbf{x})] \cdot \mathbf{q}^{(m)}. \qquad (4) $$

The displacements at any point on the fluid mesh surface can thus be written as,

$$ \boldsymbol{\mu}(\mathbf{X}) = [\eta^{(m)}(\mathbf{x})] \cdot \mathbf{q}^{(m)} - [\mathbf{X} - \mathbf{x}] \cdot [\eta_{\theta}^{(m)}(\mathbf{x}_n)] \cdot \mathbf{q}^{(m)}, \qquad (5) $$

where, the matrix $[\mathbf{X} - \mathbf{x}]$ is in cross product form. Introducing the displacement extrapolation functions $N(\mathbf{X})$ based on the global nodal displacements $\mathbf{q}$ of the structural model, Equation 5 can be rewritten using $\mathbf{q}$ because $\mathbf{q}^{(m)}$ is a part of $\mathbf{q}$.

$$ \boldsymbol{\mu}(\mathbf{X}) = [[\eta(\mathbf{x})] - [\mathbf{X} - \mathbf{x}] \cdot [\eta_{\theta}(\mathbf{x})]] \cdot \mathbf{q} = [N(\mathbf{X})] \cdot \mathbf{q} \qquad (6) $$

**Load Transfer System** For a load transfer algorithm, it is required to satisfy the conservation of load and energy expressed in Equations 1 and 2, respectively. In Equation 2, $\delta\boldsymbol{\mu}(\mathbf{X})$ can be related with $\delta\mathbf{q}$ introducing a set of assumed displacement interpolation functions $N(\mathbf{X})$ given in Equation 6. Substituting Equation 6 into Equation 2 yields,

$$ \mathbf{f} \cdot \delta\mathbf{q} = \sum_m \mathbf{f}^{(m)} \delta\mathbf{q}^{(m)} = \int_{\partial\Omega} p \cdot [N(\mathbf{X})] \delta\mathbf{q} d\mathbf{S}, \qquad (7) $$
where **f** is the global nodal force vector for the structural model. Consider the complementary projected surface area vector **S**<sub>i</sub> around the fluid mesh point **X**<sub>i</sub>, **S**<sub>i</sub> is surrounded by the corner points which are given by simply averaging the corner points of four fluid mesh faces around **X**<sub>i</sub> as shown in Figure 3. Thus the right hand side of the Equation 7 can be discretized as follows,

$$ \mathbf{f} \cdot \delta\mathbf{q} = \left(\sum_i \mathbf{f}_i\right) \cdot \delta\mathbf{q} = \sum_i p_i \mathbf{S}_i \cdot [N(\mathbf{X}_i)] \delta\mathbf{q}, \quad (8) $$

where **f**<sub>i</sub> is the nodal force vector given by the contribution of the pressure at **X**<sub>i</sub>, which components are shown in Figure 3. **f**<sup>(m)</sup> is now obtained by simply summing the all **f**<sub>i</sub> whose fluid mesh point is associated with the m-th element. Equation 8 allows the fluid pressure within each face to be taken out of the summation,

$$ \mathbf{f}_i = \{p_i \mathbf{S}_i\} [N(\mathbf{X}_i)], \quad (9) $$

Since N(**X**<sub>i</sub>) consists of the area coordinate system normalized by the entire area of the element associated with the point i, it is obvious that the norm of the function is equal to one. Therefore, the conservation of loads is validated by collecting **f**<sub>i</sub>.

$$ \sum_m \mathbf{f}^{(m)} = \sum_i \mathbf{f}_i = \sum_i \{p_i \mathbf{S}_i\} [N(\mathbf{X}_i)] = \sum_i p_i \mathbf{S}_i = \int_{\partial\Omega} p d\mathbf{S}. \quad (10) $$

Thus the algorithm satisfies the conservation of loads and energy.

**Synchronization** If a same time step in both fluid and structure system is used to solve the aeroelastic equations, the time step must be governed only by the critical time step of the explicit fluid solver. It is also noted that since the fluid system must be advanced in the physical time scale, the fast convergence methods such as the local time stepping, or multigrid can not be available. Thus, using the same time step is not efficient in advancing the aeroelastic system. This difference in allowable time step size is, however, overcome by the dual time stepping scheme (Jameson:1991) in which the fluid system is updated implicitly in physical time, and the pseudo time step is introduced to integrate the fluid equations within the physical time period. In this way, a same large physical time step is used for both solvers.

In aeroelastic calculations, a steady flow is first computed around a structure in equilibrium or in a given stationary condition with a set of initial displacements. Next the structure is allowed to respond to the aerodynamic forces determined by the

<figure><img src="image_4.png" /><figcaption>Figure 3. LOAD TRANSFER SYSTEM</figcaption></figure>

steady flow field. This should be followed by a simple and popular staggered procedure as described below.

(1) Transfer the motion of the boundary of the structure to the fluid system using the Deformation Tracking System and update the fluid mesh accordingly. At this time, the updated cell volumes, projected areas and surface normals are calculated.

(2) Advance the fluid system using the dual time stepping scheme to a steady-state and obtain new pressure fields using the solutions at the previous time levels stored in memory for the third order accurate backward time discretization.

(3) Convert the new pressure fields into a set of nodal forces for the structural model using the Load Transfer System.

(4) Advance the structural system using the Newmark scheme under the given set of nodal forces.

The procedure repeats from step (1) through (4) until the desired time-marching solutions are obtained.

## RESULTS

The goal of this study is to predict the aeroelastic response of a compressor, to examine the influence of the inter-blade phase angle and to observe how the flow field affects the instability of the blade vibrations. NASA Rotor 67 has been recognized as a popular test case for three-dimensional viscous flow prediction procedures (Chima, 1991; Jennions, 1993; Arnone, 1994) because of its detailed experimental data using a laser anemometer (Strazisar, 1989). Although the structural properties for this configuration are not available, Rotor 67 has been also used to validate aeroelastic applications by a few researchers (Chuang, 1998; He, 1994; Vahdati, 1995). Thus, Rotor 67 was chosen for the current investigation to demonstrate the capability of the
previously described nonlinear aeroelastic prediction procedure. NASA Rotor 67 consists of 22 blades. The blade aspect ratio is 1.56. At the design point, the rotational speed of the rotor is 16,043 RPM, with a tip inlet relative Mach number of 1.38, a total pressure ratio of 1.63 and a mass flow of 33.25 Kg/s.

The grid used for the present study is shown in Figure 4 with the blade-to-blade view. The grid for the single passage consists of 137 mesh points in the stream-wise direction, 65 points in the blade-to-blade direction, and 81 points in the span-wise direction. In the stream-wise direction, the grid extends one axial chord upstream and downstream from the blade row at the hub, and has 33, 73 and 33 axial points along the upstream periodic boundary, in the blade passage and along the downstream periodic boundary, respectively. A total of 17 points are used in the span-wise direction to describe the tip-clearance.

A NASTRAN finite element model for the blade consists of 11×11 nodes forming 200 varying-thickness three-node triangular linear shell elements, TRIA3. The thickness is distributed to match the actual thickness of the blade. For a boundary condition, both translational and rotational displacements are not allowed at the root of the blade. The direct transient response solution sequences, SOL129, (Herting, 1997) are used to solve the structural response of the blade with time-dependent nonlinear aerodynamic loads with initial conditions provided by the displacements and velocities of the last transient solution sequence. In modeling NASA Rotor 67, an imaginary material whose Young's modules, Poisson's ratio, and density are chosen to be 1.422e+11 Pa, 0.3 and 4539.5 kg/m³, respectively, such that the first natural frequency becomes around 400 Hz.

For the design rotation speed of the Rotor 67, the reduced frequency based on the tip section becomes 0.27 according to the first mode frequency calculated with the structural model by MSC/NASTRAN. Finally, for multi-passage simulations, the cascade is assumed to be well tuned, which means that all the blades have an identical structural property.

## Steady Flow Calculations

Prior to performing the aeroelastic unsteady calculations, validation of the steady flow solution needs to be performed such that the proper initial conditions for unsteady calculations can be provided. Steady flow computations are performed with the measured total pressure, total temperature and absolute flow angle profiles specified at the inlet and the measured static pressure profile specified at the exit.

To show the ability of the flow solver to capture the correct performance, the numerically calculated performance map of the rotor is shown in Figure 5 along with the experimentally measured values. As shown in Figure 5, the overall predicted performance of the compressor is in good agreement with the experiment except for the obvious underpredicted absolute total pressure ratio. This underprediction is likely due to an over-

prediction of shock and viscous flow losses and flow deviation. However, the steady-state solution including the structural deformation mentioned in the next section results in a higher mass flow and higher pressure ratio. Although the structural model is not based on the real structure of Rotor 67, this suggests that the underestimated pressure ratio and mass flow may be caused by the structural deformation not included in the present calculations.

Because of the availability of the measured data, detailed comparisons between the calculations and the experiment are made at the near peak efficiency point. In Figure 6, relative Mach number contours at three different spanwise locations, i.e. 10, 30, and 70 percent span measured from the casing, are shown for the near peak efficiency condition. At 10 percent of span, a bow shock at the leading edge and a normal in-passage shock near the trailing edge on the suction surface can be observed. At 30 percent of span, a lower Mach number moves the in-passage shock slightly forward but the shock structure is qualitatively similar. At 70 percent of span, the inlet relative Mach number is subsonic and a supersonic bubble appear near the leading edge on the suction surface.

## Aeroelastic Calculations

Next, the time dependent aeroelastic results are considered. The behavior of the blade oscillation can be determined by computing the aeroelastic response to a set of forced initial displacements. As previously shown, the steady-state solutions can be obtained for the rigid unloaded blade geometry. Prior to the aeroelastic response calculations, the blade is allowed to move and dynamically relax to a steady-state loaded condition. Then, the blade is deformed an additional amount predominantly corresponding to the first mode. Once the new steady solution for the blade with the initial deformations obtained, the blade is released

<figure><img src="image_5.png" /><figcaption>Figure 4. Blade-to-Blade View of Computational Grid for NASA Rotor 67</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Comparison of Rotor Performance at Design Speed</figcaption></figure>

and allowed to move according to the aeroelastic interaction.

A constant inter-blade phase angle $σ$ is determined by the initial displacement and velocity of the displacement for multi-passage computations. In addition, the inter-blade phase angle is forced by the periodic boundary condition. For example, consider the case for $σ$=180 degree. The computation needs to be performed with two passages where one blade is given an initial displacement and zero initial velocity of the displacement and the other is given the initial displacement in the opposite direction with the same amplitude and zero initial velocity of the displacement as well. In this study, cases for $σ$=0 (one passage) and $σ$=180 (two passages) degrees are examined in order to save computational costs.

Five cycles of oscillation, which is about 0.02 seconds, were found to be sufficient to distinguish whether the oscillation decays or grows. A total of 36 time steps per vibration period of the first oscillation mode was chosen to ensure the accuracy of the coupled system with the dual time stepping scheme and the Newmark scheme (Alonso, 1997). Between 40 to 60 W-cycle multigrid cycles per unit physical time step were necessary to provide a reasonable time-accurate pressure distribution to the structure.

The flutter stability is assessed from the time history of the

<figure><img src="image_7.png" /><figcaption>Figure 6. Experimental and Numerical Relative Mach Number Contour for Near Peak Efficiency</figcaption></figure>

displacements. Figure 7 shows the variation of the circumferential direction displacements of the blade with time for $σ$ = 0 and 180 degrees at the near peak efficiency point. The displacements are defined as the deflections from the original rigid shape of the structure. The response for $σ$ = 0 degree shows almost a typical response with constant amplitude of oscillation. On the other hand, the response for $σ$ = 180 degree shows decreasing amplitude with each oscillation cycle, indicating a decay in the oscillation and flutter free condition. Thus the aerodynamic damping for $σ$ = 180 degree is much greater than that for $σ$ = 0 degree.

To examine how the unsteady pressure has influence on the stability of the blade motion, the work per cycle is a parameter often used to represent the contributions of unsteady pressure fields to the stability. The work per cycle is the work done by the fluid on a given blade over one cycle of its motion. Negative work implies positive damping with which aerodynamic forces act in the opposite direction of the velocity of displacements. Figure 8 shows the calculated distribution of the aerodynamic work per cycle over the last oscillation cycle along with the distribution of the mean value of the unsteady pressure on the suction and pressure side, respectively. The results are presented for $σ$ = 0
<figure><img src="image_8.png" /><figcaption>Figure 7. DEFLECTION AT MID-CHORD OF THE TIP SECTION FOR NEAR PEAK EFFICIENCY, $\sigma = 0, 180^\circ$ deg</figcaption></figure>

degrees. As shown, significant work per cycle magnitudes exist in the region where the shocks sit on the blade. On the suction side, the shock appears along the entire span as a continuous line. However, a negative peak in the work per cycle distribution appears only at the location where the lower part of the shock sits on the suction side. Between the casing and 40 percent of the span, the passage shock is very strong and indicates that an adverse pressure gradient follows that may cause shock induced separation. On the other hand, below 40 percent of span, the shock is formed inside the supersonic bubble and does not reach the pressure side of the adjacent blade as seen in the Mach contour at the 70 percent span in Figure 6. Therefore, only the lower shock oscillation works to dampen the blade motion. On the pressure side, there is a positive peak near the location of the in-passage shock. Notice that the regions of exciting and damping forces on both sides are quite close. That means that the exciting and damping forces generated by the shock oscillations cancel each other resulting in an almost constant amplitude motion.

Figure 9 shows the mean pressure and work per cycle distributions for $\sigma = 180$ degrees. For this case, the two kinds of shocks on the suction side mentioned above do not form a continuous line. The forces generated by the in-passage shock on the suction side works to dampen the oscillation, while the peaks along the oblique shock and the lower part of the shock disappear. Is is also obvious that the in-passage shock on the pressure side produces the large damping forces toward the leading edge. Thus the physical mechanism for the decaying response of $\sigma = 180$ case is clearly explained. A similar argument about the trend with inter-blade phase angle is reported by Chuang (Chuang, 1998) for the same configuration though the reduced frequency is different.

## CONCLUSIONS

A fluid/structure coupled aeroelastic solver for turbomachinery based on a three-dimensional unsteady RANS solver, TFLO, and the finite element structural analysis package, MSC/NASTRAN, was developed for use in turbomachinery flut-ter simulations. The capabilities of the solver are demonstrated

<figure><img src="image_9.png" /><figcaption>Figure 8. MEAN PRESSURE AND WORK PER CYCLE DISTRIBUTION ON THE BLADE SURFACES FOR NEAR PEAK EFFICIENCY, $\sigma = 0$ deg</figcaption></figure>

by applying it to investigate the aeroelastic response characteristics of a transonic fan rotor, NASA Rotor 67. The aeroelastic solver successfully produces different time history of blade displacements for turbomachinery depending on conditions. Results presented for Rotor 67 revealed that the main contribution to the stability of a transonic fan is the unsteady forces generated by the shock motions. In determining these shock motions, inter-blade phase angles play the most important role.

One of the biggest concern for these fluid/structure computations is that the coupled procedure is still computationally expensive. Advanced engines, however, will likely require aeroelastic analysis for multiple blade rows, or a whole wheel system to admit all possible frequencies and inter-blade phase angles. Fluid/structure couplings such as presented here would contribute to these kind of numerical predictions for aeroelastically severe conditions. Furthermore, one of the biggest benefit of the numerical analysis is its richness in data that the experiment has been lacking due the difficulties in installing the measurement instruments.
<figure><img src="image_10.png" /><figcaption>Figure 9. MEAN PRESSURE AND WORK PER CYCLE DISTRIBUTION ON THE BLADE SURFACES FOR NEAR PEAK EFFICIENCY, $\sigma = 180^\circ \text{deg}$</figcaption></figure>

## REFERENCES

Alonso, J. J., Parallel Computations of Unsteady and Aeroelastic Flow Using an Implicit Multigrid-Driven Algorithm, Ph.D. Thesis, Princeton University, 1997

Arnone, A., Viscous Analysis of Three-Dimensional Rotor Flow Using a Multigrid Method, *Journal of Turbomachinery*, Vol.116, No.3, 1994

Bakhle, M. A., Srivastava R., Keith Jr., T. G. and Stefko, G. L., A 3D Euler/Navier-Stokes Aeroelastic Code for Propulsion Applications, AIAA-97-2749, 1997

Brown, S. A., Displacement extrapolations for CFD+CSM aeroelastic analysis, AIAA-97-1090, 1997

Chew, J. W., Marshall, J. G., Vahdati, M. and Imregun, M., Part-Speed Flutter Analysis of a Wide-Chord Fan Blade, *Unsteady Aerodynamics and Aeroelasticity of Turbomachines*, Kluwer Academic Publisher, 1998

Chima, R. V., Viscous Three-Dimensional Calculations of Transonic Fan Performance, NASA-TM-103800, 1991

Chuang, H. A. and Verdon, J. M., A Numerical Simulator for Three-Dimensional Flows Through Vibrating Blade Rows, NASA-CR-1998-208511, 1998

Gerolymos, G. A., Advances in the Numerical Integration of the Three-Dimensional Euler Equations in Vibrating Cascades, *Journal of Turbomachinery*, Vol.115, No.4, 1993

Hall, K. C. and Lorence, C. B., Calculation of Three Dimensional Unsteady Flows in Turbomachinery Using the Linearized Harmonic Euler Equations, *Journal of Turbomachinery*, Vol.115, No.4, 1993

He, L. and Denton, J. D., Three-Dimensional Time-Marching Inviscid and Viscous Solutions for Unsteady Flows Around Vibrating Blades, *Journal of Turbomachinery*, Vol.116, No.3, 1994

Herting, D. N., MSC/NASTRAN Advanced Dynamic Analysis User's Guide, MacNeal-Schwendler Corporation, 1997

Jameson, A., Time Dependant Calculations Using Multigrid with Applications to Unsteady Flows Past Airfoils and Wings, AIAA-91-1596, 1991

Jennions, I. K. and Turner, M. G., Three-Dimensional Navier-Stokes Computations of Transonic Fan Flow Using an Explicit Flow Solver and an Implicit $k-\omega$ Solver, *Journal of Turbomachinery*, Vol.115, No.1, 1993

Ji, S. and Liu, F., Flutter Computation of Turbomachinery Cascades Using a Parallel Unsteady Navier-Stokes Code, *AIAA Journal*, Vol.37, No.3, 1999

Ning, W. and He, L., Computation of Unsteady Flows Around Oscillating Blades Using Linear and Nonlinear Harmonic Euler Methods, *Journal of Turbomachinery*, Vol.120, No.3, 1998

Reuther, J. J., Jameson, A., Farmer, J., Martinelli, L. and Saunders, D., Aerodynamic Shape Optimization of Complex Aircraft Configurations via an Adjoint Formulation, AIAA-96-0094, 1996

Reuther, J. J., Alonso, J. J., Vassberg, J. C., Jameson, A. and Martinelli, L., An Efficient Multiblock Method for Aerodynamic Analysis and Design on Distributed Memory Systems, AIAA-97-1893, 1997

Reuther, J. J., Alonso, J. J., Martins, J. R. R. A. and Smith, S. C., A Coupled Aero-Structural Optimization Method For Complete Aircraft Configurations, AIAA-99-0187, 1999

Sitton, G., MSC/NASTRAN Basic Dynamic Analysis User's Guide, MacNeal-Schwendler Corporation, 1997

Strazisar, A. J., Wood, J. R., Hathaway, M. D. and Suder, K. L., Laser Anemometer Measurements in a Transonic Axial-Flow Fan Rotor, NASA-TP-2879, 1989

Vahdati, M. and Imregun, M., Non-linear Aeroelasticity Analysis Using Unstructured Dynamics Meshes, *Unsteady Aerodynamics and Aeroelasticity of Turbomachines*, Elsevier, 1995

Yao, J., Jameson, A., Alonso, J. J. and Liu, F., Development and Validation of a Massively Parallel Flow Solver for Turbomachinery Flows, AIAA-00-0882, 2000

Yao, J., Davis, R. L., Alonso, J. J. and Jameson, A., Unsteady Flow Investigations in an Axial Turbine Using the Massively Parallel Flow Solver TFLO, AIAA-2001-0529, 2001