**Fluid-Structure Interaction**

Stoia-Djeska Marius and Safta Carmen-Anca *POLITEHNICA University of Bucharest, Bucharest Romania* 

#### **1. Introduction**

194 Fluid Dynamics, Computational Modeling and Applications

[5] Davis, S. J., Allen, M. G., Kessler, W. J., McManus, K. R., Miller, M. F., and Muihall, P. A.,

[6] Davis, S. J., Kessler, W. J., Bachmann, M., and Mulhall, P.A., "Collisional broadening

[7] Davis, S. J., Kessler, W. J., and Keating, P. B., "Progress in the development of sensors for

[8] Nikolaev, V.D., Zagidullin, M. V., Svistun, M. I., Anderson, B. T., Tate, R. F., and Hager,

[9] Madden, T. J. and Miller, J. H., "Simulation of Flow Unsteadiness in Chemical Laser

[10] Fric, T. F., and Roshko, A., "Vortical structure in the wake of a transverse jet," *J. Fluid* 

[12] R. F. Heidner III, C. E. Gardner, G. I. Segal, and T. M. El-Sayed, "Chain-Reaction

[13] Azyazov, V. N.; Pichugin, S. Yu.; Heaven, M. C., "On the dissociation of I2 by O2(a 1Δ):

[16] Madden, T. J. and Solomon, W. C., AIAA 97-2387, *28th Plasmadyamics and Lasers* 

[17] Madden, T. J., *SPIE Proceedings of XIVth International Symposium On Gas Flow & Chemical Lasers and High Power Laser Conference*, Wrocław, Poland, 25-30 August, 2002. [18] Sujudi, D. and Haimes, R., "Identification of Swirling Flow in 3-D Vector Fields," AIAA

[20] Levenberg, K. "A Method for the Solution of Certain Problems in Least Squares," *Quart.* 

[21] Marquardt, D. "An Algorithm for Least-Squares Estimation of Nonlinear Parameters"

*and Lasers Conference*, Toronto, Ontario, Canada, June 6-9, 2005.

*and Intense Beam Control and Applications*, 3931, pp. 156-161, 2000.

CA, 1996.

218-226, San Jose, CA, 1998.

*Electronics*, 38, no. 5, May 2002.

[11] Perram, G. P, .*Int. J. Chem. Kinet*. 27, 817-28 (1995).

*Conference*, Atlanta, GA, June 23-25, 1997.

*SIAM J. Appl. Math*. 11, pp. 431-441, 1963. [22] Keating, P. B., private communication, Sept. 1998.

San Diego CA, 19-22 June 1995.

*Appl. Math*. 2, pp. 164-168, 1944.

[19] www.netlib.org/minpack.

Reno, NV, 5-8 Jan, 2004.

*Mech*, 279, pp. 1-47, 1994.

*Chemistry*, 87, 2348 (1983).

"Diode laser-based sensors for chemical oxygen iodine lasers," Paper 2702-18, *Proceedings of SPIE Conference on Gas and Chemical Lasers*, 2702, pp. 195-201, San Jose,

coefficients for oxygen and water absorption lines used in COIL diagnostics," Paper 3268-80, *Proceedings of SPIE Conference on Gas and ChemicalLasers*, 3268, pp.

COIL devices," *Proceedings of SPIE Conference on Gas, Chemical, and Electrical Lasers* 

G. D., "Results of Small-Signal Gain Measurements on a Supersonic Chemical Oxygen Iodine Laser with an Advanced Nozzle Bank." *IEEE Journal of Quantum* 

Flowfields," AIAA-2004-0805. *42nd AIAA Aerospace Sciences Meeting and Exhibit*,

Mechanism for I2 Dissociation In the O2(1-Delta)-I Atom Laser," *Journal of Physical* 

