**Fluid-Structure Interaction of a Radial Turbine**

Zied Driss, Sarhan Karray, Hedi Kchaou and Mohamed Salah Abid *University of Sfax (US) National Engineering School of Sfax (ENIS) Laboratory of Electro-Mechanic Systems (LASEM) Tunisia* 

### **1. Introduction**

200 Applied Computational Fluid Dynamics

Pezzella, G., Kauffmann, J., Dapra, A., Tomatis, C., *An Italian Contribution to the Next* 

Tomatis, C., Bouaziz, L., Franck, T., Kauffmann, J., *RLV Candidates for European Future* 

Viviani, A., Pezzella, G., *Heat Transfer Analysis for a Winged Reentry Flight Test Bed*,

XX AIDAA Congress, Milan, Italy, June 29-July 3, 2009.

(USA), paper AIAA-2008-2639.

2312.

2006, Valencia, Spain. IAC-06-D2.2.07.

Hypersonic Systems and Technologies Conference. Apr. 28-1, 2008. Dayton, Ohio

*Generation Launcher in the Framework of ESA Future Launchers Preparatory Programme*,

*Launchers Preparatory Programme*. 57th International Astronautical Congress, 2-6 Oct.

International Journal of Engineering (IJE), Volume (3): Issue (3). 2009. ISSN: 1985-

Fluid-structure interaction (FSI) plays a very important role in many industrial applications as diverse as aerospace, automotive, biomedical, civil and nuclear engineering. In the coupling between fluid and structure, the deforming structure modifies the boundary conditions for the fluid due to the motion of the fluid boundary. At the same time, the fluid flow causes varying loading conditions on the structure. Though efficient solvers for both the fluid and the structural dynamics exist, the development of tools for modelling various fluid-structure interaction problems remains a challenge. In order to describe FSI problems, an exchange of data has to take place between the fluid and structural fields. Based on this data exchanged, the methods for solving fluid-structure interaction problems can be divided into monolithic and partitioned couplings. For example, we can mention the researcher of Bobovnika et al. (2005). In this paper, the authors are interested to couple a finite-volume model to a finiteelement model of the straight-tube Coriolis flowmeter. Michler et al. (2004) have developed a monolithic approach to uid-structure interaction. Piperno and Farhat (2001) are interested to partitioned procedures for the transient solution of coupled aeroelastic problems. Sternel et al. (2008) are developed an efficiency and accuracy of fluid-structure interaction simulations using an implicit partitioned approach. Van Brummelen et al. (2003) are studied the energy conservation under incompatibility for uid-structure interaction problems. The monolithic coupling consists of solving both parts of fluid and structure in the same system of equations. The choice of time step is only limited by the required precision. However, this type of algorithm requires a development of a computer code. Moreover, the numerical methods employed for the fluid and the structure domains are different and are difficult to place within the same code, hence the utility of the partitioned algorithm methodology. In this method, the fluid and structural parts are resolved separately. This method has been introduced by Park and Felippa (1983), and further investigated by Wood (1990). Problems dealing with FSI involve the coupling of Computational Fluid Dynamics (CFD) and Computational Structure Dynamics (CSD) codes. Various numerical simulations study the uid behaviour, when the movement of

Fluid-Structure Interaction of a Radial Turbine 203

the post-processor. Data is interchanged between these three modules. The implementation of the computational simulation requires many hypotheses such as the geometrical concept, the choice of the meshing type, the physical models, the boundary conditions and the

The pre-processor serves as a source of input for the CFD or CSD codes by means of an operator interface. This is transformed into a suitable form for the processor. The user activities at the pre-processing stage involve definition of the geometry, mesh generation and specification of the appropriate boundary conditions. The computational model requires that the volume occupied by the fluid inside the vessel is described by a computational grid or cells. In these cells, variables are computed and stored. The computational grid must fit the contours of the vessel and its internals, even if the components are geometrically complex. In our CFD applications, a new meshing method is developed to study the impellers turbines with complex geometries. This method permits the numerical analyses of turbines with simple and complex geometries (Driss et al., 2005, 2007, 2010). In this method, the computer aided design (CAD) is used at first to construct the turbine shape. After that, a list of nodes is defined to belong to the interface separating the

discretization method (Driss et al., 2010).

Fig. 1. Schematic presentation of the FSI method

**2.1 Pre-processor** 

the structure is prescribed analytically (Ralph and Pedley, 1985; Natarajan and Mokhtarzadeh-Dehghan, 2000). Other researchers concentrated on the uid part, while a simple structural model for a rigid body was used (Cossu and Morino, 2000). Further simplications have been done by neglecting the dynamic effects and simulating static FSI. For example, Beckert et al. (2000) used a multivariate interpolation scheme for coupling uid and structural models in three dimensional spaces. They applied it to static aeroelastic problems, in order to predict the equilibrium of elastic wing models in transonic uid ow. The main advantage of the partitioned approach is that it allows already developing efficient and well validated solvers for both the fluid and structure subtasks to be combined. Some theoretical and numerical studies of partitioned coupling algorithms for one and two dimensional problems can be found (Cebral et al., 1995). Commercial codes for combining existing solvers have been developed. Seiber (2002) developed an efficient coupling algorithm combining FASTEST-3D code and FEAP code for solving various uid-structure interaction problems in three-dimensional domains for arbitrary elastic structures. Glück et al. (2003) applied a partitioned coupling between the CFD code (FASTEST-3D) and the CSD code (ASE) to thin shells and membranous structures with large displacements. The latter method have been modified and coupled by MpCCI. Bucchignani et al. (2004) presented a numerical code to study the problem of an incompressible ow in a stirred vessel. It was based on a method of a partition treatment type, with the uid and structural elds resolved by coupling two distinct models. Wang (2008) provided an effective new idea to solve aeroelasticity problems, in which the tools Fluent and ABAQUS/ANSYS are employed.

In this chapter, a computational analysis of the fluid-structure interaction in a stirred vessel equipped with a radial turbine is presented. The hydrodynamic behavior of the turbine was studied previously (Nagata, 1975; Suzukwa et al., 2006). However, they did not take into account the effect of the structural deformations. For this reason, we are interested to study the turbulent flow in a stirred vessel as well as the mechanical deformation of the structure. In a stirred vessel, the fluid can have a significant effect on the deformation of the blades during mixing depending on the turbine speed and the fluid flow. So, a coupling algorithm has been developed to predict the equilibrium of the elastic blade of the radial turbine in a turbulent flow.