Pathways involving the excited species I2(A' <sup>3</sup>Π2u,A 3Π1u), I2(X 1∑,υ), and O2(a 1Δ,υ)," Journal of Chemical Physics, Volume 130, Issue 10, pp. 104306-104306-9 (2009). [14] Waichman, K., Barmashenko, B. D. and Rosenwaks, S. "Comparing modeling and

measurements of the output power in chemical oxygen-iodine lasers: A stringent test of I2 dissociation mechanisms," Journal of Chemical Physics, 133, (2010). [15] Madden, T. J. "An Analysis of Mechanisms Underlying Flow Unsteadiness in Chemical

Oxygen-Iodine Laser Mixing Systems," AIAA-2005-5390, *36th AIAA Plasmadynamics* 

Paper 95-1715, *12th AIAA Computational Fluid Dynamics Conference and Open Forum*,

The fluid-structure interaction (FSI) phenomenon is the result of the interactions of multiple continuum fields. The fluid (gas or/and liquid) forces act on a neighboring elastic solid, which is deformed and thus influences the flow of the adhering fluid. Due to the deformation of the solid, both the fluid velocity and the fluid domain change. Usually, this interaction takes place in the presence of other mechanical fields and as important examples we mention here the external body forces and acoustic fields. The mathematical modeling of the FSI phenomenon requires mathematical models for each field implied in and also the coupling mechanisms between the fields. For instance, one can imagine a mathematical model of a FSI phenomenon where the fluid flow field is modeled with the unsteady Navier-Stokes equations, the solid elastic field with the Navier equations while for the acoustic field the inhomogeneous wave equation based on Lighthill's analogy can be used and the external body forces are due to gravity. The boundary conditions for the fluid, solid and acoustic models represent here the coupling mechanism between the fields. However, even everything looks simple and clear the accurate and efficient solution of FSI problems is still a highly complicated task and an open area of research. Furthermore, using complicated models as those enumerated above the only way to obtain a solution to a specific FSI problem is the numerical one. The Computational Fluid Dynamics (CFD), Computational Structural Dynamics (CSD), Finite Element Method (FEM) provide us with specific **mathematical models** and **numerical techniques** that can be coupled to build a numerical solver for a FSI problem. However, even for simple and without engineering relevance problems the computational effort required to solve a specific problem becomes huge if the **mathematical modeling of the physical phenomenon** is not carefully done. As in many other multiple field interaction problems the **level of accuracy** of the mathematical model is crucial for both the efficiency and quality of the results. Further, we notice also that the development of in-house codes for such purposes is a difficult task.

Today we know from the experience earned from a century of practice that different levels of approximation can be used successfully to build up a simplified mathematical model for a FSI problem. This can be done by simplifying the flow equations and by choosing engineering models for the elastic structure. With appropriate coupling mechanisms one can thus obtain mathematical models which are well-matched for practical purposes, starting from the design and up to advanced (flow and/or structural dynamics) control purposes. The use of different simplifying hypothesis (about the dimensionality of the flow and structural models, for instance) requires from the engineer a deep understanding of the multifield interaction phenomenon so that by neglecting different components of the

Fluid-Structure Interaction 197

explicit analytical form and is useful in preparing the computational approach for the real

**Flow induced vibrations** in heat exchanger tubes is another example of FSI, (Mitra, 2009). The induced vibration in these structures is caused by vortex shedding behind the bluff body of the heat exchanger tubes, the fluid-elastic instability and turbulent buffeting. The experiments shows that fully flexible arrays become unstable at a lower flow velocity and those tubes are more stable in steam-water flow as compared to air-water flow. It was found that the boiling water could have a stabilizing effect on fluid-structure instability. **Fluttering and galloping** of lightweight structures are typical results of the fluid structure interaction. The flow regime is important in the analysis of these phenomena, (Barrero-Gil, 2009). Numerical simulations of a transverse galloping phenomenon of a square cylinder at low Reynolds numbers confirm the possibility of galloping with no hysteresis for Re 159 . Large domains of technical problems are dominated by **vortex-induced vibrations** (VIV) as an effect of fluid-structure interaction phenomenon, (Galvao, 2008). The controlling of the wake behind a bluff body means to eliminate unsteady transverse loads while reducing drag. Experimental procedure and numerical simulation demonstrated that it could be obtained a potential flow in the wake behind the bluff body if combinations of flow directing hydrofoils are attached behind a circular cylinder. Two and four symmetric hydrofoils, and a triangular fairing attachment configuration are studied as a **passive control** method with the consequences of no VIV and small drag forces on the cylinder. The stability behaviour of a wake flow pattern is influenced by Reynolds number, turbulence intensity, aspect ratio, end effect, wall proximity, (Kuo, 2009). **Active and passive control** methods could be used to control the flow in the wake behind a bluff body so that could be reduced the form drag, could be suppress vortex shedding or could be changed the heat transfer characteristics. For example, two small control cylinders, with diameter ratio *d/D* = 0.25 are placed symmetrically along the separating shear layers at various stream locations to control the wake behind a circular cylinder in uniform flow at Re 80 *<sup>D</sup>* . Further, in (Pereira Gomes, 2011) the flow– structure interaction a reference test case involving the coupling of unsteady fluid flow and structure motion is studied. It was considered the structure of an aluminium front **cylindrical body** with an attached elastic thin metal plate including a rear mass at the trailing edge. The structure is fixed with one rotational degree of freedom located in the centre of the model front cylinder. The structure is such designed as to attain a self- exciting periodical swivelling movement when exposed to a uniform laminar flow with a Reynolds number up to 270. Good results in the time-phase space were obtained regarding the reproducibility of the coupled fluid–structure motion. An experimental approach of a **structure array of cantilever beams** is demonstrated that neighbouring beams interact through the fluid and their dynamic behaviour is modified, (Kimber, 2009). Aerodynamic interaction between neighbouring cantilever beams operating near their first resonance mode and vibrating at amplitudes was found to be comparable to their widths. Experimental correlations were found to be used to predict the aerodynamic damping in arrays of vibrating cantilevers. In (Facchinetti, 2004) the model of a VIV is described as a one degree of freedom system, elastically supporting a rigid circular cylinder constrained to oscillate transversally to a stationary and uniform flow of free stream velocity. The fluctuating behaviour of vortex street is modelled by a nonlinear oscillator satisfying the van der Pol equation. The wake oscillator was coupled with the structure oscillator and generic forms of coupling have been qualitatively and quantitatively analyzed. The van der Pol wake oscillator model may be extended to 3-D vortex induced vibrations.

structure of a floating airport.

governing equations the simplified mathematical model still captures the essential of the physical phenomenon under concern. Based on this approach the engineers have now at disposal simple mathematical models and even formulae, which offer an insight into a FSI phenomenon or another and also a useful tool to solve a specific FSI problem. This is evident if one takes a look in the field of aeroelasticity, see (Bishplinghoff et all, 1955; Bishplinhoff &Ashley, 1956; Dowell & Ilgamov, 1988).

Historically, the aeroelasticity was the first branch of the FSI which occurred, (Fung, 1956; Dowell, 1975). Aeroelastic and/or hydroelastic vibrations are sustained by air/water forces induced by the moving structure itself. Due to this interaction the elastic structure may suffer large elastic deformations or even divergent oscillations causing failure. This evidence and the concern for very light-weight and thus flexible aircraft and propulsion systems structures have led to the occurrence of the aeroelasticity branch of science from the time when the aeronautical industry was at the beginning. The focus in aeroelasticity was and remains on the safety of the aircraft or civil structures only. However, in the last decades and in parallel with the progress in the computational fluid and structural dynamics and in the computer technology the aeroelasticity shifted forward in the direction of FSI. Currently we are interested in the behavior of both fluid and structural systems and perhaps the best examples come from the field of hydroelasticity and bioengineering applications.

Nowadays, the FSI phenomenon covers a large area of engineering applications and the dedicated literature means thousands of papers. In what follows we present a rather coarse than comprehensive survey of the research work done in different engineering domains to investigate experimentally and numerically different FSI type phenomena of actual interest. The reader is encouraged to search forward in the list of publications for a deeper investigation of the newest results in its domain of work.

A recent application of fluid structure interaction is the **energy harvested** by converting the mechanical strain of an elastic material, under the pressure of a fluid, into electric potential using piezoelectric materials, (Doaré, 2011). For example, a flexible plate fully coupled to a simple dissipative electrical circuit through piezoelectric layers could be an attractive candidate for flow energy harvesting if an axial flow through the plate will produce selfsustained periodic oscillations of the solid body. A global analysis of fluttering modes of a finite- length plate confirmed that waves or modes destabilized by piezoelectric coupling maximize the energy conversion efficiency. Another application is the hydroelastic analysis of very large floating structures used in ocean space (as hydromechanical equipment in the wave's energy conversion process), (Karmakar, 2009). These structures consist of many **articulated elastic plates**, called modules. The articulation of the elastic plates is done by the connectors which depend on the stiffness constants known as the vertical linear spring stiffness and flexural rotational spring stiffness. It is very important to have a stable behaviour of the entire floating structure in any case of water depth (finite depth, infinite depth or shallow water). Assuming that the fluid is inviscid and incompressible, and the motion is nonrotational and simple harmonic in time with angular frequency *ω*, the fluid structure interaction analysis showed that the number of zeros in the reflection coefficient is maximum in the case of infinite water depth and minimum in case of the shallow water approximation. It was observed that the multiple articulated plates were behaved like a single continuous plate if the vertical linear springs and the flexural rotation springs are operating simultaneously. Furthermore, in (Ohkusu, 2004), the fluid structure interaction of **a large and thin floating structure** exposed to the sea wave action is analyzed. An analytical approach to predict vibration of this class of structures was developed. The plate vibration is obtained in an

governing equations the simplified mathematical model still captures the essential of the physical phenomenon under concern. Based on this approach the engineers have now at disposal simple mathematical models and even formulae, which offer an insight into a FSI phenomenon or another and also a useful tool to solve a specific FSI problem. This is evident if one takes a look in the field of aeroelasticity, see (Bishplinghoff et all, 1955;

Historically, the aeroelasticity was the first branch of the FSI which occurred, (Fung, 1956; Dowell, 1975). Aeroelastic and/or hydroelastic vibrations are sustained by air/water forces induced by the moving structure itself. Due to this interaction the elastic structure may suffer large elastic deformations or even divergent oscillations causing failure. This evidence and the concern for very light-weight and thus flexible aircraft and propulsion systems structures have led to the occurrence of the aeroelasticity branch of science from the time when the aeronautical industry was at the beginning. The focus in aeroelasticity was and remains on the safety of the aircraft or civil structures only. However, in the last decades and in parallel with the progress in the computational fluid and structural dynamics and in the computer technology the aeroelasticity shifted forward in the direction of FSI. Currently we are interested in the behavior of both fluid and structural systems and perhaps the best

examples come from the field of hydroelasticity and bioengineering applications.

Nowadays, the FSI phenomenon covers a large area of engineering applications and the dedicated literature means thousands of papers. In what follows we present a rather coarse than comprehensive survey of the research work done in different engineering domains to investigate experimentally and numerically different FSI type phenomena of actual interest. The reader is encouraged to search forward in the list of publications for a deeper

A recent application of fluid structure interaction is the **energy harvested** by converting the mechanical strain of an elastic material, under the pressure of a fluid, into electric potential using piezoelectric materials, (Doaré, 2011). For example, a flexible plate fully coupled to a simple dissipative electrical circuit through piezoelectric layers could be an attractive candidate for flow energy harvesting if an axial flow through the plate will produce selfsustained periodic oscillations of the solid body. A global analysis of fluttering modes of a finite- length plate confirmed that waves or modes destabilized by piezoelectric coupling maximize the energy conversion efficiency. Another application is the hydroelastic analysis of very large floating structures used in ocean space (as hydromechanical equipment in the wave's energy conversion process), (Karmakar, 2009). These structures consist of many **articulated elastic plates**, called modules. The articulation of the elastic plates is done by the connectors which depend on the stiffness constants known as the vertical linear spring stiffness and flexural rotational spring stiffness. It is very important to have a stable behaviour of the entire floating structure in any case of water depth (finite depth, infinite depth or shallow water). Assuming that the fluid is inviscid and incompressible, and the motion is nonrotational and simple harmonic in time with angular frequency *ω*, the fluid structure interaction analysis showed that the number of zeros in the reflection coefficient is maximum in the case of infinite water depth and minimum in case of the shallow water approximation. It was observed that the multiple articulated plates were behaved like a single continuous plate if the vertical linear springs and the flexural rotation springs are operating simultaneously. Furthermore, in (Ohkusu, 2004), the fluid structure interaction of **a large and thin floating structure** exposed to the sea wave action is analyzed. An analytical approach to predict vibration of this class of structures was developed. The plate vibration is obtained in an

Bishplinhoff &Ashley, 1956; Dowell & Ilgamov, 1988).

investigation of the newest results in its domain of work.

explicit analytical form and is useful in preparing the computational approach for the real structure of a floating airport.

**Flow induced vibrations** in heat exchanger tubes is another example of FSI, (Mitra, 2009). The induced vibration in these structures is caused by vortex shedding behind the bluff body of the heat exchanger tubes, the fluid-elastic instability and turbulent buffeting. The experiments shows that fully flexible arrays become unstable at a lower flow velocity and those tubes are more stable in steam-water flow as compared to air-water flow. It was found that the boiling water could have a stabilizing effect on fluid-structure instability. **Fluttering and galloping** of lightweight structures are typical results of the fluid structure interaction. The flow regime is important in the analysis of these phenomena, (Barrero-Gil, 2009). Numerical simulations of a transverse galloping phenomenon of a square cylinder at low Reynolds numbers confirm the possibility of galloping with no hysteresis for Re 159 .

Large domains of technical problems are dominated by **vortex-induced vibrations** (VIV) as an effect of fluid-structure interaction phenomenon, (Galvao, 2008). The controlling of the wake behind a bluff body means to eliminate unsteady transverse loads while reducing drag. Experimental procedure and numerical simulation demonstrated that it could be obtained a potential flow in the wake behind the bluff body if combinations of flow directing hydrofoils are attached behind a circular cylinder. Two and four symmetric hydrofoils, and a triangular fairing attachment configuration are studied as a **passive control** method with the consequences of no VIV and small drag forces on the cylinder. The stability behaviour of a wake flow pattern is influenced by Reynolds number, turbulence intensity, aspect ratio, end effect, wall proximity, (Kuo, 2009). **Active and passive control** methods could be used to control the flow in the wake behind a bluff body so that could be reduced the form drag, could be suppress vortex shedding or could be changed the heat transfer characteristics. For example, two small control cylinders, with diameter ratio *d/D* = 0.25 are placed symmetrically along the separating shear layers at various stream locations to control the wake behind a circular cylinder in uniform flow at Re 80 *<sup>D</sup>* . Further, in (Pereira Gomes, 2011) the flow– structure interaction a reference test case involving the coupling of unsteady fluid flow and structure motion is studied. It was considered the structure of an aluminium front **cylindrical body** with an attached elastic thin metal plate including a rear mass at the trailing edge. The structure is fixed with one rotational degree of freedom located in the centre of the model front cylinder. The structure is such designed as to attain a self- exciting periodical swivelling movement when exposed to a uniform laminar flow with a Reynolds number up to 270. Good results in the time-phase space were obtained regarding the reproducibility of the coupled fluid–structure motion. An experimental approach of a **structure array of cantilever beams** is demonstrated that neighbouring beams interact through the fluid and their dynamic behaviour is modified, (Kimber, 2009). Aerodynamic interaction between neighbouring cantilever beams operating near their first resonance mode and vibrating at amplitudes was found to be comparable to their widths. Experimental correlations were found to be used to predict the aerodynamic damping in arrays of vibrating cantilevers. In (Facchinetti, 2004) the model of a VIV is described as a one degree of freedom system, elastically supporting a rigid circular cylinder constrained to oscillate transversally to a stationary and uniform flow of free stream velocity. The fluctuating behaviour of vortex street is modelled by a nonlinear oscillator satisfying the van der Pol equation. The wake oscillator was coupled with the structure oscillator and generic forms of coupling have been qualitatively and quantitatively analyzed. The van der Pol wake oscillator model may be extended to 3-D vortex induced vibrations.

Fluid-Structure Interaction 199

casing is covered with source panels. A discrete vortex method was used to simulate the flow within a centrifugal pump. The unsteady impeller blade surface forces are estimated. The velocity and pressure fluctuations on the impeller blades are dominated by rotational frequency and its higher harmonics. Further, (Langthjem, 2004b), a hydroacoustic compressible analysis was considered to find estimate the noise in a centrifugal pump. A two-dimensional numerical method to estimate the acoustic pressure fluctuations in a centrifugal pump has been described. The unsteady surface forces which act on the rotor blade are considered to be acoustic dipoles. The discrete vortex method was used to estimate the strengths of the dipoles and to obtain the pressure distribution, and then to obtain the velocity field by applying the unsteady Bernoulli equation. A boundary element method was used to obtain the strengths of the dipoles. The frequency-domain solution is

In (Howell, 2009) it was developed a new computational model of the linear fluid–structure interaction of a **cantilevered-free flexible plate** with an ideal flow in a channel. The transient behaviour of the system is analysed by numerical simulations and a global linearstability map of the system for the infinite-time limit is obtained. The effects of shed vorticity, channel walls, a rigid-inlet surface, temporally varying inlet flow-velocity, and variable plate stiffness were investigated and the flutter instability dependence upon system configuration was found. For example, a short flexible plate (in standard configuration) is destabilised by single-mode flutter caused by an irreversible energy transfer from fluid to structure that principally occurs over the middle part of the flexible plate. A long flexible plate is destabilised by a modal-coalescence flutter and the region of the plate where most

Aeroelastic design considerations related to **long-span bridges** and VIV are also described in (Frandsen, 2004). Governing design criteria for long-span bridges involve the aeroelastic phenomena of vortex-induced oscillations, flutter and buffeting. So, an acceptable flutter limit is one of the principal design criteria for long-span bridges. Scanlan's linearized theory assumes the prescribed motions and is used to estimate the flutter derivatives. Theodorsen's inviscid flat-plate theory is used to estimate the flutter derivatives, too. This study attempts to investigate the use of a coupled fluid-structure interaction finite element solver applied to a long-span bridge. Aerodynamics effects of bridge flutter are investigated using fluid and structural two-dimension finite elements on moving non adaptive grids. The moving interface between fluid and structure is modelled through the arbitrary Lagrangian-Eulerian formulation. It is verified that the flat-plat theory of Theodorsen gives comparatively accurate solutions despite of inviscid flow hypothesis. Also, the prediction of flutter instability for sharp edge bridge deck does not appear sensitive to turbulence and three-dimensional modelling. In the field of aeronautical applications, the fundamental research in the field of FSI continues to go forward. For instance, in (Dessi, 2004) a three-degree-of-freedom of an airfoil with a **control surface** is the physical model which includes different. A 2-D incompressible potential flow has been considered in the model. A standard Runge–Kutta algorithm in conjunction with a 'shooting method' was used to numerically integrate the governing equation and a **stable and unstable limit cycle oscillation** was obtained. The amplitudes and frequencies of limit cycles dependences on the flow speed *V*∞ are obtained. The terms, from the normal-form equations, which are essentially responsible for the nonlinear system behaviour are identified. Using a discrete gust model, Dessi identified and analyzed the damped or undamped **wing oscillations** for different gust's parameters, i.e., intensity and gradient, (Dessi, 2008). Stability analysis has been carried out on a simplified aeroelastic

useful in design optimization so that the flow-noise is minimized.

destabilising energy transfer occurs mainly downstream half.

Applications **in medicine** regarding the interaction between of a rigid, spherical cancer cell with a deformable white blood cell were simulated by developing a quasi-steady technique, (Hoskins, 2009). The six-degree- of-freedom motion, adhesion kinetics, structural mechanics, and fluid dynamics governing equations were all solved, in the context of an octree-based adaptive mesh. The Lagrangian approach was used in grid generation to simulate a cellular system. Weakly coupled fluid-structure interaction models were used for the analysis of the periodic unsteady incompressible flow inside compliant vessels and to simulate **the blood flow in arteries,** (Beulen, 2009). Here it was successfully applied a time- periodic method which proved to have a far better computational stability than the weakly coupled methods based on time step-wise coupling. The method was applied to straight, curved and bifurcating vessels geometries.

Commonly, the coupled fluid–structure dynamic analysis problem in **turbomachinery** could be performed by considering the two indirect coupling methods: a) one based on the cyclic symmetry properties of both structure and fluid, or b) the uncoupled approach which assumes that there is no aerodynamic coupling between the modes and aerodynamic forces, (Tran, 2009). A multi-parameter minimum state modelling method was developed using the spline approximation and the minimum state modelling applied to a numerical model of an aircraft engine compressor disk. In this manner, the number and the cost of the aerodynamic computations in the solutions of the aeroelastic systems were reduced. Further, in (Gnesin, 2004), a fluid structure interaction problem is analyzed as the interaction of aerodynamic, inertial and elastic forces acting on **the blades of the rotor of a turbine stage**. Mathematical model and numerical results for aeroelastic behaviour of a steam turbine last stage with rotor blades of 760 mm are presented. The algorithm proposed involves the coupled solution of 3D unsteady flow through a turbine stage and the dynamics problem for rotor-blade motion by the action of aerodynamic forces, without separating the outer and inner flow fluctuations. The unsteady Euler equations described the 3-D transonic gas flow through the stator and rotor blades in relative motion. An explicit monotonous finite volume difference scheme Godunov-Kolgan scheme was used. Modal approach and a 3-D finite element model of a blade are used in the structural analysis. It was obtained the high-frequency harmonics which are corresponding to the rotor moving past one stator blade pitch. Low-frequency harmonics are caused by blade oscillations and flow non-uniformity downstream from the blade row. Based on experimental **flutter stability** data acquired and using a new model of actuator disk approach, a **reduced-order model of a compressor** was analyzed in (Copeland, 2004). Here it was considered the actuator disk model for the blade rows. The control volume approach was used to determine the aerodynamic force on the flexible blades. Numerical simulations were used to analyze the dynamic behaviour of the compressor and the predicted flutter and stall stability boundaries were proposed.

A numerical method and a procedure to calculate the flow-induced noise in a **centrifugal pump** is presented in (Langthjem, 2004a). First a hydrodynamic incompressible analysis is done to obtain the "background-flow" and noise generating fluid forces. The flow-induced noise in a centrifugal pump was calculated using a computationally method to estimate the noise generating "background-flow". The analysis is restricted to a two-dimension formulation. The fluid structure interaction between the fluid, the rotating blades and the volute tongue, and the interaction between the fluid and the rotor alone are the causes of pressure fluctuations and correspond to dipole sources. Using an acoustic analogy the geometry of the pump was assimilated with a point source for the pump's inlet, the blades of the impeller are covered with vortex elements with discrete, bound vortices and the

Applications **in medicine** regarding the interaction between of a rigid, spherical cancer cell with a deformable white blood cell were simulated by developing a quasi-steady technique, (Hoskins, 2009). The six-degree- of-freedom motion, adhesion kinetics, structural mechanics, and fluid dynamics governing equations were all solved, in the context of an octree-based adaptive mesh. The Lagrangian approach was used in grid generation to simulate a cellular system. Weakly coupled fluid-structure interaction models were used for the analysis of the periodic unsteady incompressible flow inside compliant vessels and to simulate **the blood flow in arteries,** (Beulen, 2009). Here it was successfully applied a time- periodic method which proved to have a far better computational stability than the weakly coupled methods based on time step-wise coupling. The method was applied to straight, curved and

Commonly, the coupled fluid–structure dynamic analysis problem in **turbomachinery** could be performed by considering the two indirect coupling methods: a) one based on the cyclic symmetry properties of both structure and fluid, or b) the uncoupled approach which assumes that there is no aerodynamic coupling between the modes and aerodynamic forces, (Tran, 2009). A multi-parameter minimum state modelling method was developed using the spline approximation and the minimum state modelling applied to a numerical model of an aircraft engine compressor disk. In this manner, the number and the cost of the aerodynamic computations in the solutions of the aeroelastic systems were reduced. Further, in (Gnesin, 2004), a fluid structure interaction problem is analyzed as the interaction of aerodynamic, inertial and elastic forces acting on **the blades of the rotor of a turbine stage**. Mathematical model and numerical results for aeroelastic behaviour of a steam turbine last stage with rotor blades of 760 mm are presented. The algorithm proposed involves the coupled solution of 3D unsteady flow through a turbine stage and the dynamics problem for rotor-blade motion by the action of aerodynamic forces, without separating the outer and inner flow fluctuations. The unsteady Euler equations described the 3-D transonic gas flow through the stator and rotor blades in relative motion. An explicit monotonous finite volume difference scheme Godunov-Kolgan scheme was used. Modal approach and a 3-D finite element model of a blade are used in the structural analysis. It was obtained the high-frequency harmonics which are corresponding to the rotor moving past one stator blade pitch. Low-frequency harmonics are caused by blade oscillations and flow non-uniformity downstream from the blade row. Based on experimental **flutter stability** data acquired and using a new model of actuator disk approach, a **reduced-order model of a compressor** was analyzed in (Copeland, 2004). Here it was considered the actuator disk model for the blade rows. The control volume approach was used to determine the aerodynamic force on the flexible blades. Numerical simulations were used to analyze the dynamic behaviour of the compressor and the predicted flutter and stall

A numerical method and a procedure to calculate the flow-induced noise in a **centrifugal pump** is presented in (Langthjem, 2004a). First a hydrodynamic incompressible analysis is done to obtain the "background-flow" and noise generating fluid forces. The flow-induced noise in a centrifugal pump was calculated using a computationally method to estimate the noise generating "background-flow". The analysis is restricted to a two-dimension formulation. The fluid structure interaction between the fluid, the rotating blades and the volute tongue, and the interaction between the fluid and the rotor alone are the causes of pressure fluctuations and correspond to dipole sources. Using an acoustic analogy the geometry of the pump was assimilated with a point source for the pump's inlet, the blades of the impeller are covered with vortex elements with discrete, bound vortices and the

bifurcating vessels geometries.

stability boundaries were proposed.

casing is covered with source panels. A discrete vortex method was used to simulate the flow within a centrifugal pump. The unsteady impeller blade surface forces are estimated. The velocity and pressure fluctuations on the impeller blades are dominated by rotational frequency and its higher harmonics. Further, (Langthjem, 2004b), a hydroacoustic compressible analysis was considered to find estimate the noise in a centrifugal pump. A two-dimensional numerical method to estimate the acoustic pressure fluctuations in a centrifugal pump has been described. The unsteady surface forces which act on the rotor blade are considered to be acoustic dipoles. The discrete vortex method was used to estimate the strengths of the dipoles and to obtain the pressure distribution, and then to obtain the velocity field by applying the unsteady Bernoulli equation. A boundary element method was used to obtain the strengths of the dipoles. The frequency-domain solution is useful in design optimization so that the flow-noise is minimized.

In (Howell, 2009) it was developed a new computational model of the linear fluid–structure interaction of a **cantilevered-free flexible plate** with an ideal flow in a channel. The transient behaviour of the system is analysed by numerical simulations and a global linearstability map of the system for the infinite-time limit is obtained. The effects of shed vorticity, channel walls, a rigid-inlet surface, temporally varying inlet flow-velocity, and variable plate stiffness were investigated and the flutter instability dependence upon system configuration was found. For example, a short flexible plate (in standard configuration) is destabilised by single-mode flutter caused by an irreversible energy transfer from fluid to structure that principally occurs over the middle part of the flexible plate. A long flexible plate is destabilised by a modal-coalescence flutter and the region of the plate where most destabilising energy transfer occurs mainly downstream half.

Aeroelastic design considerations related to **long-span bridges** and VIV are also described in (Frandsen, 2004). Governing design criteria for long-span bridges involve the aeroelastic phenomena of vortex-induced oscillations, flutter and buffeting. So, an acceptable flutter limit is one of the principal design criteria for long-span bridges. Scanlan's linearized theory assumes the prescribed motions and is used to estimate the flutter derivatives. Theodorsen's inviscid flat-plate theory is used to estimate the flutter derivatives, too. This study attempts to investigate the use of a coupled fluid-structure interaction finite element solver applied to a long-span bridge. Aerodynamics effects of bridge flutter are investigated using fluid and structural two-dimension finite elements on moving non adaptive grids. The moving interface between fluid and structure is modelled through the arbitrary Lagrangian-Eulerian formulation. It is verified that the flat-plat theory of Theodorsen gives comparatively accurate solutions despite of inviscid flow hypothesis. Also, the prediction of flutter instability for sharp edge bridge deck does not appear sensitive to turbulence and three-dimensional modelling.

In the field of aeronautical applications, the fundamental research in the field of FSI continues to go forward. For instance, in (Dessi, 2004) a three-degree-of-freedom of an airfoil with a **control surface** is the physical model which includes different. A 2-D incompressible potential flow has been considered in the model. A standard Runge–Kutta algorithm in conjunction with a 'shooting method' was used to numerically integrate the governing equation and a **stable and unstable limit cycle oscillation** was obtained. The amplitudes and frequencies of limit cycles dependences on the flow speed *V*∞ are obtained. The terms, from the normal-form equations, which are essentially responsible for the nonlinear system behaviour are identified. Using a discrete gust model, Dessi identified and analyzed the damped or undamped **wing oscillations** for different gust's parameters, i.e., intensity and gradient, (Dessi, 2008). Stability analysis has been carried out on a simplified aeroelastic

Fluid-Structure Interaction 201

through the boundary conditions. There are different levels of accuracy in the modeling of

If natural coordinates are used, then the left-hand side of eq. (1) has a decoupled form. Starting from Eq. (1) and using appropriate initial conditions, the dynamic response of the structure can be calculated using dedicated numerical algorithms like those known as Newmark or Hughes-Hilbert-Taylor. However, for aeroelastic calculations the right-hand side of Eq. (1) represents the generalized aerodynamic forces that depend on the generalized coordinates and can be calculated only by solving the flow problem. The structural problem

*dx Ax F U x x*

 *T*

0

*M K MC* 1 1

*<sup>f</sup>*

*M QpU xx* <sup>1</sup> 0 ,

In classical aeroelastic investigation the force term given by Eq. (5) is linearised and, in some simple but useful cases it can be expressed analytically, (Bishplinghoff et all, 1955; Bishplinhoff & Ashley, 1956; Fung, 1956). Under additional assumptions like harmonic time dependence of the displacements and forces the solution of the classical aeroelastic problem of flutter boundaries calculation is simplified from Eq. (2) to an eigenvalue problem. This is the so-called

flutter analysis in the frequency domain, which is a usual practice in aeroelasticity.

 

*dt*

*A*

*F*

Finally, the Eq. (2) is accompanied by the initial condition:

*<sup>f</sup>*

*I*

*Mq Cq Kq Q* (1)

, (2)

*x q q*, (3)

(5)

*xt x* <sup>0</sup> 0 (6)

(4)

**2.1.1 The computational structural dynamics problem and its use in aeroelasticity**  Usually the equations of structural dynamics are obtained using the Finite Element Method and represented by a system of second order differential equations, (Przeminiecki, 1968; Zienkiewiecs, 1971). The linear behaviour of a structure under the action of time dependent external loads is thus described using a finite element model which, for small deformations

both the structure behavior and the fluid flow.

can be rewritten in the phase-state space as:

where the state variable vector is:

and with the matrix

and the force term:

has the well-known linear form:

model in two cases, without and with gust excitation. In the first case, an approximation of the basin of attraction of stable limit cycles in the space of initial conditions was discussed. In the second case, was identified a critical gust intensity for a given gust gradient.

**Flow control techniques** are currently investigated as a primary tool for the **control of aeroelastic vibrations**. **Active control** of flow separation over an air foil using synthetic jets is an already operational technique, see (You, 2008). An unstructured-grid finite-volume large-eddy solver was used to simulate the behaviour of flow inside the synthetic-jet actuator and the synthetic-jet/cross-flow interaction. Numerical simulations confirm 70% increase in lift coefficient if the flow separation is delayed with a synthetic-jet actuation. In (Levasseur, 2008) it was proposed the **control of cavity** flows with two passive acoustic oscillation suppression devices: the rod-in-cross flow and the flat-top spoiler located in the upstream boundary layer. The average reduction of the total pressure level was similar for both devices and was obtained 3–4 dB reduction. The selected test-case is a cavity of length/depth ratio equal to 5, at Mach *M* 0.85 and Reynolds number of *<sup>L</sup>* <sup>6</sup> Re 7 10 .

Finally, a comprehensive and utilitarian review of experimental and numerical modelling on **vortex-induced vibrations** in the last twenty years, (Sarpkaya, 2004) and on open FSI problems (Païdoussis, 2005), are useful references to understand the consequences of the fluid structures interaction in technical applications.

One of the conclusions coming off from the examples above is that the numerical solution of the FSI problems plays an essential role in any kind of work in this field: research, design, decision make.

#### **2. Numerical solution of FSI problems**

In the classical linear approach the behavior of a fluid-structure system is investigated in the frequency domain, (Bishplinghoff et all, 1955; Fung, 1956). In this case, the calculation of either stability boundaries (flutter problem) or dynamic response of the linear system is separated from the computation of linearized unsteady aerodynamic forces. However, in order to investigate nonlinear fluid-structure interactions the governing equations of structural and fluid motion have to be solved in the time domain, (Alonso &Jameson, 1994; Carstens et al., 2003; Liu et al., 2001). This can be done either with time-staggered algorithms (Piperno et al., 1995; Lesoinne & Farhat, 1996) or with coupled algorithms, (Alonso & Jameson, 1994). In the first case optimized numerical integration methods for each of the structural and fluid dynamics models can be used and thus the two sets of equations are treated separately and integrated in a leap-frog fashion. The drawback of this approach is that one can never have a fully converged fluid-structure system at any one time step. This may cause significant energy errors for large times and thus the prediction of the dynamic response and mainly of dynamic instabilities is uncertain. The full coupling of the flow equations with the structural model can be achieved efficiently using an implicit time advancement scheme based on the dual-time approach. The aeroelastic system is thus fully coupled and at each time step a fully converged solution is obtained. In order to avoid the grid generation at each time step the transpiration type boundary conditions are used for the flow equations.

#### **2.1 The mathematical model and numerical solution of the fluid-structure Interaction problem**

The numerical simulation of Fluid-Structure Interaction (FSI) phenomena is based on a mathematical model that describes the structural and fluid dynamics and their coupling

model in two cases, without and with gust excitation. In the first case, an approximation of the basin of attraction of stable limit cycles in the space of initial conditions was discussed.

**Flow control techniques** are currently investigated as a primary tool for the **control of aeroelastic vibrations**. **Active control** of flow separation over an air foil using synthetic jets is an already operational technique, see (You, 2008). An unstructured-grid finite-volume large-eddy solver was used to simulate the behaviour of flow inside the synthetic-jet actuator and the synthetic-jet/cross-flow interaction. Numerical simulations confirm 70% increase in lift coefficient if the flow separation is delayed with a synthetic-jet actuation. In (Levasseur, 2008) it was proposed the **control of cavity** flows with two passive acoustic oscillation suppression devices: the rod-in-cross flow and the flat-top spoiler located in the upstream boundary layer. The average reduction of the total pressure level was similar for both devices and was obtained 3–4 dB reduction. The selected test-case is a cavity of

Finally, a comprehensive and utilitarian review of experimental and numerical modelling on **vortex-induced vibrations** in the last twenty years, (Sarpkaya, 2004) and on open FSI problems (Païdoussis, 2005), are useful references to understand the consequences of the

One of the conclusions coming off from the examples above is that the numerical solution of the FSI problems plays an essential role in any kind of work in this field: research, design,

In the classical linear approach the behavior of a fluid-structure system is investigated in the frequency domain, (Bishplinghoff et all, 1955; Fung, 1956). In this case, the calculation of either stability boundaries (flutter problem) or dynamic response of the linear system is separated from the computation of linearized unsteady aerodynamic forces. However, in order to investigate nonlinear fluid-structure interactions the governing equations of structural and fluid motion have to be solved in the time domain, (Alonso &Jameson, 1994; Carstens et al., 2003; Liu et al., 2001). This can be done either with time-staggered algorithms (Piperno et al., 1995; Lesoinne & Farhat, 1996) or with coupled algorithms, (Alonso & Jameson, 1994). In the first case optimized numerical integration methods for each of the structural and fluid dynamics models can be used and thus the two sets of equations are treated separately and integrated in a leap-frog fashion. The drawback of this approach is that one can never have a fully converged fluid-structure system at any one time step. This may cause significant energy errors for large times and thus the prediction of the dynamic response and mainly of dynamic instabilities is uncertain. The full coupling of the flow equations with the structural model can be achieved efficiently using an implicit time advancement scheme based on the dual-time approach. The aeroelastic system is thus fully coupled and at each time step a fully converged solution is obtained. In order to avoid the grid generation at each time step the transpiration

**2.1 The mathematical model and numerical solution of the fluid-structure Interaction** 

The numerical simulation of Fluid-Structure Interaction (FSI) phenomena is based on a mathematical model that describes the structural and fluid dynamics and their coupling

<sup>6</sup> Re 7 10 .

In the second case, was identified a critical gust intensity for a given gust gradient.

length/depth ratio equal to 5, at Mach *M* 0.85 and Reynolds number of *<sup>L</sup>*

fluid structures interaction in technical applications.

type boundary conditions are used for the flow equations.

**2. Numerical solution of FSI problems** 

decision make.

**problem** 

through the boundary conditions. There are different levels of accuracy in the modeling of both the structure behavior and the fluid flow.

#### **2.1.1 The computational structural dynamics problem and its use in aeroelasticity**

Usually the equations of structural dynamics are obtained using the Finite Element Method and represented by a system of second order differential equations, (Przeminiecki, 1968; Zienkiewiecs, 1971). The linear behaviour of a structure under the action of time dependent external loads is thus described using a finite element model which, for small deformations has the well-known linear form:

$$
\underline{M}\ddot{\underline{q}} + \underline{C}\dot{\underline{q}} + \underline{K}\underline{q} = \underline{Q} \tag{1}
$$

If natural coordinates are used, then the left-hand side of eq. (1) has a decoupled form. Starting from Eq. (1) and using appropriate initial conditions, the dynamic response of the structure can be calculated using dedicated numerical algorithms like those known as Newmark or Hughes-Hilbert-Taylor. However, for aeroelastic calculations the right-hand side of Eq. (1) represents the generalized aerodynamic forces that depend on the generalized coordinates and can be calculated only by solving the flow problem. The structural problem can be rewritten in the phase-state space as:

$$\frac{d\underline{\chi}}{dt} = \underline{A}\underline{\chi} + \underline{E}\left(\underline{U}\_f\left(\underline{\chi}, \dot{\underline{\chi}}\right)\right) \tag{2}$$

where the state variable vector is:

$$\underline{\mathbf{x}} = \left\{ \underline{q}, \dot{\underline{q}} \right\}^T \tag{3}$$

and with the matrix

$$\underline{A} = \begin{bmatrix} 0 & \underline{I} \\ -\underline{M}^{-1}\underline{K} & -\underline{M}^{-1}\underline{C} \end{bmatrix} \tag{4}$$

and the force term:

$$\underline{F} = \begin{cases} 0\\ \underline{M}^{-1} \underline{Q} \left( p \left( \underline{L}\_f \left( \underline{x} \,\dot{\underline{x}} \right) \right) \right) \end{cases} \tag{5}$$

Finally, the Eq. (2) is accompanied by the initial condition:

$$
\underline{\mathbf{x}}(t=0) = \underline{\mathbf{x}}\_0 \tag{6}
$$

In classical aeroelastic investigation the force term given by Eq. (5) is linearised and, in some simple but useful cases it can be expressed analytically, (Bishplinghoff et all, 1955; Bishplinhoff & Ashley, 1956; Fung, 1956). Under additional assumptions like harmonic time dependence of the displacements and forces the solution of the classical aeroelastic problem of flutter boundaries calculation is simplified from Eq. (2) to an eigenvalue problem. This is the so-called flutter analysis in the frequency domain, which is a usual practice in aeroelasticity.