### **2. FSI method**

The FSI problems involve the coupling between the CFD and the CSD codes. This coupling is used in the partitioned method (Bobovnika et al., 2005 ; Driss et al., 2007 ; Sternel et al., 2008) when the fluid and structural parts are resolved separately. The coupling interface allows the exchange of the pressures at interface, when the mesh in each fluid and structure field is different and is difficult to place within the same code, hence the utility of the partitioned algorithm methodology. The coupling interface organizes the transfer of information and makes an effective coupling between the fluid and the structure fields. At the fluid-structure interface, the results from the CFD code are applied as a load to the CSD code. The displacements are calculated using the CSD code to characterize the turbine blade deformation. To predict numerically the FSI characteristics of the radial turbine, we have developed a numerical method to simulate the flow depending of the structure displacements. This method involves three modules: the pre-processor, the processor and

the structure is prescribed analytically (Ralph and Pedley, 1985; Natarajan and Mokhtarzadeh-Dehghan, 2000). Other researchers concentrated on the uid part, while a simple structural model for a rigid body was used (Cossu and Morino, 2000). Further simplications have been done by neglecting the dynamic effects and simulating static FSI. For example, Beckert et al. (2000) used a multivariate interpolation scheme for coupling uid and structural models in three dimensional spaces. They applied it to static aeroelastic problems, in order to predict the equilibrium of elastic wing models in transonic uid ow. The main advantage of the partitioned approach is that it allows already developing efficient and well validated solvers for both the fluid and structure subtasks to be combined. Some theoretical and numerical studies of partitioned coupling algorithms for one and two dimensional problems can be found (Cebral et al., 1995). Commercial codes for combining existing solvers have been developed. Seiber (2002) developed an efficient coupling algorithm combining FASTEST-3D code and FEAP code for solving various uid-structure interaction problems in three-dimensional domains for arbitrary elastic structures. Glück et al. (2003) applied a partitioned coupling between the CFD code (FASTEST-3D) and the CSD code (ASE) to thin shells and membranous structures with large displacements. The latter method have been modified and coupled by MpCCI. Bucchignani et al. (2004) presented a numerical code to study the problem of an incompressible ow in a stirred vessel. It was based on a method of a partition treatment type, with the uid and structural elds resolved by coupling two distinct models. Wang (2008) provided an effective new idea to solve aeroelasticity problems, in

In this chapter, a computational analysis of the fluid-structure interaction in a stirred vessel equipped with a radial turbine is presented. The hydrodynamic behavior of the turbine was studied previously (Nagata, 1975; Suzukwa et al., 2006). However, they did not take into account the effect of the structural deformations. For this reason, we are interested to study the turbulent flow in a stirred vessel as well as the mechanical deformation of the structure. In a stirred vessel, the fluid can have a significant effect on the deformation of the blades during mixing depending on the turbine speed and the fluid flow. So, a coupling algorithm has been developed to predict the equilibrium of the elastic blade of the radial turbine in a

The FSI problems involve the coupling between the CFD and the CSD codes. This coupling is used in the partitioned method (Bobovnika et al., 2005 ; Driss et al., 2007 ; Sternel et al., 2008) when the fluid and structural parts are resolved separately. The coupling interface allows the exchange of the pressures at interface, when the mesh in each fluid and structure field is different and is difficult to place within the same code, hence the utility of the partitioned algorithm methodology. The coupling interface organizes the transfer of information and makes an effective coupling between the fluid and the structure fields. At the fluid-structure interface, the results from the CFD code are applied as a load to the CSD code. The displacements are calculated using the CSD code to characterize the turbine blade deformation. To predict numerically the FSI characteristics of the radial turbine, we have developed a numerical method to simulate the flow depending of the structure displacements. This method involves three modules: the pre-processor, the processor and

which the tools Fluent and ABAQUS/ANSYS are employed.

turbulent flow.

**2. FSI method** 

the post-processor. Data is interchanged between these three modules. The implementation of the computational simulation requires many hypotheses such as the geometrical concept, the choice of the meshing type, the physical models, the boundary conditions and the discretization method (Driss et al., 2010).

Fig. 1. Schematic presentation of the FSI method

#### **2.1 Pre-processor**

The pre-processor serves as a source of input for the CFD or CSD codes by means of an operator interface. This is transformed into a suitable form for the processor. The user activities at the pre-processing stage involve definition of the geometry, mesh generation and specification of the appropriate boundary conditions. The computational model requires that the volume occupied by the fluid inside the vessel is described by a computational grid or cells. In these cells, variables are computed and stored. The computational grid must fit the contours of the vessel and its internals, even if the components are geometrically complex. In our CFD applications, a new meshing method is developed to study the impellers turbines with complex geometries. This method permits the numerical analyses of turbines with simple and complex geometries (Driss et al., 2005, 2007, 2010). In this method, the computer aided design (CAD) is used at first to construct the turbine shape. After that, a list of nodes is defined to belong to the interface separating the

Fluid-Structure Interaction of a Radial Turbine 205

is used to calculate the forces applied to the structure. The displacements of the structure are

The CFD code resolves the incompressible Navier-Stokes equations in conjunction with the

<sup>0</sup> *<sup>f</sup> div V* (1)

<sup>2</sup> ( ) div pI div <sup>2</sup>

*<sup>k</sup> d d div V k grad k <sup>G</sup> tD D*

*ff f R f f*

2 2 2 1 2 1 <sup>1</sup>

2 2

2 22 22 2

Re Re

2 2 Re

*k*

2 1 2 1 <sup>1</sup>

*d d div V grad C G C*

*U V W V U WV UW U V <sup>G</sup> r r r z rrr r z z r*

2 *<sup>t</sup>*

The k-ε model is based on the concept of turbulent viscosity where k and ε are calculated by solving the equation above. The value of the turbulent viscosity is deduced from the

> *D k <sup>C</sup> d*

Numerical methods consist on an approximation of the unknown flow variables by means of simple functions and a discretisation by substitution of the approximation into the governing flow equations and subsequent mathematical manipulations (Figure 3). The finite volume method consists on the formal integration of the governing equations of fluid flow over all the control volumes of the solution domain. The transport equations are integrated over its own control volume using the hybrid scheme discretization method. The pressure-velocity coupling is handled by the SIMPLE algorithm of Patankar (1980). The algebraic equations solutions are obtained in reference to the fundamental paper published by Douglas and Gunn (1964). The control volume integration, distinguishes the finite volume method from all other computers techniques. The resulting statements express the conservation of relevant properties for each finite size