Fluid-Structure Interaction 203

It is known that for flows at low Mach numbers the Eqs. (7, 8) become difficult to be solved numerically, because of the ill-conditioning that imposes very small time steps to be used. For practical purposes it is therefore preferable to replace the inviscid compressible

> *i i <sup>i</sup> <sup>V</sup> f V gV h V t xyy* <sup>0</sup>

> > 2

and the compressibility factor

*u n v n w n w zt x y zn* , (13)

*dt* , (14)

0, 1, 1, 1 , , 1, 1, 1 ,

,

,

.

be provided. This conservative formulation of the incompressible flow equations allows a numerical treatment similar to that used for the compressible equations. The initial

For both inviscid flow models used in this work, **the coupling of the structural and fluid fields** is assured by the statement that the fluid normal velocity on the moving contact

This physical boundary condition at the moving solid wall is used for all the inviscid flow equations, accompanied by other appropriate boundary conditions imposed on all the boundaries of the computational domain. **The pressure field on the solid boundary** is obtained from the solution of the flow equations and is used in Eq. (2) through the force

There are a lot of numerical methods that can be used to solve the flow equations. Starting from the conservative of the unsteady flow equations, the most common way is to use the finite volume method, (Batina, 1991; Ferziger & Peric, 1999; Frink, 1992; Nkonga & Guillard, 1994). The spatial discretization of the flow equations using a cell-centred Godunov type finite volume method leads to a system of differential equations written in semidiscrete form:

> *<sup>f</sup> f*

where the right hand side is the residual vector that contains the discretised fluxes and includes implicitly the boundary conditions given by Eq. (13). The left hand side contains

*RU xx*

*dU*

the cell-averaged values of the fluid states on a current cell:

(11)

(12)

has to

equations with the inviscid incompressible flow model in the Chorin's formulation:

where the unknowns and the flux terms are:

In Eqs. (11, 12) the pressure is defined by *P p*

condition for the Eqs. (11, 12) is similar to Eq. (10).

term given by Eq. (5).

*V Puvw*

*T*

,

*diag diag*

*<sup>T</sup> <sup>i</sup>*

*f u u P uv uw*

2

*g v uv v P vw*

*h w uw vw w P*

surface *S* equals the wall normal velocity *w zt z S xt <sup>n</sup>* , , , i.e.:

*<sup>T</sup> <sup>i</sup>*

2

*<sup>T</sup> <sup>i</sup>*

2

If the linearization is not possible, than Eq. (5) shows that the structural problem cannot be separated from the flow problem and the FSI has to be investigated in the time-domain. This means that the solution of Eq. (2) is obtained numerically, using an explicit or implicit solver that advances in time the structure state.

However, in the phase-space the numerical solution of Eq. (2) with (6) can be done using the so-called time staggered algorithms, (Piperno et al., 1995; Lesoinne & Farhat, 1996). These allow a formal decoupling of the structural and flow problems through a time step. The main advantage of the time staggered algorithms is that there are now available optimised numerical solvers for each of the two problems, taken individually. Their drawback is that the structure and fluid states are not determined exactly at the same time and thus the energy transfer from fluid to structure and reverse cannot be accurately predicted (Alonso &Jameson, 1994; Carstens et al., 2003). This has a negative impact on the accurate timedomain calculation of flutter like instabilities of the aeroelastic systems.

#### **2.1.2 The computational fluid dynamics problem**

In the last decades, unsteady Computational Fluid Dynamics (CFD) has emerged as the basic approach available for predicting the flow behaviour the case of FSI problems, (Anderson et al., 1984; Ferziger & Peric, 1999; Löhner, 2008).

In the present work the fluid is considered inviscid, in order to simplify the mathematical formulations. For compressible flows with shocks, the Euler equations of gas dynamics written in conservative form and with standard notations are:

$$\frac{\partial \underline{\mathsf{U}}}{\partial t} + \frac{\partial \underline{f}(\underline{\mathsf{U}})}{\partial \mathsf{x}} + \frac{\partial \underline{g}(\underline{\mathsf{U}})}{\partial y} + \frac{\partial \underline{h}(\underline{\mathsf{U}})}{\partial z} = 0 \tag{7}$$

where the unknowns and fluxes are given by:

$$\begin{aligned} \underline{\mathbf{U}} &= \begin{Bmatrix} \rho & \rho u & \rho v & \rho w & \rho E \end{Bmatrix}^T, \\ \underline{f} &= \begin{Bmatrix} \rho u & \rho u^2 + p & \rho uv & \rho uw & \rho uw & \rho u \left(E + \frac{p}{\rho}\right) \end{Bmatrix}^T, \\ \underline{g} &= \begin{Bmatrix} \rho v & \rho uv & \rho v^2 + p & \rho wv & \rho v \left(E + \frac{p}{\rho}\right) \end{Bmatrix}^T, \\ \underline{h} &= \begin{Bmatrix} \rho w & \rho uw & \rho vw & \rho w^2 + p & \rho w \left(E + \frac{p}{\rho}\right) \end{Bmatrix}^T \end{aligned} \tag{8}$$

The pressure is determined through the state equation for a thermodynamic perfect gas:

$$p = \left(\kappa - 1\right) \left[\rho E - \frac{\left(\rho u\right)^2 + \left(\rho v\right)^2 + \left(\rho w\right)^2}{2\rho}\right] \tag{9}$$

The initial condition necessary for the solution of Eqs. (7, 8) states that in the entire computational domain the fluid state is known at the beginning of the simulations, i.e.:

$$
\underline{\mathsf{LI}}(t=0) = \underline{\mathsf{LI}}\_0 \tag{10}
$$

If the linearization is not possible, than Eq. (5) shows that the structural problem cannot be separated from the flow problem and the FSI has to be investigated in the time-domain. This means that the solution of Eq. (2) is obtained numerically, using an explicit or implicit solver

However, in the phase-space the numerical solution of Eq. (2) with (6) can be done using the so-called time staggered algorithms, (Piperno et al., 1995; Lesoinne & Farhat, 1996). These allow a formal decoupling of the structural and flow problems through a time step. The main advantage of the time staggered algorithms is that there are now available optimised numerical solvers for each of the two problems, taken individually. Their drawback is that the structure and fluid states are not determined exactly at the same time and thus the energy transfer from fluid to structure and reverse cannot be accurately predicted (Alonso &Jameson, 1994; Carstens et al., 2003). This has a negative impact on the accurate time-

In the last decades, unsteady Computational Fluid Dynamics (CFD) has emerged as the basic approach available for predicting the flow behaviour the case of FSI problems,

In the present work the fluid is considered inviscid, in order to simplify the mathematical formulations. For compressible flows with shocks, the Euler equations of gas dynamics

> *U f U gU h U tx y z* <sup>0</sup>

*<sup>p</sup> f u u p uv uw u E*

 

 

*<sup>p</sup> g v uv v p wv v E*

2

 

> 

*<sup>p</sup> h w uw vw w p w E*

The pressure is determined through the state equation for a thermodynamic perfect gas:

 

*uvw*

The initial condition necessary for the solution of Eqs. (7, 8) states that in the entire computational domain the fluid state is known at the beginning of the simulations, i.e.:

2

*T*

 

  ,

 

 

> 

22 2

2 

*U uvwE*

 

2

(7)

*T*

,

*T*

*Ut U* <sup>0</sup> 0 (10)

*T*

,

(8)

(9)

domain calculation of flutter like instabilities of the aeroelastic systems.

**2.1.2 The computational fluid dynamics problem** 

where the unknowns and fluxes are given by:

(Anderson et al., 1984; Ferziger & Peric, 1999; Löhner, 2008).

written in conservative form and with standard notations are:

*p E*

1

 

that advances in time the structure state.

It is known that for flows at low Mach numbers the Eqs. (7, 8) become difficult to be solved numerically, because of the ill-conditioning that imposes very small time steps to be used. For practical purposes it is therefore preferable to replace the inviscid compressible equations with the inviscid incompressible flow model in the Chorin's formulation:

$$\underline{\Theta} \stackrel{\scriptstyle \mathcal{O}}{=} \frac{\partial \underline{V}}{\partial t} + \underline{\Phi} \left( \frac{\left\| \underline{f} \underline{f}^{i} (\underline{V}) \right\|}{\partial \underline{\text{x}}} + \frac{\left\| \underline{g} \underline{f}^{i} (\underline{V}) \right\|}{\partial \underline{y}} + \frac{\left\| \underline{h} \underline{f}^{i} (\underline{V}) \right\|}{\partial \underline{y}} \right) = 0 \tag{11}$$

where the unknowns and the flux terms are:

$$\begin{aligned} \underline{V} &= \begin{bmatrix} P & u & v & w \end{bmatrix}^T, \\ \underline{\Theta} &= \text{diag}\begin{bmatrix} 0, & 1, & 1, & 1 \end{bmatrix}; \underline{\Phi} = \text{diag}\begin{bmatrix} \beta^2, & 1, & 1, & 1 \end{bmatrix}, \\ \underline{f}^i &= \begin{bmatrix} u & u^2 + P & uv & uw \end{bmatrix}^T, \\ \underline{g}^i &= \begin{bmatrix} v & uv & v^2 + P & vw \end{bmatrix}^T, \\ \underline{h}^i &= \begin{bmatrix} w & uw & vw & w^2 + P \end{bmatrix}^T. \end{aligned} \tag{12}$$

In Eqs. (11, 12) the pressure is defined by *P p* and the compressibility factor has to be provided. This conservative formulation of the incompressible flow equations allows a numerical treatment similar to that used for the compressible equations. The initial condition for the Eqs. (11, 12) is similar to Eq. (10).

For both inviscid flow models used in this work, **the coupling of the structural and fluid fields** is assured by the statement that the fluid normal velocity on the moving contact surface *S* equals the wall normal velocity *w zt z S xt <sup>n</sup>* , , , i.e.:

$$
\mu \cdot n\_x + \upsilon \cdot n\_y + w \cdot n\_z = w\_n(\underline{z}, t) \tag{13}
$$

This physical boundary condition at the moving solid wall is used for all the inviscid flow equations, accompanied by other appropriate boundary conditions imposed on all the boundaries of the computational domain. **The pressure field on the solid boundary** is obtained from the solution of the flow equations and is used in Eq. (2) through the force term given by Eq. (5).

There are a lot of numerical methods that can be used to solve the flow equations. Starting from the conservative of the unsteady flow equations, the most common way is to use the finite volume method, (Batina, 1991; Ferziger & Peric, 1999; Frink, 1992; Nkonga & Guillard, 1994). The spatial discretization of the flow equations using a cell-centred Godunov type finite volume method leads to a system of differential equations written in semidiscrete form:

$$\frac{d\underline{\Pi}\_f}{dt} = \mathcal{R}\left(\underline{\Pi}\_f\left(\underline{\mathbf{x}}\_\prime \dot{\underline{x}}\right)\right) \tag{14}$$

where the right hand side is the residual vector that contains the discretised fluxes and includes implicitly the boundary conditions given by Eq. (13). The left hand side contains the cell-averaged values of the fluid states on a current cell:

Fluid-Structure Interaction 205

*n nn <sup>n</sup>*

1 1 3 4 <sup>1</sup> <sup>0</sup>

For typical FSI problems Eq. (19) is a huge nonlinear algebraic system that must be solved at each time step. A better way to obtain the new values *<sup>n</sup> w* 1 is to build up and solve a

*not not n n <sup>n</sup>*

The residuals contain the boundary conditions and the force acting on the structure that effectively couple the two problems. The boundary conditions and the force are actualized at every internal step. Thus, the solution represents the fluid and structure states at the new

A four time steps Runge-Kutta scheme in the dual time provides the solution of Eq. (20)

physical time step *t* used in Eq. (19) and from numerical stability reasons should obey a

*t t* \* <sup>2</sup> 3

Due to the fully implicit character of Eq. (19), which is also A-stable, the physical time step is determined from accuracy considerations only. The main advantage of this approach is the explicit advancement in the dual time space, which eliminates the solution of huge algebraic

This work deals with the computation of unsteady flows in moving geometries. Because the position of the solid boundaries, i.e. of the structure determines at least partially the fluid domain boundaries, it becomes necessary to perform the numerical calculation of the fluid flow on a moving mesh. The flow equations must therefore be re-written in an Arbitrary Lagrangian-Eulerian (ALE) formulation and the solution procedure of the FSI problem is coupled to the grid dynamics, (Batina, 1991; Ferziger & Peric, 1999; Lesoinne & Farhat, 1996;

It is generally recognized that it is difficult to generate a structured grid about complex configurations. Furthermore, this difficulty is magnified when the aeroelastic deformation of the aircraft or other displacements of the solid boundaries are considered, since the grid must move to conform to the instantaneous shape of the computational domain. As an alternative, algorithms have been developed that make use of unstructured grids. The unstructured grid methods, therefore, have the advantage (over structured grid methods in) that they can easily

**2.1.5 Problems related to moving grids and the geometric conservation law** 

treat complex geometric configurations as well as complicated flow physics.

\* for each physical time step. The dual time step *t*

3 4 , 2 2 

*w w w wRw w R w*

*w ww R w t*

*w*

*w*

(19)

*t*

(21)

1

(20)

\* depends on

\* :

A second order accurate fully implicit time discretization of Eq. (17) is given by:

2

0,

1 \*\*\* \* \*

system of ordinary differential equations in the dual time *t*

 

*dt*

physical time.

condition like:

systems.

very efficiently for *t*

Trepanier et al., 1993).

\* \* \* \* \* \* \*

*w t w where*

0 ,

*dw Rw t*

*n*

$$\begin{aligned} \underline{\underline{U}}\_f(t) &= \left\langle \underline{\underline{U}}\_e(t) \right| e = 1, \,\underline{M} \right\rangle^T, \,\text{with} \\ \underline{\underline{U}}\_e(t) &= \frac{1}{\underline{\Omega}\_e} \int\_{\underline{\Omega}\_e} \underline{\underline{U}}(t) \cdot d\underline{\Omega} \end{aligned} \tag{15}$$

CFD now makes available a variety of explicit and/or implicit algorithms dedicated to the solution of the unsteady flow problems modelled by Eqs. (14, 15) and with appropriate boundary and initial conditions, (Aftosmis et al., 1995; Alonso & Jameson, 1994; Ferziger & Peric, 1999; Frink, 1992; Singh et al., 1995; Trepanier et al., 1993).

#### **2.1.3 The fluid-structure interaction model**