*V V Nr N V*

(2)

Re Re

(4)

 

(3)

 

(6)

     

(5)

1 2

standard k-ε turbulence model. The continuity equation is written as follows:

The standard k-ε turbulence equations are given in the following form:

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

The turbulent kinetic energy production takes the following form:

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

*t D k D*

used to update the control volume and the boundary conditions.

The momentum equation is given in the following form:

t *f f*

2 *<sup>t</sup>*

following equation:

cell.

*V*

**2.2.1 CFD code** 

solid domain from the fluid domain passing through the computational structure dynamics (CSD) code. Using this list, the mesh in the flow domain is automatically generated for the three-dimensional simulations. After that, a staggered mesh is used in such a way that four different control volumes are defined for a given node point, one for each of the three vector components and one for the scalar variable. Therefore, the region to be modeled is subdivided into a number of control volumes defined on a cylindrical coordinates system. As application to this method, we are interested to study a stirred tank equipped by a radial turbine (Fig. 2). The tank is a vertical cylindrical vessel with a height-to-diameter ratio H/D of 1. The shaft is placed concentrically with a diameter ratio s/D of 0.04. The tip to tip impeller diameter ratio d/D is a 0.5. The turbine has four blades characterised by the length equal to d/2 and the height equal to h. The turbine is placed in a mid-height tank (z=0.5 H). The geometry system resembles that already studied by Suzukawa et al. (2006) and Nagata (1975). For the numerical computations, the Young's modulus and the Poisson's ratio are equal to E=2.1 1011 MPa and ν=0.3 respectively. In the following investigations, the fluid is assumed to be incompressible and Newtonian. The structure is assumed to be isotropic linear and elastic material law is applied.

Fig. 2. Stirred tank equipped with a radial turbine

#### **2.2 Processor**

In the processor, the coupling interface is used to allow the transfer of the pressure and the displacement between the CFD and the CSD codes. The pressure provided by the CFD code is used to calculate the forces applied to the structure. The displacements of the structure are used to update the control volume and the boundary conditions.

#### **2.2.1 CFD code**

204 Applied Computational Fluid Dynamics

solid domain from the fluid domain passing through the computational structure dynamics (CSD) code. Using this list, the mesh in the flow domain is automatically generated for the three-dimensional simulations. After that, a staggered mesh is used in such a way that four different control volumes are defined for a given node point, one for each of the three vector components and one for the scalar variable. Therefore, the region to be modeled is subdivided into a number of control volumes defined on a cylindrical coordinates system. As application to this method, we are interested to study a stirred tank equipped by a radial turbine (Fig. 2). The tank is a vertical cylindrical vessel with a height-to-diameter ratio H/D of 1. The shaft is placed concentrically with a diameter ratio s/D of 0.04. The tip to tip impeller diameter ratio d/D is a 0.5. The turbine has four blades characterised by the length equal to d/2 and the height equal to h. The turbine is placed in a mid-height tank (z=0.5 H). The geometry system resembles that already studied by Suzukawa et al. (2006) and Nagata (1975). For the numerical computations, the Young's modulus and the Poisson's ratio are equal to E=2.1 1011 MPa and ν=0.3 respectively. In the following investigations, the fluid is assumed to be incompressible and Newtonian. The structure is assumed to be isotropic

linear and elastic material law is applied.

H

Fig. 2. Stirred tank equipped with a radial turbine

In the processor, the coupling interface is used to allow the transfer of the pressure and the displacement between the CFD and the CSD codes. The pressure provided by the CFD code

D

d

s

z

h

**2.2 Processor** 

The CFD code resolves the incompressible Navier-Stokes equations in conjunction with the standard k-ε turbulence model. The continuity equation is written as follows:

$$\operatorname{div}\vec{V}\_f = 0\tag{1}$$

The momentum equation is given in the following form:

$$\frac{\partial(\rho\_f \vec{V}\_f)}{\partial \mathbf{t}} + \text{div}\left(\rho\_f \vec{V}\_f \otimes \vec{V}\_f + \mathbf{p}\vec{\mathbf{I}}\right) = \text{div}\left(\overline{\mathbf{r}} + \overline{\mathbf{r}}\_R\right) + \rho\_f \left(N^2 \,\vec{r} - 2\overline{N} \wedge \vec{V}\_f\right) \tag{2}$$

The standard k-ε turbulence equations are given in the following form:

$$\frac{\partial k}{\partial t} + \operatorname{div} \left[ \vec{V}\_f \, k - \frac{2}{\pi} \left( \frac{d}{D} \right)^2 \frac{1}{\operatorname{Re}} \left( 1 + \frac{\nu\_t}{\sigma\_k} \right) \overline{grad} \, k \right] = \frac{2}{\pi} \left( \frac{d}{D} \right)^2 \frac{1}{\operatorname{Re}} G - \varepsilon \tag{3}$$

$$\frac{\partial \mathcal{E}}{\partial t} + \operatorname{div} \left[ \vec{V}\_f \text{ \( \vec{r} - \frac{2}{\pi} \left( \frac{d}{D} \right)^2 \frac{1}{\operatorname{Re}} \left( 1 + \frac{\nu\_t}{\sigma\_\varepsilon} \right) \overline{g \operatorname{rad}}^\circ \text{ \( \mathcal{E} \)} \right] = \frac{\mathcal{E}}{k} \left[ \mathcal{C}\_{1\varepsilon} \frac{2}{\pi} \left( \frac{d}{D} \right)^2 \frac{1}{\operatorname{Re}} \mathcal{G} - \mathcal{C}\_{2\varepsilon} \mathcal{E} \right] \tag{4}$$

The turbulent kinetic energy production takes the following form:

$$\mathbf{G} = \boldsymbol{\mu}\_t \left[ 2 \left[ \left( \frac{\partial \boldsymbol{U}}{\partial r} \right)^2 + \left( \frac{\partial \boldsymbol{V}}{r \partial \theta} + \frac{\boldsymbol{U}}{r} \right)^2 + \left( \frac{\partial \boldsymbol{W}}{\partial z} \right)^2 \right] + \left[ \frac{\partial \boldsymbol{V}}{\partial r} - \frac{\boldsymbol{V}}{r} + \frac{\partial \boldsymbol{U}}{r \partial \theta} \right]^2 + \left[ \frac{\partial \boldsymbol{W}}{r \partial \theta} + \frac{\partial \boldsymbol{V}}{\partial z} \right]^2 + \left[ \frac{\partial \boldsymbol{U}}{\partial z} + \frac{\partial \boldsymbol{W}}{\partial r} \right]^2 \right] \tag{5}$$

The k-ε model is based on the concept of turbulent viscosity where k and ε are calculated by solving the equation above. The value of the turbulent viscosity is deduced from the following equation:

$$
\omega\_t = \mathbf{C}\_{\mu} \frac{\pi}{2} \left(\frac{D}{d}\right)^2 \text{Re } \frac{k^2}{\varepsilon} \tag{6}
$$

Numerical methods consist on an approximation of the unknown flow variables by means of simple functions and a discretisation by substitution of the approximation into the governing flow equations and subsequent mathematical manipulations (Figure 3). The finite volume method consists on the formal integration of the governing equations of fluid flow over all the control volumes of the solution domain. The transport equations are integrated over its own control volume using the hybrid scheme discretization method. The pressure-velocity coupling is handled by the SIMPLE algorithm of Patankar (1980). The algebraic equations solutions are obtained in reference to the fundamental paper published by Douglas and Gunn (1964). The control volume integration, distinguishes the finite volume method from all other computers techniques. The resulting statements express the conservation of relevant properties for each finite size cell.

Fluid-Structure Interaction of a Radial Turbine 207

To achieve this goal, the CSD code must read the input file already created. After that, the implementation edition can be preceded. At this stage, the new parameters of our application are introduced. Particularly, the analysis type, the mesh shape, the boundary conditions and the result field of the problem are redefined. After the program execution, different results can

In this chapter, we have used a partitioned coupling algorithm between the fluid and structure solver. The fluid-structure interface represents the contact mechanics problem between an elastic structure and the fluid flow. This requires the characterization of boundary conditions exchanged and describe the interaction between the fluid and the structure. Two conditions on

*Vn V n <sup>s</sup> <sup>f</sup>*

(9)

the level of the interface are used. These conditions are given in the following form:

be visualized like the nodes displacement, the strain and the Von Mises stress (Fig. 4).

Fig. 4. Schematic presentation of the CSD code

**2.2.3 Coupling algorithm** 

Fig. 3. Schematic presentation of the CFD code

#### **2.2.2 CSD code**

The CSD code allows the static calculation of the structure and solves the partial differential equations of the structure model. In this study, elastic isotropic structures were considered. The continuity equation and the balance of the momentum equation are given in the following form:

$$\frac{D\rho\_s}{Dt} + \operatorname{div} \left(\rho\_s \vec{V}\_s\right) = 0 \tag{7}$$

$$
\rho\_s \frac{D\vec{V}\_s}{Dt} - \overline{div}\,\overline{\sigma} = \vec{f}\_{ex} \tag{8}
$$

The CSD code allows the static calculation of the structure and solves the partial differential equations of the structure model. In this study, elastic isotropic structures were considered. The continuity equation and the balance of the momentum equation are given in the

> 0 *<sup>s</sup> s s <sup>D</sup> div V*

*<sup>s</sup> <sup>s</sup> ex*

*<sup>ρ</sup> div f Dt* 

(7)

(8)

*Dt* 

*DV*

Fig. 3. Schematic presentation of the CFD code

**2.2.2 CSD code** 

following form:

To achieve this goal, the CSD code must read the input file already created. After that, the implementation edition can be preceded. At this stage, the new parameters of our application are introduced. Particularly, the analysis type, the mesh shape, the boundary conditions and the result field of the problem are redefined. After the program execution, different results can be visualized like the nodes displacement, the strain and the Von Mises stress (Fig. 4).

Fig. 4. Schematic presentation of the CSD code

#### **2.2.3 Coupling algorithm**

In this chapter, we have used a partitioned coupling algorithm between the fluid and structure solver. The fluid-structure interface represents the contact mechanics problem between an elastic structure and the fluid flow. This requires the characterization of boundary conditions exchanged and describe the interaction between the fluid and the structure. Two conditions on the level of the interface are used. These conditions are given in the following form:

$$
\vec{V}\_s \vec{n} = \vec{V}\_f \vec{n} \tag{9}
$$

Fluid-Structure Interaction of a Radial Turbine 209

processing techniques and display. These include geometry domain, grid display, vector plots, surface plots and view manipulation. Several results can be gotten such as the velocity field, the viscous dissipation rate and the turbulence characteristics. Also, it's possible to exchange

The hydrodynamic results are obtained after execution of our CFD code. Particulary, we are interested to the velocity fields, the turbulent kinetic energy, the dissipation rate of the turbulent kinetic energy and the turbulent viscosity. In these conditions, the flow simulation

Figures 7 and 8 show the velocity vector plot (U,W) in two r-z planes defined respectively by the angular coordinate equal to θ=35° and θ=57°. These presentation planes are situated respectively in the upstream and in the downstream of the blade plane. The distribution of the velocity field shows the presence of a radial jet on the level of the turbine which changes against the walls of the tank with two axial flows thus forming two recirculation zones on the two sides of the turbine. Far from the region swept by the turbine, we note a deceleration of the flow. Even if the turbine attains a shape deformation due to the flexion of the blade, the velocity field slightly differs from iteration to another. In fact, during the first three iterations the hydrodynamic structure of the flow varies from iteration to another. From the fourth iteration, the field velocity starts to be preserved and the turbine preserves

are defined by the Reynolds number Re=100000 and the Froude number Fr=0.19.

data through the post-processor with other codes or other interfaces [4].

Fig. 6. Schematic presentation of the post-processor

**3. Fluid-structure results** 

**3.1 CFD code results** 

**3.1.1 Flow patterns** 

his deformed shape.

$$\mathbf{With:}$$

$$
\vec{m} = \vec{n}\_f = -\vec{n}\_s\tag{11}
$$

(10)

In the case of no-slip walls, the compatibility condition for partitioned coupling states that the velocity at the fluid-structure interface must be shared by both the fluid and the structure at the interface. The coupling interface ensures the transfer of the pressure and displacement mentioned above, when the mesh in each fluid and structure codes is different. It organizes the transfer of information in time and makes an effective coupling between the fluid and structure codes. Forces equilibrium must also exist at the interface. The force applied by the fluid to the structure and vice versa should be of the same magnitude but in opposing directions. The structural and fluid equations are solved independently. Loads and boundary conditions are exchanged after each converged increment. The quantities exchanged from the CFD and the CSD codes are respectively the pressure and displacement modes. Then, the coupling algorithm is employed in which the CFD and CSD solvers run concurrently and exchange solution quantities at the end of each converged time step. The advantage here is that in order to solve the basic equations, different numerical solution schemes and analysis software can be applied. The coupling algorithm is based on the sequential process (Fig. 5). In each iteration step, the pressure, calculated in CFD code, is used to determine the force exerted by fluid flow. The geometry is then updated after the calculation of the structural displacement by the CSD code. Therefore, the role of the coupling algorithm is to exchange nearly any kind of data between the coupled codes. The iteration steps are stopped when the convergence is reached.

 *n n* 

*nn n <sup>f</sup> <sup>s</sup>*

Fig. 5. Coupling algorithm

#### **2.3 Post-processor**

The post-processor defines the third module presenting the results obtained after the execution of the numerical codes. A huge resource has been devoted to the development post-

*nn n <sup>f</sup> <sup>s</sup>*

In the case of no-slip walls, the compatibility condition for partitioned coupling states that the velocity at the fluid-structure interface must be shared by both the fluid and the structure at the interface. The coupling interface ensures the transfer of the pressure and displacement mentioned above, when the mesh in each fluid and structure codes is different. It organizes the transfer of information in time and makes an effective coupling between the fluid and structure codes. Forces equilibrium must also exist at the interface. The force applied by the fluid to the structure and vice versa should be of the same magnitude but in opposing directions. The structural and fluid equations are solved independently. Loads and boundary conditions are exchanged after each converged increment. The quantities exchanged from the CFD and the CSD codes are respectively the pressure and displacement modes. Then, the coupling algorithm is employed in which the CFD and CSD solvers run concurrently and exchange solution quantities at the end of each converged time step. The advantage here is that in order to solve the basic equations, different numerical solution schemes and analysis software can be applied. The coupling algorithm is based on the sequential process (Fig. 5). In each iteration step, the pressure, calculated in CFD code, is used to determine the force exerted by fluid flow. The geometry is then updated after the calculation of the structural displacement by the CSD code. Therefore, the role of the coupling algorithm is to exchange nearly any kind of data between the coupled codes. The iteration steps are stopped when the convergence is

The post-processor defines the third module presenting the results obtained after the execution of the numerical codes. A huge resource has been devoted to the development post-

(10)

(11)

 *n n* 

With:

reached.

Fig. 5. Coupling algorithm

**2.3 Post-processor** 

processing techniques and display. These include geometry domain, grid display, vector plots, surface plots and view manipulation. Several results can be gotten such as the velocity field, the viscous dissipation rate and the turbulence characteristics. Also, it's possible to exchange data through the post-processor with other codes or other interfaces [4].

Fig. 6. Schematic presentation of the post-processor

### **3. Fluid-structure results**

### **3.1 CFD code results**

The hydrodynamic results are obtained after execution of our CFD code. Particulary, we are interested to the velocity fields, the turbulent kinetic energy, the dissipation rate of the turbulent kinetic energy and the turbulent viscosity. In these conditions, the flow simulation are defined by the Reynolds number Re=100000 and the Froude number Fr=0.19.

### **3.1.1 Flow patterns**

Figures 7 and 8 show the velocity vector plot (U,W) in two r-z planes defined respectively by the angular coordinate equal to θ=35° and θ=57°. These presentation planes are situated respectively in the upstream and in the downstream of the blade plane. The distribution of the velocity field shows the presence of a radial jet on the level of the turbine which changes against the walls of the tank with two axial flows thus forming two recirculation zones on the two sides of the turbine. Far from the region swept by the turbine, we note a deceleration of the flow. Even if the turbine attains a shape deformation due to the flexion of the blade, the velocity field slightly differs from iteration to another. In fact, during the first three iterations the hydrodynamic structure of the flow varies from iteration to another. From the fourth iteration, the field velocity starts to be preserved and the turbine preserves his deformed shape.

Fluid-Structure Interaction of a Radial Turbine 211

(a) Flat blade (b) 1st iteration (c) 2nd iteration

(d) 3rd iteration (e) 4th iteration (f) 5th iteration

Fig. 8. Flow patterns induced in r-z plane (θ=57°)

Fig. 7. Flow patterns induced in r-z plane (θ=33°)

(a) Flat blade (b) 1st iteration (c) 2nd iteration

(d) 3rd iteration (e) 4th iteration (f) 5th iteration

Fig. 7. Flow patterns induced in r-z plane (θ=33°)

Fig. 8. Flow patterns induced in r-z plane (θ=57°)

Fluid-Structure Interaction of a Radial Turbine 213

blade plane. During five iterations of the coupling algorithm, it's noted that the area of the maximum values are located in the wake which develops an upstream of blades. In all fields swept by the turbine blades, the turbulent kinetic energy remains rather high. Out of this field, the turbulent kinetic energy becomes very weak. The comparison of these results

(a) Flat blade (b) 1st iteration (c) 2nd iteration

(d) 3rd iteration (e) 4th iteration (f) 5th iteration

Fig. 10. Distribution of the turbulent kinetic energy k in r-z plane (θ=35°)

confirms that the maximum values are reached during the second iteration.

### **3.1.2 Axial profile of the radial velocity**

Figure 9 presents the axial profiles of the dimensionless radial velocity component U(z). These profiles are defined in the flat-blade downstream. Under the case, the dimensionless radial coordinates are equal to r=0.33 (Figure 9.a), r=0.5 (Figure 9.b) and r=0.75 (Figure 9.c). These figures adequately portray the swirling radial jet character. Within a flat blade, the axial profile shows a parabolic pace. The extremer is defined at the mid-height of the blade. The same fact is observed within a retreated blade in the 5th iteration. But, the radial velocity component value decreases. These results are also confirmed in the research paper results published by Driss et al. (2005). In this case, we superimposed the experimental results founded within a flat blade by Nagata (1975). The good agreement between the experimental results and the numerical results confirms the validity of the analysis method.

Fig. 9. Axial profiles of the radial velocity component

#### **3.1.3 Turbulent kinetic energy**

Figure 10 presents the distribution of the turbulent kinetic energy in r-z plane defined by the angular coordinate equal to θ=35°. The presentation plane is situated in the upstream of the

Figure 9 presents the axial profiles of the dimensionless radial velocity component U(z). These profiles are defined in the flat-blade downstream. Under the case, the dimensionless radial coordinates are equal to r=0.33 (Figure 9.a), r=0.5 (Figure 9.b) and r=0.75 (Figure 9.c). These figures adequately portray the swirling radial jet character. Within a flat blade, the axial profile shows a parabolic pace. The extremer is defined at the mid-height of the blade. The same fact is observed within a retreated blade in the 5th iteration. But, the radial velocity component value decreases. These results are also confirmed in the research paper results published by Driss et al. (2005). In this case, we superimposed the experimental results founded within a flat blade by Nagata (1975). The good agreement between the experimental results and the numerical results confirms the validity of the analysis method.

> Flat blade (Numerical Model) Retreated blade (5th iteration) Flat blade [Nagata, 1975]

> > z

Radial velocity component U

Radial velocity component U

(a) Radial position r=0.33 (b) Radial position r=0.5 (c) Radial position r=0.66

Figure 10 presents the distribution of the turbulent kinetic energy in r-z plane defined by the angular coordinate equal to θ=35°. The presentation plane is situated in the upstream of the

z

**3.1.2 Axial profile of the radial velocity** 

z

Radial velocity component U

**3.1.3 Turbulent kinetic energy** 

Fig. 9. Axial profiles of the radial velocity component

blade plane. During five iterations of the coupling algorithm, it's noted that the area of the maximum values are located in the wake which develops an upstream of blades. In all fields swept by the turbine blades, the turbulent kinetic energy remains rather high. Out of this field, the turbulent kinetic energy becomes very weak. The comparison of these results confirms that the maximum values are reached during the second iteration.

Fig. 10. Distribution of the turbulent kinetic energy k in r-z plane (θ=35°)

Fluid-Structure Interaction of a Radial Turbine 215

Figure 12 presents the evolution of the turbulent viscosity distribution in r-z plane defined by the angular coordinate equal to θ=35° during iterations of the coupling algorithm. In the field swept by the turbine blades, the turbulent viscosity remains rather high. Near the walls and around the axis of the turbine, the turbulent viscosity undergoes a very fast fall. This is due to the deceleration of the flow. The comparison of these results confirms that the

(a) Flat blade (b) 1st iteration (c) 2nd iteration

(d) 3rd iteration (e) 4th iteration (f) 5th iteration

Fig. 12. Distribution of the turbulent viscosity in r-z plane (θ=35°)

maximum value is reached in the second iteration (Fig. 12-c).

**3.1.5 Turbulent viscosity** 

### **3.1.4 Dissipation rate of the turbulent kinetic energy**

Figure 11 presents the evolution of the distribution of dissipation rate in r-z plane defined by the angular coordinate equal to θ=35° during iterations of the coupling algorithm. Globally, we observe a distribution similar to that already obtained with the turbulent kinetic energy. Indeed, it is noted that the area of the maximum values are located in the wake which develops nearly to the blades end. However, the dissipation rate becomes very weak outside the field swept by the blades. The greatest rate is reached in the fourth and the fifth iterations (Fig. 11-d and 11-e).

Fig. 11. Distribution of dissipation rate in r-z plane (θ=35°)

#### **3.1.5 Turbulent viscosity**

214 Applied Computational Fluid Dynamics

Figure 11 presents the evolution of the distribution of dissipation rate in r-z plane defined by the angular coordinate equal to θ=35° during iterations of the coupling algorithm. Globally, we observe a distribution similar to that already obtained with the turbulent kinetic energy. Indeed, it is noted that the area of the maximum values are located in the wake which develops nearly to the blades end. However, the dissipation rate becomes very weak outside the field swept by the blades. The greatest rate is reached in the fourth and the

(a) Flat blade (b) 1st iteration (c) 2nd iteration

(d) 3rd iteration (e) 4th iteration (f) 5th iteration

Fig. 11. Distribution of dissipation rate in r-z plane (θ=35°)

**3.1.4 Dissipation rate of the turbulent kinetic energy** 

fifth iterations (Fig. 11-d and 11-e).

Figure 12 presents the evolution of the turbulent viscosity distribution in r-z plane defined by the angular coordinate equal to θ=35° during iterations of the coupling algorithm. In the field swept by the turbine blades, the turbulent viscosity remains rather high. Near the walls and around the axis of the turbine, the turbulent viscosity undergoes a very fast fall. This is due to the deceleration of the flow. The comparison of these results confirms that the maximum value is reached in the second iteration (Fig. 12-c).

Fig. 12. Distribution of the turbulent viscosity in r-z plane (θ=35°)

Fluid-Structure Interaction of a Radial Turbine 217

Figure 14 presents the evolution of the deformed turbine shape during iterations of the coupling algorithm. These results are obtained after the execution of the CSD code. The comparison between the results from iteration to another, it's noted that the maximum value is reached in the third iteration (Fig. 14-c). In this case, the retreated angle, defined between the tangent plane at the end of the blade and the radial plane passing by the same end and

(a) 1st iteration: ξ=0° (b) 2nd iteration: ξ=11° 12'' (c) 3rd iteration: ξ=16° 40''

(d) 4th iteration: ξ=5° 11'' (e) 5th iteration: ξ=5° 6'' (f) 6th iteration: ξ=4° 37''

Fig. 14. Deformed shape of the radial turbine

the axis of rotation to the turbine (Driss et al., 2007), is equal to ξ=16° 40''.

**3.2.2 Deformed turbine shape** 

### **3.2 CSD code results**

In this part, we are interested to the radial turbine structure in order to determine the displacement field, the deformed turbine shape and the Von Mises Stress. These results are conducted from the CSD code. The discretization of the domain is ensured by a finite elements method.

### **3.2.1 Displacement field**

Figure 13 presents the evolution of the displacement field of the blades during various iterations of our coupling algorithm. Particularly, we are interested to the displacement of the nodes belonging to the blade. This step is fundamental in the coupling algorithm process to can deduce the new meshing domain. This operation is repeated several times until obtaining the equilibrium. Under these conditions, it is supposed that the deformation of the turbine axis due to torsion is negligible.

Fig. 13. Displacement field of the radial turbine

#### **3.2.2 Deformed turbine shape**

216 Applied Computational Fluid Dynamics

In this part, we are interested to the radial turbine structure in order to determine the displacement field, the deformed turbine shape and the Von Mises Stress. These results are conducted from the CSD code. The discretization of the domain is ensured by a finite

Figure 13 presents the evolution of the displacement field of the blades during various iterations of our coupling algorithm. Particularly, we are interested to the displacement of the nodes belonging to the blade. This step is fundamental in the coupling algorithm process to can deduce the new meshing domain. This operation is repeated several times until obtaining the equilibrium. Under these conditions, it is supposed that the deformation of the

(a) 1st iteration (b) 2nd iteration (c) 3rd iteration

(d) 4th iteration (e) 5th iteration (f) 6th iteration

**3.2 CSD code results** 

**3.2.1 Displacement field** 

turbine axis due to torsion is negligible.

Fig. 13. Displacement field of the radial turbine

elements method.

Figure 14 presents the evolution of the deformed turbine shape during iterations of the coupling algorithm. These results are obtained after the execution of the CSD code. The comparison between the results from iteration to another, it's noted that the maximum value is reached in the third iteration (Fig. 14-c). In this case, the retreated angle, defined between the tangent plane at the end of the blade and the radial plane passing by the same end and the axis of rotation to the turbine (Driss et al., 2007), is equal to ξ=16° 40''.

Fig. 14. Deformed shape of the radial turbine

Fluid-Structure Interaction of a Radial Turbine 219

In this chapter, we have coupled a computational fluid dynamics (CFD) code with a computational structural dynamics (CSD) code by using an efficient coupling interface for solving fluid-structure interaction (FSI) problems in the stirred tank equipped by a radial turbine. Specific techniques of rearrangement of grid and treatment of the boundary conditions are used to follow the behavior of the system. This method takes advantage of the parallel process involved within each analysis code. This allows both parts of the fluidstructure interaction problem to be solved in the best possible way: a Finite Volume Method for the fluid dynamics and a finite element method for the structure. The CFD results obtained allow a visualizing of the velocity field, the turbulent kinetic energy, the dissipation rate of the turbulent kinetic energy, the turbulent viscosity. The CSD results permit to obtain the displacement field, the deformed shape and the Von Mises stress. These results proof that the fluid flow can have a significant effect on the deformation of the turbine blades. So, we can conclude that it's fundamental to take consideration of this phenomenon in the design of the industrial process. In the future, we propose to develop

the numerical method to can study instationary FSI problem.

constant in the standard k-ε model

constant in the standard k-ε model

constant in the standard k-ε model

G turbulent kinetic energy production

U radial velocity components, dimensionless V tangential velocity components, dimensionless W axial velocity components, dimensionless

ε dissipation rate of the turbulent kinetic energy

*<sup>t</sup>* turbulent viscosity, dimensionless

**4. Conclusion** 

**5. Nomenclature** 

d turbine diameter, m D tank diameter, m E Young's modulus, MPa Fr Froude number

h blade height, m

p pressure r radial position Re Reynods number s shaft diameter, m

t time

z axial position

ν Poisson's ratio

density

k turbulent kinetic energy N velocity of the turbine, s-1

*C*1

*C*2

*C*

#### **3.2.3 Von Mises stress**

Figure 15 presents the distribution of the Von Mises stress during iterations of the coupling algorithm. These results are obtained after the execution of the CSD code. The comparison of the results from iteration to another, it's noted that the maximum value is reached in the second iteration (Fig. 15-b).

Fig. 15. Distribution of the Von Mises stress

### **4. Conclusion**

218 Applied Computational Fluid Dynamics

Figure 15 presents the distribution of the Von Mises stress during iterations of the coupling algorithm. These results are obtained after the execution of the CSD code. The comparison of the results from iteration to another, it's noted that the maximum value is reached in the

(a) 1st iteration (b) 2nd iteration (c) 3rd iteration

(d) 4th iteration (e) 5th iteration (f) 6th iteration

Fig. 15. Distribution of the Von Mises stress

**3.2.3 Von Mises stress** 

second iteration (Fig. 15-b).

In this chapter, we have coupled a computational fluid dynamics (CFD) code with a computational structural dynamics (CSD) code by using an efficient coupling interface for solving fluid-structure interaction (FSI) problems in the stirred tank equipped by a radial turbine. Specific techniques of rearrangement of grid and treatment of the boundary conditions are used to follow the behavior of the system. This method takes advantage of the parallel process involved within each analysis code. This allows both parts of the fluidstructure interaction problem to be solved in the best possible way: a Finite Volume Method for the fluid dynamics and a finite element method for the structure. The CFD results obtained allow a visualizing of the velocity field, the turbulent kinetic energy, the dissipation rate of the turbulent kinetic energy, the turbulent viscosity. The CSD results permit to obtain the displacement field, the deformed shape and the Von Mises stress. These results proof that the fluid flow can have a significant effect on the deformation of the turbine blades. So, we can conclude that it's fundamental to take consideration of this phenomenon in the design of the industrial process. In the future, we propose to develop the numerical method to can study instationary FSI problem.

### **5. Nomenclature**


Fluid-Structure Interaction of a Radial Turbine 221

Driss, Z.; Karray, S.; Kchaou, H.; Abid, M.S. (2007). Computer Simulations of Fluid-Structure

Driss, Z.; Bouzgarrou, G.; Chtourou, W.; Kchaou, H.; Abid, M.S. (2010).

Douglas, J. & Gunn, J. E. (1964). A general formulation of alternating direction implicit

Glück, M.; Breuer, M.; Durst, F.; Hlfmann A. & Rank, E. (2003). Computation of wind-

Karray, S.; Driss, Z.; Kchaou, H. & Abid M.S. (2011). Hydromechanics characterization of the

Karray, S.; Driss, Z.; Kchaou, H. & Abid M.S. (2011). Numerical simulation of fluid-structure

Michler, C.; Hulshoff, S. J., Van Brummelen E. H. & De Borst R. (2004). A monolithic approach to uid-structure interaction, Computers & Fluids, 33, 839-848

Natarajan S. & Mokhtarzadeh-Dehghan M. R. (2000). Numerical prediction of a (potential)

Patankar, S. V. (1980). Numerical heat transfer and fluid flow. Series in Computational Methods in Mechanics and Thermal Sciences, McGraw Hill, New York. Piperno, S. & Farhat C. (2001). Partitioned procedures for the transient solution of coupled

Ralph, M.E. & Pedley T.J. (1989). Viscous and inviscid ows in a channel with a moving

Sternel, D. C.; Schäfer, M.; Heck M. & Yigit, S. (2008). Efficiency and accuracy of fluid-

Sieber, G. (2002). Numerical Simulation of Fluid-Structure Interaction Using Loose Coupling

Suzukawa, K.; Mochizukib, S. & Osaka, H. (2006). Effect of the attack angle on the roll and

Van Brummelen, E.H.; Hulshoff S.J. & De Borst, R. (2003). Energy conservation under

Review of Mechanical Engineering, 1, 608-617

Computational Fluid Mechanics, (in press), 5, 1-14

Mechanical Science and Technology, (in press), 25, 1-12

indentation, Journal of Fluid Mechanics. 209, 543-566

Engineering, Darmstadt University of Technology

Applied Mechanics and Engineering, 192, 2727-2748

Engineering Science, 61, 2791- 2798

Nagata, S. (1975). Mixing: principles and applications, Halstead press, Japan.

methods, Num. Math. 6, 428-453

and Structures 17 739-765

Fluids, 32, 711-724

Mech., 43, 103-113

3170

29, 236-245

Interaction Generated by a Flat-Blade Paddle in a Vessel Tank, International

Computational studies of the pitched blade turbines design effect on the stirred tank flow characteristics, European Journal of Mechanics B/Fluids,

induced vibrations of flexible shells and membranous structures. Journal of Fluids

turbulent flow generated by anchors impellers, Engineering Applications of

interaction in a stirred vessel equipped with an anchor impeller, Journal of

soft acting peristaltic blood pump. International Journal for Numerical Methods in

aeroelastic problems, Part II: energy transfer analysis and three-dimensional applications, Computer Methods in Applied Mechanics and Engineering. 190, 3147-

structure interaction simulations using an implicit partitioned approach, Comput

Methods, PhD thesis, at the Department of Numerical Methods in Mechanical

trailing vortex structures in an agitated vessel with a paddle impeller, Chemical

incompatibility for uid-structure interaction problems, Computer Methods in


#### Indices


s structure

### **6. References**


Bobovnika, G.; Mole, N.; Kutina, J.; Stokb, B. & Bajsic, I. (2005). Coupled finite-

Bucchignani, E.; Stella F. & Paglia F. (2004). A partition method for the solution of a

Cebral, C.; Piperno, S. & Larrouturou B. (1995). Partitioned procedures for the transient

Cossu C. & Morino L., On the instability of a spring-mounted circular cylinder

Driss, Z.; Kchaou H., Baccar, M. & Abid M.S. (2005). Numerical investigation of

Driss, Z. (2008). Contribution in studies of the turbines in an agitated vessel, PhD thesis,

Driss, Z.; Bouzgarrou, G.; Kchaou, H.; Abid, M.S. (2011). Computer simulation of the

National School of Engineers of Sfax, University of Sfax, Tunisia

volume/finite-element modeling of the straight-tube Coriolis flowmeter, Journal of

coupled liquid-structure interaction problem, Applied Numerical Mathematics,

solution of coupled aeroelastic problems. Part 1: Model problem, theory and twodimensional application, Journal of Computer Methods in Applied Mechanics and

in a viscous flow at low Reynolds numbers, Journal of Fluids and Structures, 14,

internal laminar flow generated by a retreated-blade paddle and a flat-blade paddle in a vessel tank, International Journal of Engineering Simulation,

laminar flow in stirred tanks generated by the proximity impellers of a mono and double screws type with simple and modified profiles, Mechanics & Industries, 12,

 

*n*

*r* 

*V* 

*ex f*

Indices

f fluid s structure

**6. References** 

*R* 

*<sup>k</sup>* constant in the standard k-ε model

θ angular coordinate, rad

velocity vector, m.s-1

Reynolds tensor

51, 463-475.

183-196.

6, 10-16

109-121

stress tensor of the structure

Fluids and Structures, 20, 785-800.

Engineering, 124, 79-112.

ξ retreated angle

normal vector

stress tensor

I identity tensor

exterior force

vector position

constant in the standard k-ε model


**11**

*1India 2USA* 

**Industrial Sprays: Experimental** 

*1Dow Chemical International Pvt. Ltd., Pune 2The Dow Chemical Company, Freeport, TX* 

**Characterization and Numerical Modeling** 

Sprays in industrial processes are becoming more common as well as more challenging for design. Spray formation processes are critical to the performance of a number of technologies and applications. The spray applications can be quite varied and include combustion systems (gas turbine engines, internal combustion engines, incinerators, furnaces, rocket engines), agriculture (pesticide and herbicide treatments), paints and coatings (furniture, automobiles), consumer products (cleaners, personal care products), fire suppression systems, spray cooling (materials processing, computer chip cooling), gas treatment (thermal quench, acid neutralization), medicinal/pharmaceutical, reactor feed injections, liquid distribution into chemical units (defoaming, uniform flux over column

The process requirement varies depending upon the application. For example, making small drops are important for internal combustion engines and other applications where rapid mass transfer into the gas phase is critical. Large drops may be desired where entrainment in a counter-current gas flow must be minimized. In another case, creating a narrow drop size distribution, neither too small or too large, is desired for agricultural sprays when trying to

Controlling the spread of drop distribution plays very important role in paints and coatings, consumer products, medicinal delivery, and spray drying applications. On the other hand, distributing liquid mass more uniformly throughout the spray and increasing the fraction of liquid that impacts a target are key requirements in paints and coatings, spray cooling, and

Injector designers face challenges in designing the atomization process for achieving the desired process output. Despite immense progress in laser diagnostics for spray flows, the reliable experimental data on atomization characteristics are still scarce due to measurement difficulties associated with the dense spray zones formed around the liquid jet core. As a result, application of empirical correlations for atomization flow characteristics has a very narrow range of applicability. It is therefore, essential to develop and apply new tools to enhance our understanding of the atomization process. Such an understanding will be

achieve low driftable fines while at the same time maintaining good target coverage.

useful in devising cost effective and reliable scale-up of spraying operations.

packing), and spray drying (foods, drugs, materials processing).

**1. Introduction** 

fire suppression applications.

Avinash Khopkar1, Michael D. Cloeter2 and Quan Yuan2

Wang, Y. (2008). Combination of CFD and CSD packages for fluid-structure interaction, Journal of Hydrodynamics, 20, 756-761

Wood, W.M. (1990). Practical Time-stepping Schemes, Oxford University Press, New York