From Eqs. (2) and (14) the coupled system of FSI equations in semidiscrete form is obtained as:

$$
\frac{d}{dt} \begin{Bmatrix} \underline{\mathbf{x}} \\ \underline{\mathbf{U}}\_f \end{Bmatrix} - \begin{Bmatrix} \underline{\mathbf{A}} \underline{\mathbf{x}} + \underline{\mathbf{F}} \left( \underline{\mathbf{U}}\_f \left( \underline{\mathbf{x}}, \dot{\underline{\mathbf{x}}} \right) \right) \\ R \left( \underline{\mathbf{U}}\_f \left( \underline{\mathbf{x}}, \dot{\underline{\mathbf{x}}} \right) \right) \end{Bmatrix} = \mathbf{0} \tag{16}
$$

The coupling in Eq. (16) is generally nonlinear. This equation clearly shows that the fluid and structure states should be evaluated exactly at the same time. However, there are classical **decoupled** solvers using **time-staggered algorithms**. The most known example is due to (Piperno et al., 1995), who introduced a difference of *t* / 2 between the structure and fluid state calculations.

#### **2.1.4 A fully-implicit numerical scheme**

The direct solution of Eq. (16) by an explicit or implicit time advancement algorithm eliminates the already mentioned drawbacks of the time staggered algorithms.

The coupled FSI equations have at least two important features. The first is the huge number of unknowns to be solved simultaneously. It is therefore a necessity the development and implementation of efficient and less time consuming algorithms for the solution of the coupled equations. The second aspect to be taken into account is that the typical time integration steps imposed by accuracy and/or stability reasons usually differ by order of magnitude between the CSD and CFD problem. For instance, in aeroelastic problems the number of unknowns of the CFD problem considerably overrides the number of unknowns of the CSD problem. Further, the temporal resolution of the CSD problem, where usually only the low frequency vibration modes of the structure are investigated, is much greater than that required by the accurate solution of the CFD problem with a numerically stable method. Therefore, the best choice for the solution of Eqs. (16) seem to be the fully implicit time advancement, coupled with an efficient dual-time stepping for the solution of the nonlinear discrete equations, (Alonso &Jameson, 1994). One can rewrite the coupled fluid-structure system given by Eq. (16) in a more convenient form of a system of ordinary differential equations:

$$\frac{d\underline{w}}{dt} + \underline{R}\_w(\underline{w}, \dot{\underline{w}}) = 0 \tag{17}$$

where the unknowns and residuals are:

$$\underline{\mathbf{u}} = \begin{Bmatrix} \underline{\mathbf{x}} \\ \underline{\mathbf{U}}\_f \end{Bmatrix}, \quad \underline{R}\_w = -\begin{Bmatrix} \underline{\mathbf{A}} \underline{\mathbf{x}} + \underline{\mathbf{F}} \left( \underline{\mathbf{U}}\_f \left( \underline{\mathbf{x}}, \dot{\underline{\mathbf{x}}} \right) \right) \\ \mathcal{R} \left( \underline{\mathbf{U}}\_f \left( \underline{\mathbf{x}}, \dot{\underline{\mathbf{x}}} \right) \right) \end{Bmatrix} \tag{18}$$

*T*

(15)

1, ,

,

0

, 0 (17)

(18)

,

 *f*

,

(16)

 *f*

,

*U t U t e M with*

CFD now makes available a variety of explicit and/or implicit algorithms dedicated to the solution of the unsteady flow problems modelled by Eqs. (14, 15) and with appropriate boundary and initial conditions, (Aftosmis et al., 1995; Alonso & Jameson, 1994; Ferziger &

From Eqs. (2) and (14) the coupled system of FSI equations in semidiscrete form is obtained as:

The coupling in Eq. (16) is generally nonlinear. This equation clearly shows that the fluid and structure states should be evaluated exactly at the same time. However, there are classical **decoupled** solvers using **time-staggered algorithms**. The most known example is due to (Piperno et al., 1995), who introduced a difference of *t* / 2 between the structure

The direct solution of Eq. (16) by an explicit or implicit time advancement algorithm

The coupled FSI equations have at least two important features. The first is the huge number of unknowns to be solved simultaneously. It is therefore a necessity the development and implementation of efficient and less time consuming algorithms for the solution of the coupled equations. The second aspect to be taken into account is that the typical time integration steps imposed by accuracy and/or stability reasons usually differ by order of magnitude between the CSD and CFD problem. For instance, in aeroelastic problems the number of unknowns of the CFD problem considerably overrides the number of unknowns of the CSD problem. Further, the temporal resolution of the CSD problem, where usually only the low frequency vibration modes of the structure are investigated, is much greater than that required by the accurate solution of the CFD problem with a numerically stable method. Therefore, the best choice for the solution of Eqs. (16) seem to be the fully implicit time advancement, coupled with an efficient dual-time stepping for the solution of the nonlinear discrete equations, (Alonso &Jameson, 1994). One can rewrite the coupled fluid-structure system given by Eq. (16)

eliminates the already mentioned drawbacks of the time staggered algorithms.

in a more convenient form of a system of ordinary differential equations:

*dt*

,

*<sup>w</sup> dw R ww*

*<sup>w</sup> <sup>f</sup> <sup>f</sup>*

*w R <sup>U</sup> RU xx*

*x Ax F U x x*

 

*<sup>f</sup> <sup>f</sup>*

*x Ax F U x x d dt <sup>U</sup> RU xx*

 *e*

*U t Ut d*

1

*e*

*f e*

*e*

Peric, 1999; Frink, 1992; Singh et al., 1995; Trepanier et al., 1993).

**2.1.3 The fluid-structure interaction model** 

and fluid state calculations.

**2.1.4 A fully-implicit numerical scheme** 

where the unknowns and residuals are:

A second order accurate fully implicit time discretization of Eq. (17) is given by:

$$\frac{2\underline{\mathbf{w}}^{n+1} - 4\underline{\mathbf{w}}^{n} + \underline{\mathbf{w}}^{n-1}}{2\Delta t} + \underline{R}\_{w} \left( \underline{\underline{w}}^{n+1} \right) = 0 \tag{19}$$

For typical FSI problems Eq. (19) is a huge nonlinear algebraic system that must be solved at each time step. A better way to obtain the new values *<sup>n</sup> w* 1 is to build up and solve a system of ordinary differential equations in the dual time *t* \* :

$$\begin{aligned} \frac{d\overline{w}^\*}{dt} + \underline{R}^\* \left(\underline{w}^\*\right) &= 0, \; t^\* \to \infty \\ \underline{w}^\* \left(t^\* = 0\right) &= \underline{w}^\*, \quad \text{where} \\ \underline{w}^{n+1} \stackrel{\text{net}}{=} \underline{w}^\*, \underline{R}^\* \left(\underline{w}^\*\right) &\stackrel{\text{net}}{=} \underline{b}^\* \underline{w}^\* + R\_w \left(\underline{w}^\*\right) - \frac{4\underline{w}^n - \underline{w}^{n-1}}{2\Delta t} \end{aligned} \tag{20}$$

The residuals contain the boundary conditions and the force acting on the structure that effectively couple the two problems. The boundary conditions and the force are actualized at every internal step. Thus, the solution represents the fluid and structure states at the new physical time.

A four time steps Runge-Kutta scheme in the dual time provides the solution of Eq. (20) very efficiently for *t* \* for each physical time step. The dual time step *t* \* depends on physical time step *t* used in Eq. (19) and from numerical stability reasons should obey a condition like:

$$
\Delta t^\* \le \frac{2}{3} \Delta t \tag{21}
$$

Due to the fully implicit character of Eq. (19), which is also A-stable, the physical time step is determined from accuracy considerations only. The main advantage of this approach is the explicit advancement in the dual time space, which eliminates the solution of huge algebraic systems.

#### **2.1.5 Problems related to moving grids and the geometric conservation law**

This work deals with the computation of unsteady flows in moving geometries. Because the position of the solid boundaries, i.e. of the structure determines at least partially the fluid domain boundaries, it becomes necessary to perform the numerical calculation of the fluid flow on a moving mesh. The flow equations must therefore be re-written in an Arbitrary Lagrangian-Eulerian (ALE) formulation and the solution procedure of the FSI problem is coupled to the grid dynamics, (Batina, 1991; Ferziger & Peric, 1999; Lesoinne & Farhat, 1996; Trepanier et al., 1993).

It is generally recognized that it is difficult to generate a structured grid about complex configurations. Furthermore, this difficulty is magnified when the aeroelastic deformation of the aircraft or other displacements of the solid boundaries are considered, since the grid must move to conform to the instantaneous shape of the computational domain. As an alternative, algorithms have been developed that make use of unstructured grids. The unstructured grid methods, therefore, have the advantage (over structured grid methods in) that they can easily treat complex geometric configurations as well as complicated flow physics.

Fluid-Structure Interaction 207

The integral conservative form of the flow equations can be written in ALE formulation for a

*UdV F U ndS <sup>t</sup>*

Here *n nnn <sup>x</sup> y y* , , is the exterior unit normal on the boundary *t* , the fluid velocity is

is the velocity of the volume boundary. The projections of these velocities on the

*dV ndS*

*<sup>n</sup>* , respectively.

*FU n V U p n n n V n n xyz* 

The equation (26) describes a relationship where the time rate of change of the volume of any cell must be exactly balanced by the volumes swept by its boundaries. In order to avoid errors induced by the moving mesh, the condition (26) needs to be satisfied numerically, in addition to the flow equations. Furthermore, this condition must be solved numerically using the same scheme that is used to advance the conservation laws of the fluid to provide self-consistent solution for the local cell volumes. Thus, once the new positions of the moving nodes of the mesh have been calculated, the geometric conservation law is discretized and used to correlate the local cell volumes at the current time level and the normal velocity of the faces which is used in the normal numerical flux calculation. The GCL and the appropriate calculation of the normal velocities of the faces completely define the

*t*

**2.1.6 Numerical implementation of the physical boundary condition for ALE** 

On the solid moving walls, the boundary conditions like Eq. (13) are enforced by using the classical idea of image cells, (Frink, 1992). In these imaginary cells, a fictitious fluid state is defined by setting the flow variables according to the type of the reflective boundary. Then, the Riemann problem so created is solved with the approximate Riemann solver to compute the normal flux across these boundaries. For a moving solid wall, the relative normal velocity is reflected to ensure the impermeability condition. This gives for the normal

The numerical boundary conditions for the pressure, density and tangent velocity fields are taken as those of the adjacent boundary cell. For a static solid wall, the normal velocity is

The time-dependent domain *t* in (24) can be the entire computational domain as well as an arbitrarily specified control cell. When applied to a moving cell, the condition that (24) preserves the trivial solution of a uniform flow field leads to the integral condition called the

0

0, , , , (25)

(24)

*T*

(26)

*imag V Vu n n wall* 2 (27)

time-dependent bounded domain *t* with a boundary *t* as

where the normal inviscid flux looks like:

boundary normal are denoted with *Vn* and

numerical problems related to the moving mesh.

velocity component in the imaginary cell:

**Geometric Conservation Law**:

**formulations** 

simply negated.

*V* and 

The grid management discussed here follows the work pioneered by (Batina, 1991), which has proposed a dynamic grid algorithm employable for problems where the solid body deforms or performs small displacements. The original unstructured grid is generated corresponding to the initial (undeformed) position of the boundaries. The edges of the grid are considered as a system of interconnected springs. In two dimensions, the unstructured grids typically are made up of triangles. Each edge of each triangle is modeled by a tension spring. The spring stiffness for a given edge *i-j* is taken to be inversely proportional to the length of the edge raised to a power, which is written as

$$K\_{ij} = \frac{1}{\left[\left(x\_i - x\_j\right)^2 + \left(y\_i - y\_j\right)^2\right]^{p\_{ij}/2}}\tag{22}$$

where *pij* is a parameter used to control the stiffness depending on its position relative to the moving solid wall. Usually, *ij p* 1 , but we mention here the fact that we are interested to be able to maintain the grid quality my controlling the nodes displacements via this parameter. Then, by assembling the contributions of all edges, the unknown displacements *x* of the free internal nodes and the known displacements *mb x* of the border corresponding to the moving body or other fixed boundaries are coupled by the finite-element type matrix homogeneous equation:

$$
\begin{bmatrix}
\K & K\_{comp} \\
\K\_{comp}^T & K\_{mb}
\end{bmatrix} \cdot \begin{Bmatrix}
\delta \vec{\mathcal{X}} \\
\delta \vec{\mathcal{X}}\_{mb}
\end{Bmatrix} = 0
\tag{23}
$$

The rigidity matrix in (23) is singular and the indices *coup* indicate coupling terms due to the springs connecting internal nodes with border nodes. The solution *x* of the non-singular algebraic system obtained from (23) can be obtained iteratively, by using several Jacobi iterations. The initial guesses of the displacements are the displacements at the previous time level. An example of deforming grids is given in Figure 1.

Fig. 1. The undeformed (a) and deformed (b) grid around a NACA 0012 airfoil.

The grid management discussed here follows the work pioneered by (Batina, 1991), which has proposed a dynamic grid algorithm employable for problems where the solid body deforms or performs small displacements. The original unstructured grid is generated corresponding to the initial (undeformed) position of the boundaries. The edges of the grid are considered as a system of interconnected springs. In two dimensions, the unstructured grids typically are made up of triangles. Each edge of each triangle is modeled by a tension spring. The spring stiffness for a given edge *i-j* is taken to be inversely proportional to the

> *ij ij <sup>p</sup> ij ij*

<sup>1</sup>

Then, by assembling the contributions of all edges, the unknown displacements

*T*

Fig. 1. The undeformed (a) and deformed (b) grid around a NACA 0012 airfoil.

springs connecting internal nodes with border nodes. The solution

time level. An example of deforming grids is given in Figure 1.

<sup>2</sup> 2 2

*x*

0

(22)

of the border corresponding to the

(23)

 *x*   *x* 

of the non-singular

of the

*xx yy*

 

where *pij* is a parameter used to control the stiffness depending on its position relative to the moving solid wall. Usually, *ij p* 1 , but we mention here the fact that we are interested to be able to maintain the grid quality my controlling the nodes displacements via this parameter.

moving body or other fixed boundaries are coupled by the finite-element type matrix

*coup*

*coup mb mb K K x K K x*

The rigidity matrix in (23) is singular and the indices *coup* indicate coupling terms due to the

algebraic system obtained from (23) can be obtained iteratively, by using several Jacobi iterations. The initial guesses of the displacements are the displacements at the previous

length of the edge raised to a power, which is written as

free internal nodes and the known displacements *mb*

homogeneous equation:

*K*

The integral conservative form of the flow equations can be written in ALE formulation for a time-dependent bounded domain *t* with a boundary *t* as

$$\frac{\partial}{\partial t} \iiint\_{\Omega} \underline{\underline{L}} \, dV + \iint\_{\partial \Omega} \vec{F}(\underline{L}) \cdot \vec{n} \, dS = 0 \tag{24}$$

where the normal inviscid flux looks like:

$$
\vec{F} \left( \underline{\mathsf{U}} \right) \cdot \vec{n} = \left( V\_n - \alpha\_n \right) \underline{\mathsf{U}} + p \cdot \left( \mathbf{0}, n\_{x'}, n\_{y'}, n\_{z'} \, V \right)^T \tag{25}
$$

Here *n nnn <sup>x</sup> y y* , , is the exterior unit normal on the boundary *t* , the fluid velocity is *V* and is the velocity of the volume boundary. The projections of these velocities on the boundary normal are denoted with *Vn* and *<sup>n</sup>* , respectively.

The time-dependent domain *t* in (24) can be the entire computational domain as well as an arbitrarily specified control cell. When applied to a moving cell, the condition that (24) preserves the trivial solution of a uniform flow field leads to the integral condition called the **Geometric Conservation Law**:

$$\frac{\partial}{\partial t} \iiint\_{\Omega} dV = \iint\_{\partial \Omega} \vec{\partial} \cdot \vec{n} dS \tag{26}$$

The equation (26) describes a relationship where the time rate of change of the volume of any cell must be exactly balanced by the volumes swept by its boundaries. In order to avoid errors induced by the moving mesh, the condition (26) needs to be satisfied numerically, in addition to the flow equations. Furthermore, this condition must be solved numerically using the same scheme that is used to advance the conservation laws of the fluid to provide self-consistent solution for the local cell volumes. Thus, once the new positions of the moving nodes of the mesh have been calculated, the geometric conservation law is discretized and used to correlate the local cell volumes at the current time level and the normal velocity of the faces which is used in the normal numerical flux calculation. The GCL and the appropriate calculation of the normal velocities of the faces completely define the numerical problems related to the moving mesh.

#### **2.1.6 Numerical implementation of the physical boundary condition for ALE formulations**

On the solid moving walls, the boundary conditions like Eq. (13) are enforced by using the classical idea of image cells, (Frink, 1992). In these imaginary cells, a fictitious fluid state is defined by setting the flow variables according to the type of the reflective boundary. Then, the Riemann problem so created is solved with the approximate Riemann solver to compute the normal flux across these boundaries. For a moving solid wall, the relative normal velocity is reflected to ensure the impermeability condition. This gives for the normal velocity component in the imaginary cell:

$$V\_{\text{n}}^{\text{imag}} = -V\_{\text{n}} + \mathbf{2} \cdot \boldsymbol{\mu}\_{\text{wall}} \tag{27}$$

The numerical boundary conditions for the pressure, density and tangent velocity fields are taken as those of the adjacent boundary cell. For a static solid wall, the normal velocity is simply negated.

Fluid-Structure Interaction 209

In order to overcome the difficulties due to the use of moving grids and to be able to solve at low costs fluid-structure interaction problems at moderate deformations, the transpiration method represents a useful engineering approach. Its most important advantage is that it does not require to update the computational grid or the flux solvers routines. It only involves modifications of the interface boundary conditions. However, this works well only

The modeling procedure of a given fluid-structure system is required to analyze its dynamic behavior under different functional conditions. Furthermore, the mathematical modeling is a necessary step in the designing a control system with a given purpose. Usually, the FSI mathematical model is nonlinear and comprises independent and dependent variables, the dependent being related to the independent ones through a system of equations and relationships. By linearizing this system around an operating point a linear model can be obtained. This linear model describes the behavior of the FSI system around the equilibrium

Mathematical models include uncertainties, because of the parameters whose actual values are known only approximately or could vary around some reference values, (Cacuci et al., 1980; Cacuci, 1981, 1981). As an example taken from aeroelasticity, we enumerate here the flow speed, the air density, the elastic and inertial characteristics of the elasto-dynamic

The purpose of the local sensitivity analysis is to determine quantitatively the behavior of the system responses locally around a chosen point of the trajectory in the phase-space of parameters and dependent variables. The sensitivities are the derivatives of functional type responses with respect to the parameters. These sensitivities are useful for a better understanding of the behavior of the aeroelastic system. There are three methods that can be

The first method is the most common. Running the code that calculates the system response a second time with an increment of a single parameter provides the sensitivity to that parameter. This method consumes large CPU times because for each sensitivity, one complete new run of the code is required. The second method that is available is called the forward sensitivity analysis method. It requires the construction of the linear forward sensitivity equations that may be derived directly from the original system in a consistent manner. Obviously, the sensitivity of the numerical solution to a parameter variable requires a complete solution of the linear forward sensitivity problem (FSP). The advantage is that this calculation is cheaper than the previous one, because of the linearity of the FSP. The disadvantage is that the conclusions are limited to small variations in parameters, (Appel & Gunzburger, 1997; Cacuci et al., 1980). The third method is called the adjoint sensitivity analysis method (ASAM). It requires the construction of the adjoint sensitivity equations, which may be obtained from the forward sensitivity equations and are also linear. Thus, it provides only linear answers. The advantage is that the adjoint function method is largely cheaper, because it provides the sensitivities to all parameters by solving only one time the adjoint equations. The adjoint method has been originally developed in (Cacuci, 1981), who introduced the basic concepts of a comprehensive sensitivity theory. This theory gives precise mathematical and physical meaning to the concept of the sensitivity of a response to parameters that are functions. The adjoint method has been

**2.2 The sensitivity analysis for fluid-structure Interaction problems** 

and thus can be used for stability as well as for dynamic response analysis.

system and last but not least the aerodynamic derivatives.

for small deformations of the structure.

used to determine the sensitivities.

The ALE formulations studied up to now has two practical drawbacks. First, at each time step anew grid must be built - in some way – on the fluid domain, and the associated grid velocity must be computed. Both grid deformations and velocities fields must follow the structure deformation and be smooth in time and space. Second, the flux vectors are modified by the ALE formulation and thus the corresponding solvers must be changed in depth.

**The transpiration method** is, quite simply, a means by which to trick the flow solver into seeing some sort of deflection in the mesh that is not actually there. If a change in surface normal is known, from structural dynamics solver for example, then this change in normal could be applied directly to the existing fixed grid through a slight modification of the existing surface normal. With transpiration, the nodes and faces affected by a surface deflection simply require a modification of its existing normal. Even though the surface is not actually deflected, all the flow sees is the normal at that particular nodal location, it does not matter what that normal is. Let's consider the case of a moving boundary represented in Figure 2.

Fig. 2. The displacement of a surface point from the initial position (1) to the final position (2) during a time step

From the first order approximation of the normal velocity on the solid surface in the point (2) one obtains:

$$
\overrightarrow{V\_2} \cdot \overrightarrow{n\_2} = \overrightarrow{V\_1} \cdot \overrightarrow{n\_1} + \frac{d\left(\overrightarrow{V\_1} \cdot \overrightarrow{n\_1}\right)}{d\overrightarrow{X}} d\overrightarrow{X} + \mathcal{O}\left(d\overrightarrow{X}^2\right) \tag{28}
$$

Then, since *n*<sup>2</sup> and thus the normal velocity *V n* 2 2 are known from the structural problem, one can derive directly the transpiration velocity that modifies the normal velocity at the fixed point (1):

$$
\overrightarrow{V\_1} \cdot \overrightarrow{n\_1} = V\_T = \overrightarrow{V\_2} \cdot \overrightarrow{n\_2} - \frac{d\left(\overrightarrow{V\_1} \cdot \overrightarrow{n\_1}\right)}{d\overrightarrow{X}} d\overrightarrow{X} \tag{29}
$$

Usually, the right hand side term in Eq. (21) is approximated by:

$$V\_T \approx \overrightarrow{V\_2} \cdot \overrightarrow{n\_2} \tag{30}$$

The ALE formulations studied up to now has two practical drawbacks. First, at each time step anew grid must be built - in some way – on the fluid domain, and the associated grid velocity must be computed. Both grid deformations and velocities fields must follow the structure deformation and be smooth in time and space. Second, the flux vectors are modified by the

**The transpiration method** is, quite simply, a means by which to trick the flow solver into seeing some sort of deflection in the mesh that is not actually there. If a change in surface normal is known, from structural dynamics solver for example, then this change in normal could be applied directly to the existing fixed grid through a slight modification of the existing surface normal. With transpiration, the nodes and faces affected by a surface deflection simply require a modification of its existing normal. Even though the surface is not actually deflected, all the flow sees is the normal at that particular nodal location, it does not matter what that

Fig. 2. The displacement of a surface point from the initial position (1) to the final position

From the first order approximation of the normal velocity on the solid surface in the point

*Vn Vn dX dX d X*

one can derive directly the transpiration velocity that modifies the normal velocity at the

*Vn V Vn d X*

22 11

Usually, the right hand side term in Eq. (21) is approximated by:

*T*

11 22

*dV n*

and thus the normal velocity *V n* 2 2 are known from the structural problem,

1 1 2

*d X* 1 1

*dV n*

(28)

(29)

*V Vn <sup>T</sup>* 2 2 (30)

(2) during a time step

(2) one obtains:

Then, since *n*<sup>2</sup>

fixed point (1):

ALE formulation and thus the corresponding solvers must be changed in depth.

normal is. Let's consider the case of a moving boundary represented in Figure 2.

In order to overcome the difficulties due to the use of moving grids and to be able to solve at low costs fluid-structure interaction problems at moderate deformations, the transpiration method represents a useful engineering approach. Its most important advantage is that it does not require to update the computational grid or the flux solvers routines. It only involves modifications of the interface boundary conditions. However, this works well only for small deformations of the structure.

#### **2.2 The sensitivity analysis for fluid-structure Interaction problems**

The modeling procedure of a given fluid-structure system is required to analyze its dynamic behavior under different functional conditions. Furthermore, the mathematical modeling is a necessary step in the designing a control system with a given purpose. Usually, the FSI mathematical model is nonlinear and comprises independent and dependent variables, the dependent being related to the independent ones through a system of equations and relationships. By linearizing this system around an operating point a linear model can be obtained. This linear model describes the behavior of the FSI system around the equilibrium and thus can be used for stability as well as for dynamic response analysis.

Mathematical models include uncertainties, because of the parameters whose actual values are known only approximately or could vary around some reference values, (Cacuci et al., 1980; Cacuci, 1981, 1981). As an example taken from aeroelasticity, we enumerate here the flow speed, the air density, the elastic and inertial characteristics of the elasto-dynamic system and last but not least the aerodynamic derivatives.

The purpose of the local sensitivity analysis is to determine quantitatively the behavior of the system responses locally around a chosen point of the trajectory in the phase-space of parameters and dependent variables. The sensitivities are the derivatives of functional type responses with respect to the parameters. These sensitivities are useful for a better understanding of the behavior of the aeroelastic system. There are three methods that can be used to determine the sensitivities.

The first method is the most common. Running the code that calculates the system response a second time with an increment of a single parameter provides the sensitivity to that parameter. This method consumes large CPU times because for each sensitivity, one complete new run of the code is required. The second method that is available is called the forward sensitivity analysis method. It requires the construction of the linear forward sensitivity equations that may be derived directly from the original system in a consistent manner. Obviously, the sensitivity of the numerical solution to a parameter variable requires a complete solution of the linear forward sensitivity problem (FSP). The advantage is that this calculation is cheaper than the previous one, because of the linearity of the FSP. The disadvantage is that the conclusions are limited to small variations in parameters, (Appel & Gunzburger, 1997; Cacuci et al., 1980). The third method is called the adjoint sensitivity analysis method (ASAM). It requires the construction of the adjoint sensitivity equations, which may be obtained from the forward sensitivity equations and are also linear. Thus, it provides only linear answers. The advantage is that the adjoint function method is largely cheaper, because it provides the sensitivities to all parameters by solving only one time the adjoint equations. The adjoint method has been originally developed in (Cacuci, 1981), who introduced the basic concepts of a comprehensive sensitivity theory. This theory gives precise mathematical and physical meaning to the concept of the sensitivity of a response to parameters that are functions. The adjoint method has been

Fluid-Structure Interaction 211

definition of the sensitivity of a response to variations in the system parameters is the Gateaux (G) differential. The G-differential *DP S s* , of an operator *P S* at *S* with an

*<sup>d</sup> DP S s P S s*

Applying Eq. (34) to Eq. (32) yields the sensitivity *DR X G x g* , ;, of the response as:

*t t f f X G X G E E DR X G x <sup>g</sup> xdt <sup>g</sup> dt*

0 0

The first and the second terms of the Eq. (35) are customarily called the indirect effect term and the direct effect term, respectively. The exact value of the perturbed response is predicted by the sensitivity to first-order accuracy in parameters and solution variations, i.e.:

*R X x G g R X G DR X G x g x g* <sup>2</sup> <sup>2</sup>

the indirect effect term and thus the desired sensitivity of the response functional can be evaluated only after determining the solution perturbations *x*. Up to first order, the relationship between *x* and *g* is obtained by taking G-differentials of Eqs. (1). This leads to

*oxgXG*

*<sup>f</sup>*

Once the solution of the base-case problem represented by Eq. (31) has been determined, then for a given vector of parameter changes the FSP could be solved to determine the variations in the flow problem solution, *x* . The solution of the FSP is advantageous to employ only if the number of different responses exceeds the number of system parameters. The alternative way, however, is to eliminate the explicit appearance of *x* in Eq. (35). This elimination process relies on the possibility to construct an adjoint operator corresponding to Eq. (37). Using the inner product given by Eq. (33) and integration by parts, one obtains

where *o XG* \* ; , is the adjoint operator to *o x* ,; , *g X G* and *B x* , is the associated

Further, **the explicit appearance** of the unknown *x* in Eq. (35) **is eliminated** if the adjoint

vector function is determined so that it satisfies the linear equation:

0 , 0,

, ; , 0,

*x xt t* <sup>0</sup>

,

, ;,

*d* <sup>0</sup>

*X G*

, , ,; , (36)

*oxgXG xo XG Bx* \* , ,; , , ; , , (38)

*xo XG Bx* \* , ;, , (39)

, ,

(34)

(35)

(37)

increment *s* is defined as

the FSP given by:

from Eqs. (7) for the adjoint vector of *x* :

bilinear form. Eqs. (37) and (38) yield:

extensively used mainly in the field of internal flows simulations, (Cacuci, 1981, 1981; Ounsy et al., 1994; Toomarian et al., 1988).

**Remark.** In parallel with the sensitivity analysis of a response, for stability analysis purposes it is also useful to derive a computational frame for the calculation of the sensitivities of the eigenvalues resulting from eigenvalue problems.

In this paragraph we present the theoretical developments necessary for the application of both the forward and the adjoint sensitivity analysis methods for a typical FSI mathematical model. Starting from the state-space formulation we present the formalism that uses adjoint functions associated to state variables to determine sensitivities of the response to all the parameters. The response considered here is of functional-type and the parameters are represented by the initial conditions of the base-case problem and by the inertial, elastic and aerodynamic coefficients of the model. The adjoint system satisfied by the adjoint functions is determined and shown to be a well-posed linear system of ordinary differential equations subjected to specific initial.

#### **2.2.1 The basics of the adjoint sensitivity analysis method**

For the purposes of sensitivity analysis, let the spatial-discretized FSI mathematical model be represented in operator form as a system of nonlinear time-dependent differential equations:

$$\begin{aligned} \bullet \, \mathcal{O}\left(X, \underline{\mathcal{G}}\right) &= 0, \, t \in \left[0, t\_f\right] \\ X\left(0\right) &= X\_0 \end{aligned} \tag{31}$$

where the state variables are *X t* and *GG G* <sup>1</sup> ,.., *<sup>M</sup>* is the vector containing all the problem parameters. We assume that the well-posed initial value problem (IVP) represented by Eqs. (31) is solved for a set of nominal (or base-case) parameter values *G* and has a unique solution, denoted here by *X*. The response considered in this analysis is a nonlinear functional depending simultaneously on *X* and *G* and generally defined as:

$$R\left(X,\underline{\mathbf{C}}\right) = \int\_0^t E\left[X\left(t\right),\underline{\mathbf{C}}\right]dt\tag{32}$$

where *tf* is some final time value and *E* is a nonlinear function depending on time and the system solution and parameters. The response of the aeroelastic system is the total mechanical energy at the final time, for instance. It is worthwhile to notice here that the integrals appearing in Eq. (32) may be considered as an inner product in the Hilbert spaces of all square integrable vector functions to whom the solution of Eq. (31) belongs, for instance:

$$\left< \alpha, \beta \right> = \bigcup\_{0}^{t\_f} \alpha(t)\beta(t)dt \tag{33}$$

It is obvious that small variations *gg g* <sup>1</sup> ,.., *<sup>M</sup>* in the system parameters induce variations in the solution, *x* , so that the perturbed solution *X x* is satisfying a system similar to Eq. (31). The objective of local sensitivity analysis is to analyze the behavior of the system responses locally around a chosen point *X G*, of trajectory in the combined phase-space of state variables *Q* and parameters *G* . The most general and fundamental concept for the

extensively used mainly in the field of internal flows simulations, (Cacuci, 1981, 1981; Ounsy

**Remark.** In parallel with the sensitivity analysis of a response, for stability analysis purposes it is also useful to derive a computational frame for the calculation of the

In this paragraph we present the theoretical developments necessary for the application of both the forward and the adjoint sensitivity analysis methods for a typical FSI mathematical model. Starting from the state-space formulation we present the formalism that uses adjoint functions associated to state variables to determine sensitivities of the response to all the parameters. The response considered here is of functional-type and the parameters are represented by the initial conditions of the base-case problem and by the inertial, elastic and aerodynamic coefficients of the model. The adjoint system satisfied by the adjoint functions is determined and shown to be a well-posed linear system of ordinary differential equations

For the purposes of sensitivity analysis, let the spatial-discretized FSI mathematical model be represented in operator form as a system of nonlinear time-dependent differential

*XG t tf*

where the state variables are *X t* and *GG G* <sup>1</sup> ,.., *<sup>M</sup>* is the vector containing all the problem parameters. We assume that the well-posed initial value problem (IVP) represented by Eqs. (31) is solved for a set of nominal (or base-case) parameter values *G* and has a unique solution, denoted here by *X*. The response considered in this analysis is a nonlinear

> *f t R X G E X t G dt* 0

> > *f t*

0

 

It is obvious that small variations *gg g* <sup>1</sup> ,.., *<sup>M</sup>* in the system parameters induce variations in the solution, *x* , so that the perturbed solution *X x* is satisfying a system similar to Eq. (31). The objective of local sensitivity analysis is to analyze the behavior of the system responses locally around a chosen point *X G*, of trajectory in the combined phase-space of state variables *Q* and parameters *G* . The most general and fundamental concept for the

where *tf* is some final time value and *E* is a nonlinear function depending on time and the system solution and parameters. The response of the aeroelastic system is the total mechanical energy at the final time, for instance. It is worthwhile to notice here that the integrals appearing in Eq. (32) may be considered as an inner product in the Hilbert spaces of all square

*t t dt*

, 0, 0,

(31)

, , (32)

, (33)

 

*X X*<sup>0</sup>

0

functional depending simultaneously on *X* and *G* and generally defined as:

integrable vector functions to whom the solution of Eq. (31) belongs, for instance:

sensitivities of the eigenvalues resulting from eigenvalue problems.

**2.2.1 The basics of the adjoint sensitivity analysis method** 

et al., 1994; Toomarian et al., 1988).

subjected to specific initial.

equations:

definition of the sensitivity of a response to variations in the system parameters is the Gateaux (G) differential. The G-differential *DP S s* , of an operator *P S* at *S* with an increment *s* is defined as

$$DP\{\underline{S}, \mathbf{s}\} = \left[\frac{d}{d\mathcal{E}}P\{\underline{S} + \mathcal{E} \cdot \mathbf{s}\}\right]\_{\mathcal{E} = 0} \tag{34}$$

Applying Eq. (34) to Eq. (32) yields the sensitivity *DR X G x g* , ;, of the response as:

$$DR\left(X,\underline{G};x,\underline{g}\right) = \int\_{0}^{t} \left[\frac{\partial E}{\partial X}\right]^{(X,\underline{G})} \cdot xdt + \int\_{0}^{t} \left[\frac{\partial E}{\partial \underline{G}}\right]^{(X,\underline{G})} \cdot \underline{g}\,dt\tag{35}$$

The first and the second terms of the Eq. (35) are customarily called the indirect effect term and the direct effect term, respectively. The exact value of the perturbed response is predicted by the sensitivity to first-order accuracy in parameters and solution variations, i.e.:

$$R\left(X+\mathbf{x},\underline{\mathsf{C}}+\underline{\mathsf{g}}\right) - R\left(X,\underline{\mathsf{C}}\right) = DR\left(X,\underline{\mathsf{C}};\mathbf{x},\underline{\mathsf{g}}\right) + \Theta\left(\left\|\mathbf{x}\right\|^2 + \left\|\underline{\mathbf{g}}\right\|^2\right) \tag{36}$$

the indirect effect term and thus the desired sensitivity of the response functional can be evaluated only after determining the solution perturbations *x*. Up to first order, the relationship between *x* and *g* is obtained by taking G-differentials of Eqs. (1). This leads to the FSP given by:

$$\begin{aligned} \sigma\left(\mathbf{x}, \underline{\underline{\mathbf{G}}}; \mathbf{X}, \underline{\underline{\mathbf{G}}}\right) &= \mathbf{0}, \\ \mathbf{x}\left(\mathbf{0}\right) &= \mathbf{x}\_0, \mathbf{t} \in \left[\mathbf{0}, \mathbf{t}\_f\right] \end{aligned} \tag{37}$$

Once the solution of the base-case problem represented by Eq. (31) has been determined, then for a given vector of parameter changes the FSP could be solved to determine the variations in the flow problem solution, *x* . The solution of the FSP is advantageous to employ only if the number of different responses exceeds the number of system parameters. The alternative way, however, is to eliminate the explicit appearance of *x* in Eq. (35). This elimination process relies on the possibility to construct an adjoint operator corresponding to Eq. (37). Using the inner product given by Eq. (33) and integration by parts, one obtains from Eqs. (7) for the adjoint vector of *x* :

$$\left\langle \Phi, o\left(\mathbf{x}, g; X, \underline{\mathbf{C}}\right) \right\rangle = \left\langle \mathbf{x}, o^\*\left(\Phi; X, \underline{\mathbf{C}}\right) \right\rangle + B\left(\mathbf{x}, \Phi\right) \tag{38}$$

where *o XG* \* ; , is the adjoint operator to *o x* ,; , *g X G* and *B x* , is the associated bilinear form. Eqs. (37) and (38) yield:

$$\left< \chi, \rho^\* \left( \Phi; X, \underline{\mathbb{C}} \right) \right> = -B \left( \chi, \Phi \right) \tag{39}$$

Further, **the explicit appearance** of the unknown *x* in Eq. (35) **is eliminated** if the adjoint vector function is determined so that it satisfies the linear equation:

Fluid-Structure Interaction 213

The reader is encouraged to obtain this formula as an exercise. The calculation of the derivative of the state matrix with respect to the parameters is not a simple task and for many situations it must be done numerically. However, the last equation is useful for

In this chapter we first outline the general principles of FSI followed by some examples taken from aeronautical, civil and biomedical engineering. One of the conclusions coming off from these examples is that the numerical solution of the FSI problems plays an essential

Then, we focused on the numerical simulation of the FSI problems using intensive computational techniques. These are based on a mathematical model that describes the coupling of structural and fluid dynamics models through the boundary conditions. There are different levels of accuracy in the modeling of both the structure behavior and the fluid flow. For simplicity reasons only we have chosen the case of inviscid flows interacting with linear structures. In the phase-space the numerical solution can be done using the so-called time staggered algorithms. These allow a formal decoupling of the structural and flow problems through a time step. The main advantage of the time staggered algorithms is that there are now available optimised numerical solvers for each of the two problems, taken individually. Their major drawback is that the structure and fluid states are not determined exactly at the same time and thus the energy transfer from fluid to structure and reverse cannot be accurately predicted. This is why a better way for the solution of the FSI equations is the fully implicit time advancement, coupled with an efficient dual-time stepping for the solution of the nonlinear discrete equations. The main advantage of this approach is the explicit advancement in the dual time space, which eliminates the solution of huge algebraic systems. Aspects related to the treatment of moving grids and the related Geometric Conservation Law and to the imposition of the physical boundary condition on the moving

The last part of the chapter presents in some details the adjoint sensitivity analysis method, which is useful for a better understanding of FSI problems. The purpose of the local sensitivity analysis is to determine quantitatively the behavior of the system responses locally around a chosen point of the trajectory in the phase-space of parameters and dependent variables. The sensitivities are the derivatives of functional type responses with respect to the parameters. The local sensitivity analysis offers the possibility to quantify the role of the parameters on the response of an aeroelastic system. The adjoint sensitivity analysis method is numerically efficient but requires important theoretical developments and these are summarized in this work. We notice that beyond the general form of the equations, the use of the adjoint sensitivity analysis method needs also a completely

For stability analysis purposes the sensitivities of the eigenvalues of a state matrix to the problem's parameters are useful and we presented the way to follow to calculate them. The derivatives of an eigenvalue with respect to the problem's parameters show the importance of those parameters on the stability of the FSI system. The last equation is useful for

different and new computer code to be developed and tested.

understanding, analysis and design of the aeroelastic system.

understanding, analysis and design of the aeroelastic system.

**3. Conclusion** 

role in any kind of work in this field.

wall are also mentioned.

$$\rho^\* \left( \Phi; X, \underline{\mathcal{G}} \right) = \left[ \frac{\partial E}{\partial X} \right]^{(X, \underline{\mathcal{G}})} \tag{40}$$

The solution of the Eq. (40) can be done successfully only if one imposes appropriate initial conditions. These must be set so that all the terms involving unknown values of *x* disappear from Eq. (39). If this is done, then the sensitivity of the response is determined directly using the relationship:

$$DR\left(X,\underline{G};\mathbf{x},\underline{\mathbf{g}}\right) = \tilde{B}\left(\mathbf{x},\boldsymbol{\Phi}\right) + \int\_{0}^{t\_f} \left[\frac{\partial E}{\partial \underline{G}}\right]^{\left(X,\underline{G}\right)} \underline{\mathbf{g}}\,dt\tag{41}$$

In Eq.(41) *B x* , is the remaining part of the bilinear form *B x* , after the elimination of the unknown terms *x t <sup>f</sup>* through the use of the final time conditions for , i.e. *tf* 0 . This condition is transformed into an initial condition using the time variable change *<sup>f</sup> t t* . This leads to a new IVP for the adjoint function :

$$\log^\*\left(\Phi; X, \underline{\mathsf{G}}\right) = \left[\frac{\partial E}{\partial X}\right]^{(X, \underline{\mathsf{G}})}, \quad \Phi\left(\tau\right) \equiv 0 \tag{42}$$

Thus, the adjoint functions are independent of parameter variations. The source term in Eq. (42) depends on the choice of the response. Finally, we notice that the IVP for the adjoint functions is linear and depends on the base-case solution *X*.

The local sensitivity analysis method offers therefore the possibility to quantify the influence of all the parameters on the response of an aeroelastic system. Compared with the direct and forward methods, the adjoint sensitivity analysis method is numerically efficient but requires important theoretical developments and these are summarized in this work. We notice that beyond the general form of the equations, the use of the adjoint sensitivity analysis method needs also a completely different and new computer code to be developed and tested.

#### **2.2.2 The sensitivity of the eigenvalues in local stability analysis**

For local stability analysis in the state-space form the forward linear system (37) or any other linear system has a homogeneous part written like:

$$\begin{cases} \dot{\mathbf{x}} - A \left( X, \underline{\mathbf{G}} \right) \mathbf{x} = \mathbf{0}, \\ \mathbf{x}(\mathbf{0}) = \mathbf{x}\_0 \end{cases} \tag{43}$$

where *A* is the state matrix. This depends on the base case solution and also on the problem's parameters, which are supposed to have small variations around the base-case values. Using the right eigenvectors problems for both the forward and adjoint homogeneous systems, one can proof that the sensitivity of the *r*th eigenvalue *r* with respect the *p*th parameter *gp* can be determined using the following formula:

$$\frac{\partial \Lambda\_r}{\partial \mathbf{g}\_p} = \Psi\_r^T \frac{\partial \mathbf{A}}{\partial \mathbf{g}\_p} \Phi\_r \tag{44}$$

*X G <sup>E</sup> o XG*

The solution of the Eq. (40) can be done successfully only if one imposes appropriate initial conditions. These must be set so that all the terms involving unknown values of *x* disappear from Eq. (39). If this is done, then the sensitivity of the response is determined

*tf X G <sup>E</sup> DR X G x <sup>g</sup> B x <sup>g</sup> dt*

the unknown terms *x t <sup>f</sup>* through the use of the final time conditions for , i.e. *tf* 0 . This condition is transformed into an initial condition using the time variable

*X*

Thus, the adjoint functions are independent of parameter variations. The source term in Eq. (42) depends on the choice of the response. Finally, we notice that the IVP for the adjoint

The local sensitivity analysis method offers therefore the possibility to quantify the influence of all the parameters on the response of an aeroelastic system. Compared with the direct and forward methods, the adjoint sensitivity analysis method is numerically efficient but requires important theoretical developments and these are summarized in this work. We notice that beyond the general form of the equations, the use of the adjoint sensitivity analysis method needs also a completely different and new computer code to

For local stability analysis in the state-space form the forward linear system (37) or any other

where *A* is the state matrix. This depends on the base case solution and also on the problem's parameters, which are supposed to have small variations around the base-case values. Using the right eigenvectors problems for both the forward and adjoint homogeneous systems, one can proof that the sensitivity of the *r*th eigenvalue *r* with

*r T*

*p p g g* <sup>A</sup>

, 0,

*r r*

*x AXGx x x*<sup>0</sup>

respect the *p*th parameter *gp* can be determined using the following formula:

0 

, \* ; , , 0

*X*

,

*G*

0

, is the remaining part of the bilinear form *B x* , after the elimination of

,

(40)

(41)

(42)

(43)

(44)

\* ; ,

, ;, ,

*t t* . This leads to a new IVP for the adjoint function :

*X G <sup>E</sup> o XG*

**2.2.2 The sensitivity of the eigenvalues in local stability analysis** 

linear system has a homogeneous part written like:

functions is linear and depends on the base-case solution *X*.

directly using the relationship:

In Eq.(41) *B x*

change *<sup>f</sup>* 

be developed and tested.

The reader is encouraged to obtain this formula as an exercise. The calculation of the derivative of the state matrix with respect to the parameters is not a simple task and for many situations it must be done numerically. However, the last equation is useful for understanding, analysis and design of the aeroelastic system.

#### **3. Conclusion**

In this chapter we first outline the general principles of FSI followed by some examples taken from aeronautical, civil and biomedical engineering. One of the conclusions coming off from these examples is that the numerical solution of the FSI problems plays an essential role in any kind of work in this field.

Then, we focused on the numerical simulation of the FSI problems using intensive computational techniques. These are based on a mathematical model that describes the coupling of structural and fluid dynamics models through the boundary conditions. There are different levels of accuracy in the modeling of both the structure behavior and the fluid flow. For simplicity reasons only we have chosen the case of inviscid flows interacting with linear structures. In the phase-space the numerical solution can be done using the so-called time staggered algorithms. These allow a formal decoupling of the structural and flow problems through a time step. The main advantage of the time staggered algorithms is that there are now available optimised numerical solvers for each of the two problems, taken individually. Their major drawback is that the structure and fluid states are not determined exactly at the same time and thus the energy transfer from fluid to structure and reverse cannot be accurately predicted. This is why a better way for the solution of the FSI equations is the fully implicit time advancement, coupled with an efficient dual-time stepping for the solution of the nonlinear discrete equations. The main advantage of this approach is the explicit advancement in the dual time space, which eliminates the solution of huge algebraic systems. Aspects related to the treatment of moving grids and the related Geometric Conservation Law and to the imposition of the physical boundary condition on the moving wall are also mentioned.

The last part of the chapter presents in some details the adjoint sensitivity analysis method, which is useful for a better understanding of FSI problems. The purpose of the local sensitivity analysis is to determine quantitatively the behavior of the system responses locally around a chosen point of the trajectory in the phase-space of parameters and dependent variables. The sensitivities are the derivatives of functional type responses with respect to the parameters. The local sensitivity analysis offers the possibility to quantify the role of the parameters on the response of an aeroelastic system. The adjoint sensitivity analysis method is numerically efficient but requires important theoretical developments and these are summarized in this work. We notice that beyond the general form of the equations, the use of the adjoint sensitivity analysis method needs also a completely different and new computer code to be developed and tested.

For stability analysis purposes the sensitivities of the eigenvalues of a state matrix to the problem's parameters are useful and we presented the way to follow to calculate them. The derivatives of an eigenvalue with respect to the problem's parameters show the importance of those parameters on the stability of the FSI system. The last equation is useful for understanding, analysis and design of the aeroelastic system.

Fluid-Structure Interaction 215

Dowell, E.H. (1975). *Aeroelasticity of plates and shells*, Noordhoff International Publishing,

Dowell, E.H., Ilgamov, M. (1988). *Studies in nonlinear aeroelasticity*, Springer-Verlag, New-York Facchinetti, M.L., de Langre, E., Biolley, F. (2004). Coupling of structure and wake oscillators

Ferziger, J.H., Peric, M. (1999), *Computational Methods for Fluid Dynamics* (2nd edition),

Frandsen, J.B. (2004). Numerical bridge deck studies using finite elements. Part I: flutter.

Frink, N.T. (1992). Upwind Scheme for Solving the Euler Equations on Unstructured

Fung, Y.C. (1956). *An introduction to the theory of aeroelasticity*, John Wiley and Sons, New-York Galvao, R., Lee, E., Farrell, D., Hover, F., Triantafyllou, M., Kitney, N., Beynet, P. (2008).

Gnesin, V.I., Kolodyazhnaya, L.V., Rzadkowski, R. (2004). A numerical modelling of stator–

Hoskins, M.H., Kunz, R.F., Bistline, J.E., Dong, C. (2009). Coupled flow–structure–biochemistry

method. *Journal of Fluids and Structures*, Vol. 25, pp. 936–953, ISSN 0889-9746 Howell, R.M., Lucey, A.D., Carpenter, P.W., Pitman, M.W. (2009). Interaction between a

Karmakar, D., Bhattacharjee, J., Sahoo, T. (2009). Wave interaction with multiple articulated

Kimber, M., Lonergan, R., Garimella, S.V. (2009). Experimental study of aerodynamic

Kuo, C.-H., Chen, C.-C. (2009). Passive control of wake flow by two small control cylinders

Langthjem, M.A., Olhoff, N. (2004). A numerical study of flow-induced noise in a two-

Langthjem, M.A., Olhoff, N. (2004). A numerical study of flow-induced noise in a two-

Lesoinne, M., Farhat, C. (1996). Geometric conservation laws for flow problems with moving

Levasseur, V., Sagaut, P., Mallet, M., Chalot, F. (2008). Unstructured Large Eddy Simulation

*Journal of Fluids and Structures,* Vol. 19, pp. 171–191, ISSN 0889-9746

Tetrahedral Meshes, *AIAA Journal*, Vol. 30, No. 1, pp. 70-77

*Structures,* Vol. 19, pp. 1141–1153, ISSN 0889-9746

in vortex-induced vibrations. *Journal of Fluids and Structures,* Vol. 19, pp, 123–140,

Flow control in flow–structure interaction. *Journal of Fluids and Structures,* Vol. 24,

rotor interaction in a turbine stage with oscillating blades. *Journal of Fluids and* 

simulations of dynamic systems of blood cells using an adaptive surface tracking

cantilevered-free flexible plate and ideal flow. *Journal of Fluids and Structures,* Vol.

floating elastic plates. *Journal of Fluids and Structures,* Vol. 25, pp. 1065–1078, ISSN

damping in arrays of vibrating cantilevers. *Journal of Fluids and Structures*, Vol. 25,

at Reynolds number 80. Brief Communication. *Journal of Fluids and Structures*, Vol.

dimensional centrifugal pump. Part I. Hydrodynamics. *Journal of Fluids and* 

dimensional centrifugal pump. Part II. Hydroacoustics. *Journal of Fluids and* 

boundaries and deformable meshes, and their impact on aeroelastic computations.

of the passive control of the flow in a weapon bay. *Journal of Fluids and Structures,*

Leyden

ISSN 0889-9746

Springer, ISBN 3-540-65373-2, Germany

pp. 1216–1226, ISSN 0889-9746

25, pp. 544–566, ISSN 0889-9746

25, pp. 1021–1028, ISSN 0889-9746

*Structures,* Vol. 19, pp. 349–368, ISSN 0889-9746

*Structures,* Vol. 19, pp. 369–386, ISSN 0889-9746

*Comput. Methods Appl. Mech. Engrg.,* 134, pp. 71-90

Vol. 24, pp. 1204–1215, ISSN 0889-9746

0889-9746

pp. 1334–1347

#### **4. Acknowledgment**

This paper was partially supported by the Ministry of Education and Research, Romania, project CNCSIS ID\_919/2007 and project 82-082/2008.

#### **5. References**


This paper was partially supported by the Ministry of Education and Research, Romania,

Aftosmis, M., Gaitonde, D., Tavares, T.S. (1995). Behavior of Linear reconstruction Techniques on Unstructured Meshes, *AIAA Journal*, Vol. 33, No. 11, pp. 2038-2049 Alonso, J.J., Jameson, A. (1994). Fully-Implicit Time-Marching Aeroelastic Solutions, *AIAA* 

Anderson, D. A., Tannehill, J.C., Pletcher, R.H. (1984). *Computational Fluid Mechanics and Heat* 

Appel, J. R., Gunzburger, M. D. (1997). Difficulties in Sensitivity Calculations for Flows with

Batina, J.T. (1991). Unsteady Euler Algorithm with Unstructured Dynamic Mesh for

Beulen, B.W.A.M.M., Rutten, M.C.M., Van de Vosse, F.N. (2009). A time-periodic approach

Bishplinghoff, R.L., Ashley, H., Halfman, R.L. (1955). *Aeroelasticity*, Addison-Wesley

Bishplinghoff, R.L., Ashley, H. (1956). *Principles of Aeroelasticity*, John Wiley and Sons, New-

Cacuci, D. G., Weber, C. F., Oblow, E. M., Marable, J. H. (1980). Sensitivity Theory for General Systems of Nonlinear Equations, *Nucl. Sci. Eng.*, Vol. 75, 88. Cacuci, D. G. (1981). Sensitivity Theory for Nonlinear Systems: I. Nonlinear Functional

Cacuci, D. G. (1981) Sensitivity Theory for Nonlinear Systems: II. Extensions to Additional

Carstens, V., Kemme, R., Schmitt, St. (2003). Coupled simulation of flow-structure interaction in turbomachinery, *Aerospace Science and Technology*, 7, pp. 298-306 Copeland, G.S., Rey, G.J. (2004). Comparison of experiments and reduced-order models for

Dessi, D., Mastroddi, F. (2008). A nonlinear analysis of stability and gust response of

Dessi, D., Mastroddi, F. (2004). Limit-cycle stability reversal via singular perturbation and wing-flap flutter. *Journal of Fluids and Structures,* Vol. 19, pp. 765–783, ISSN 0889-9746

Doaré, O., Michelin, S. (2011). Piezoelectric coupling in energy-harvesting fluttering flexible

turbomachinery high-incidence flutter. *Journal of Fluids and Structures*, Vol. 19, pp.

aeroelastic systems. Brief Communication. *Journal of Fluids and Structures*, Vol. 24,

plates: linear stability analysis and conversion efficiency. *Journal of Fluids and* 

Complex-Aircraft Aerodynamic Analysis, *AIAA Journal*, Vol. 29, No. 3, pp. 327-333 Barrero-Gil, A., Sanz-Andrés, A., Roura, M. ( 2009). Transverse galloping at low Reynolds

numbers. Brief Communication. *Journal of Fluids and Structure*s, Vol. 25, pp. 1236–

for fluid-structure interaction in distensible vessels. *Journal of Fluids and Structures*,

*Transfer*, McGraw-Hill Book Company, ISBN 0-07050328-1, USA

Discontinuities, *AIAA Journal*, Vol. 35, pp. 842-848

*Paper 94-005*6, 32nd AIAA Aerospace Sciences Meeting and Exhibit, January 10-13,

**4. Acknowledgment** 

Reno, NV.

York.

1242, ISSN 0889-9746

713–727, ISSN 0889-9746

pp. 436–445, ISSN 0889-9746

Vol. 25, pp. 954–966, ISSN 0889-9746

Company, Cambridge, Massachussets

Analysis Approach, *J. Math. Phys*., Vol. 22, 2794

Classes of Responses, *J. Math. Phys.*, Vol. 22, 2803

*Structures*, doi: 10.1016/j.jfluidsstructs.2011.04.008

**5. References** 

project CNCSIS ID\_919/2007 and project 82-082/2008.


**10** 

*China* 

**Study on Multi-Phase Flow Field in** 

*National Engineering Research Center for Integrated Utilization of Salt Lake Resources* 

Magnesium, the 8th most abundant element in the earth's crust, was discovered and isolated by Sir Humphrey in 1808. Magnesium is classified as an alkaline earth metal. It is found in Group 3 of the periodic table. It thus has a similar electronic structure to Be, Ca, Sr, Ba and Rd. The density of magnesium at 20°C is 1.738 g/cm3. At the melting point of 650°C this is reduced to 1.65 g/cm3. On melting there is an expansion in volume of 4.2%. Magnesium is the lightest metal in large-scale commercial use. After World War II, the magnesium industry attempted to develop magnesium for a number of applications. Most of its peacetime uses take advantage of the light weight or other chemical properties. The uses of magnesium as a structural material were, however, very few. The bulk was used as an alloying element in aluminium alloys. Other uses, such as deoxidation of steel, chemical and pyrotechnics, outweighed structural uses, especially in energy-saving and environmental protection applications was wide because of its contribution to reduce energy consumption and greenhouse gas emissions. The most successful peacetime application for magnesium was in the original German Volkswagen car that was designed by Ferdinand Porsche. The VW Beetle used large magnesium alloy die castings for the crankcase and the transmission housing (both cast in halves) plus a number of smaller castings. Each Beetle contained more than 20 kg of magnesium alloy. Many of the other applications developed during the World War II, could not be quickly converted to civilian uses. Some of the uses such as aircraft wheels and aircraft engine castings and troop carrying buses were modified and then used for the basis of civilian industries. With more attention to energy and environment, magnesium will hold greater promise as a new weight-saving replacement for denser steel and aluminum alloys, and demands for magnesium will increase sharply in the future. In recent decades, demands for magnesium and its productivity increased sharply shown in Figure 1, the world

There are six sources of raw materials for the production of magnesium: magnesite, dolomite, bischofite, carnallite, serpentine and sea water. These sources differ in the magnesium content, in production methods, and in their origins. Some are mined from mines, some in open

productivity of magnesium reached 860,000 tons at 2007.[1-3]

**1. Introduction** 

**Electrolysis Magnesium Industry** 

Ze Sun, Guimin Lu, Xingfu Song, Shuying Sun,

*East China University of Science and Technology, Shanghai* 

Yuzhu Sun, Jin Wang and Jianguo Yu *State Key Laboratory of Chemical Engineering* 

