A methodology and computational system for the simulation of fluid-structure interaction problem

TECHNICAL PAPERS

A methodology and computational system for the simulation of fluid-structure interaction problem

A. R. E. Antunes^I; P. R. M. Lyra^II; R. B. Willmersdorf^II

^IDepartamento de Engenharia Civil – UFPE; Av. Acadêmico Hélio Ramos, S/N; 50740-530 Recife, PE. Brazil; [e-mail]

^IIDepartamento de Engenharia Mecânica – UFPE; Av. Acadêmico Hélio Ramos, S/N; 50740-530 Recife, PE. Brazil; [e-mail]; [e-mail]

ABSTRACT

In this paper a flexible finite element computational tool developed to investigate fluid-structure interaction applications in two dimensions is described. We consider problems which can be modelled as a viscous incompressible fluid flow and a rigid body-spring system interacting nonlinearly with each other. The coupling is dealt with the use of an interface approach, in which each physics involved is solved with different schemes and the required information is transferred through the interface of both systems. This approach is, at least in principle, very flexible and computationally efficient as the best available scheme can be adopted to solve each physics. Here, a stabilized FEM considering the "ALE" (Arbitrary Lagrangian-Eulerian) formulation with Crank-Nicholson time-integration is employed for the fluid-dynamics analysis, and the Newmark Method is used for the structural dynamics. Several important tools were incorporated into our system including different possibilities for the mesh movement algorithm, the computational domain decomposition into regions with and without mesh deformation, and the remeshing strategy (either global or local) to keep the necessary mesh quality. As application we present a study of the forced lock-in phenomena and a preliminary investigation on the suppression (or at least the reduction) of the vortex induced vibrations (VIV) on a solid circular cylinder using an idealization of a periodic acoustic excitation.

Keywords: Fluid-structure interaction, vortex induced vibrations (VIV), finite element method (FEM), arbitrary Lagrangian-Eulerian (ALE) formulation, lock-in phenomena, suppression of structural vibration

Introduction

Several practical structures, in different engineering fields, are subjected to vibration as a result of flow induced phenomena. Such behavior can compromise the integrity of the structure or make it uncomfortable for human use. The analysis of these problems involve the study of a coupled fluid-structure interaction model and can be done using computational modeling. The numerical simulation of such applications is most commonly performed using an interface approach and involves the modification of the computational domain as the geometry under consideration is moving with time. In order to avoid updating the computational discretization too frequently, an Arbitrary Lagrangian Eulerian (ALE) formulation (Nomura and Hughes, 1992; Childs, 1999) is normally adopted together with a mesh movement algorithm. In such approach, generally, the reference frame at the moving interface between the structure and the fluid has a Lagrangian description, and regions away from the interface have a mixed Lagrangian and Eulerian description to accommodate the arbitrary movement of the frame of reference.

In this paper we briefly describe the finite element procedure developed to simulate the two dimensional fluid-structure interaction of a rigid circular cylinder, supported by elastic springs, immerse in an incompressible viscous fluid flow. The analysis of such model application gives insight on many problems of industrial interest, for instance the study of "VIV" (Vortex Induced Vibrations) (Blevins, 1986) on offshore platform legs. The adopted procedure uses a stabilized Petrov-Galerkin/Generalized Least Squares ALE finite element formulation with Crank-Nicholson time-integration for the fluid-dynamics analysis (Sampaio, 1991 and Sampaio et al., 1993). This scheme represents an SUPG-like algorithm (Streamline Upwind Petrov-Galerkin (Brooks and Hughes, 1982)) with the optimal upwind parameter implicitly determined through the formulation and a timescale analysis (Sampaio et al., 1993). For the structural analysis a simple lumped model with three degrees-of-freedom and the Newmark Method (Hughes, 1987) is used. The fluid-structure coupling is solved through forced coupling of the variables at the interface and implemented in a segregated approach, using an algorithm to control errors due to the existing time delay between the fluid and structural analysis (Blom and Leyland, 1998). Several alternatives for the subdivision of the domain into subdomains, where the different descriptions (Eulerian, Lagrangian or ALE) are adopted, were incorporated in our computational system. Different types of mesh movement and mesh smoothing were implemented in order to reduce the distortion of the computational meshes over the fluid domain. The capabilities of automatically generating and adapting the mesh (using local or global remeshing) (Lyra & Carvalho, 2000) were also incorporated into the computational system, allowing the study of problems with large displacements. All these features will be briefly described. Further details and the results to many applications dealing with different levels of difficulty and interaction between fluid and structure were analyzed to validate the computational system developed and can be seen in Antunes et al. (2002), Antunes (2002) and Lyra et al. (2003). The obtained results are in good agreement with the experimental, theoretical and numerical data available in the literature. The analysis of the lock-in phenomena is described and a preliminary numerical study on the suppression of the vortex-induced vibration of a circular cylinder by acoustic excitation is presented. The results shows the influence of the acoustic excitation frequency on the vortex shedding frequency and on the amplitude of the structural vibration. Finally, we draw the most important conclusions and some on going and future extension of this research.

Numerical Formulation

Most fluid structure interaction problems where there is a strong coupling between the displacement of the structure and the flow field are characterized by large displacements of some of the boundaries of the domain. The regions close to these moving boundaries are more naturally discretized with a Lagrangean approach. The fluid regions away from the moving boundaries, however, are more naturally treated with a conventional Eulerian formulation, with a fixed reference frame. We use an Arbitrary Lagrangean Eulerian framework to combine these two approaches in a single numerical technique. The differential equations that described the dynamics of the fluid and the structure therefore must be written in this framework.

Fluid Dynamics

Standard Eulerian Formulation

The flow of incompressible fluids can be described by a specialization of the general Navier-Stokes Equations, where we will also consider that the viscosity is constant and that the fluid is Newtonian. Unless otherwise stated, in the following we will use indicial notation with the summation convention. Within an Eulerian framework, i.e., using a fixed frame of reference and fixed control volumes, the incompressible Navier-Stokes equations in non-conservative form reduce to:

In the above equation, i,j = 1,2, x_i, are the spatial coordinates, t is the time variable, r is the density of the fluid, u_i are the components of the velocity of the fluid, b_i are the external body forces, and t_ij is the stress tensor. Equation (1) is subjected to the incompressibility restriction

and the stress tensor is given by

where µ is the dynamic viscosity, p is the pressure and d_ij is the Kroenecker's delta. Equations (1), (2) and (3) are written for a fixed geometric domain W_f and for a time interval I. For a well posed problem, it is also necessary to impose boundary conditions on G, the boundary of the domain W_f, and initial conditions on W_f. The boundary conditions are known velocities on G_u and known surface tractions at G_t, with G = G_u È G_t and G_u Ç G_t = 0. The boundary conditions associated with the mass balance are given in terms of known pressure at G_p and known mass flux at G_G, with G = G_p È G_G and G = G_p Ç G_G = 0. Here, G = ru_in_i with n_ibeing the outward unit normal vector to G. The initial conditions are known velocities and pressure on W_f in the initial time of the analysis.

ALE Formulation

To develop a finite element discretization applicable to deformable domains, we use an ALE formulation, as proposed in Hughes et al (1981). We define the three domains shown in Fig. 1: the spatial domain W_Y which is the physical space defined by the material particles at time t; the referential domain W_X, which is a fixed domain whose image at time t, subjected to a transformation , is the spatial domain; and the material domain W_Z, which is the domain occupied at the time t = 0 by the material particles that occupy the spatial domain at time t. If we define the transformation from the material domain to the spatial domain as f, i.e., f:W_Z ® W_Y, then the transformation from the material domain to the referential domain is given by y:W_Z ® W_X, where y = ^-1ºf (º is the functional composition operator). Now let us consider that z_i and x_i represent particles in W_Z and W_X whose image at time t is y_iin W_Y.

By adopting the following notation: S for spatial, R for referential and M for material domain, respectively, and the partial differentiation operators will be indicated by the following notations:

It follows directly from the previous definitions that

Considering the transformation f,

where d_Si is the displacement, n_Sij is the velocity and g_sij is the deformation gradient. In a similar manner, considering the transformation , we can write:

where d_Ri is the displacement,n_Ri is the velocity and g_Rij is the deformation gradient. Finally, considering the transformation y, we can write:

where d_Mi is the displacement, n_Mi is the velocity and g_Mij is the deformation gradient. Equations (6) to (8) are the kinematic relationships for the different descriptions of a continuum. The equations (8) are the classical kinematic relations of the Lagrangian description of a continuum.

With the above mappings, we can write (Hughes et al. 1981) that:

where c_i is defined as the convective velocity. If f is a scalar function of the flow field, the material derivative of this function is given by (Childs, 1999):

where is the deformation gradient between the reference and the spatial domain, and c_i is the convective velocity defined in Eq. (9). In this work, we used a single step Crank-Nicholson scheme for the time integration, therefore throughout the duration of each and every time step, the referential and spatial domains are coincident. The deformation gradient between these two domains is the identity transformation, so Eq. (10) reduces to:

It is worth noting that if n_R = 0 in Eq. (9) and (11), we have the conventional Eulerian description, and if n_R = n_Swe have the Lagrangean description. With this definition of the material derivative, the incompressible Navier-Stokes equations in the ALE formulation, subject to the divergent free velocity, reduce to:

Finite Element Discretization

Initially, the momentum equations are discretised in time using a q-averaged finite difference scheme. The Crank-Nicolson scheme, q=1/2, was adopted because it gives second order time accuracy. It leads to

where b_i = rf_i denote the body forces, and

with the superscript n+1/2 refering to the arithmetic average of the variable values at time-levels n and n+1, and a suitable linearization for the convective term has been achieved by the use of the approximation .

Suppose that the spatial domain W has been discretized using linear triangular elements. An approximate solution which belongs to the subspace of the trial functions Â^(p) is sought. The subspace of trial functions has dimension p, equal to the number of discrete nodes, and is defined by

where the superscript hat (^) define the approximated function, and with denoting the initial condition in domain W at time t=0, and N_J are trial or shape functions. The introduction of the discretization of U represented by Eq. (15) into expressions (13)-(14) leads to the approximated quantities and a least-squares discrete formulation of the problem given in Eq. (13) can be stated as

Here it should be observed that all terms are defined on the elements interiors, where the shape functions are of class C^¥, thus avoiding C¹ the inter-element continuity requirements of standard least-squares formulations when used to solve second-order partial differential equations. The solution of the equivalent problem described in (16) is achieved by minimizing the functional P with respect to the free parameter , i.e.

This expression can be re-written as

which can also be regarded as a Petrov-Galerkin weighted approximation of the momentum balances, where

where, N₁are the linear shape functions defined on the element.

Note that the second-order term in W_I vanishes and that the resulting weighting function is discontinuous across element boundaries for the adopted C⁰ linear shape function N_I. The parameter Dt/2 acts to bias the integral in Eq. (18) in favour of the upwind term of the trial function, i.e. first term of Eq. (19), and so upwinding is incorporated into the finite element framework in the streamline direction. The resultant weighting functions N_I + W_I have an equivalent structure to those employed in the standard SUPG (Brooks and Hughes, 1982). The time step is chosen according to

where h is the element size, and (Re)^E is the element Reynolds number. For linear triangles h is the smallest triangle height. Sampaio (1991) presented some heuristic arguments, concerning the time-scales of convection-diffusion processes which are representable in a given mesh, to justify the use of the expression (20) for a multidimensional incompressible Navier-Stokes algorithm.

A term which refers to the viscous flux boundary conditions and their compatibility across element interfaces, G_E, must be added to the Petrov-Galerkin method given in Eq. (18) in order to have a well-posed formulation. Here this is enforced, in a weak form, according to

Assuming the boundary conditions, where I refers to time interval

where G_D refers to Dirichlet boundary condition and G_N refers to Neumann boundary, and the compatibility conditions for the exact viscous flux

where G_E represents the internal inter-elements boundaries.

Considering an homogeneous problem, i.e., without body forces, the approximate weak form of the Eqs. (12), for the momentum equations becomes

The continuity equation has not been considered and a discretization following the standard mixed formulation is dismissed, because the same interpolation for all variables has been used and this would imply the violation of Babuska-Brezzi condition. Following Sampaio et al. 1993, the incompressibility condition is used to obtain Eq. (25) when the sum of the squared residuals is minimized with respect to the pressure parameters, and after some manipulations, results in,

which represents an elliptic-type equation. The solution of this equation is sought over the domain W for all t > t⁰, subject to the boundary conditions

where, if , at least one pressure reference value must be prescribed to determine a unique pressure field.

In the Eq. (25) the boundary integral is computed only at boundaries with prescribed velocity, being different from zero just for moving boundaries. This term is responsible for a contribution into the pressure-continuity equation producing a reduction or an increase in the pressure near the moving boundary, which in turn is responsible for the lift, that induces transversal oscillations of the solid body.

The integrals in Eq. (24) and (25) can be evaluated by summing individual contributions and, by adopting a compact matrix notation, can be summarized as

The matrices H₁, H₂, J, D, B^p, B^v represents the various terms of the Eq. (24) and (25).

These equations are solved in a segregated manner, the pressure-continuity equation are first solved to obtain the nodal values of the pressure at t^n+1/2, followed by the solution of the momentum equations, using the updated pressure p^n+1/2, to obtain the velocities u_iⁿ⁺¹. The matrices K^pand K^v are symmetric positive definite and the conjugate gradient method, with Jacobi preconditioner, is employed to solve these algebraic systems of equations. See Sampaio et al. (1993) for further information when using an Eulerian formulation, which has been extended in this work for ALE formulation.

Structural Dynamics

In this work, we only consider dynamics of rigid bodies. The movement of the body is obtained with a straightforward application of the Newmark's method (Hughes, 1987) which results in the following system of equations:

where d, v, and a are the displacement, velocity and acceleration vectors, and b, g are parameters of the method, M, C, and K are the mass, damping and stiffness matrices of the structure. In this work we used g = 1/2 and b = 1/4 , which leads to an implicit, second order accurate and unconditionally stable time integration scheme. The stability of this scheme is important because the time step increment for the structural time evolution is taken as the same as the time increment chosen for the CFD solution. This time increment is determined by the stability requirements of the CFD algorithm, and therefore its time scale is completely unrelated to the dynamic behavior of the structure.

Fluid Structural Coupling

The coupling between the fluids and structural field is imposed in a segregated manner, and the compatibility conditions are imposed a posteriori, at the end of each time step, to enforce the consistency of the interface between both fields. This approach has the advantage, in contrast to a monolithic approach, that the most efficient numerical solution technique can be used for each particular field.

At the start of each time step, we assume that the interfaces are consistent, i.e., that points lying on the interface between the structural and fluid domain have the same velocities, when considered belonging to either domain. Each domain is then advanced in time independently, according to its own physics. The other domain is considered static, and used only as a source of initial and boundary conditions for the current time step. At the end of the time step, the interfaces will therefore no longer be consistent, therefore a series of predictor-corrector steps is repeated until satisfactory agreement between the two fields is reached. Usually incompressible fluid flow analysis requires smaller time steps than the structural analysis in a couple fluid-structural applications. This fact could be exploited by advancing the flow problem several time steps before interchange data within the two sub-problems. However, as our structural model is very simple, and so very cheap computationally, we do not exploit this fact.

As the interfaces between domains move along the time, the mesh is distorted by the movement of the boundaries of the domain, and when this distortion becomes excessive, the inadequate elements are removed and the mesh is recreated in these regions. The predictor-corrector technique adopted is adapted from the one proposed by Bloom and Leyland (1998), and the general procedure is summarized below:

a. For all time step do:

b. Estimate a predictor velocity v _p

c. Move the structure and the mesh in the computational domain

d. If the quality of the mesh is not satisfactory, then:

e. Delete distorted elements and recreate the mesh

f. Interpolate the solution to the new mesh

g. Solve the CFD problem

h. Solve the structural dynamics problem

i. Compute a corrector velocity v _c

j. If v _p and v _c converged to each other, then:

k. Advance to next time step

l. Else:

m. v_p = v_c

n. Repeat from step c

The predictor velocity v_p at the start of each time increment (step b) is taken as an estimate of the current velocity of the structure. This velocity is computed by a simple linear extrapolation from the structural velocity and acceleration computed at the previous time step. As the interpolation of the solutions between meshes (step f) is not cheap and can introduce errors, we try to minimize it by restricting the regions where the mesh is allowed to move. This will be further described later. The CFD solution (step g) uses as boundary conditions for the moving boundaries the current velocity of the interface, v_p. The corrector velocity v_c (step i) is the interface velocity that results from the solution of the structural dynamics problem (step h). The convergence test (step j) verifies if the difference between the corrector velocity v_c and v_p are less than some prescribed tolerance. If this is so, the algorithm advances to the next time step (step k), otherwise the current corrector velocity is taken as the new predictor velocity and the predictor-corrector loop is repeated. In this work, we only dealt with rigid bodies, therefore the velocity of the interface nodes can be easily computed from the velocity of the body with simple transformation matrices (Nomura and Hughes, 1992).

Numerical Tools

Many practical issues must be dealt with to obtain a successful computational implementation of the procedures described above. Clearly, facilities for dealing with deformable domains, which involve automatic mesh generation, assessment of mesh quality, automatic mesh movement and regeneration are all important aspects. The most interesting aspects of our implementation will be described below.

Subdomain Partitioning

We adopted the technique of partitioning the domain into subdomains in which different descriptions can be employed for the discretization on each subdomain. We have implemented three distinct possibilities: a conventional Eulerian formulation, an ALE formulation with a deformable mesh and an ALE formulation where the mesh is movable but not deformable. The Eulerian formulation is more appropriate for regions of fluid away from the moving boundaries. The ALE formulation with non-deformable but moving meshes is used in the regions very close to the moving interfaces. These meshes are attached to the moving bodies, and are not deformable because, in general, complex flow phenomena that require very fine and high quality meshes, such as boundary layer formation or separation bubbles, occur in these regions. Deformation of the mesh would quickly destroy the quality of these meshes and compromise the quality of the solution. The ALE formulation with deformable meshes is used to couple the two types of domain just described. A simple sketch example can be seen in Fig. 2. For the structure, of course, we use a Lagrangian description.

The use of different formulations for different subdomains is also important because it minimizes the needed interpolation of solutions between different meshes. Clearly, in the subdomains where an Eulerian formulation is used, the interpolation is completely unnecessary because the mesh never varies. Interpolation between unstructured meshes, even with the use of adequate data structures to speed up searches, is a somewhat costly procedure. The interpolation can also reduce the accuracy of the numerical solution (Sampaio et al. 1993), therefore it is important that it remains restricted to the smallest possible regions. In the current implementation, the subdomain decomposition and the assignment of formulations is done a priori by the user, therefore some knowledge of the expected amplitudes of displacement is necessary. In practice, however, this has not proven itself to be a problem. The alternate digital tree (ADT) (Bonet and Peraire, 1990) is used to find out within which element, on the mesh prior to its movement or to a remeshing, each node of the "new" mesh falls inside.

The mesh generation inside each subdomain is fully automatic, using an advancing front mesh generator (Lyra & Carvalho, 2000). The mesh density is controlled by a background mesh, and of course the generator enforces compatibility between the boundaries of the subdomains. This generator has been modified to produce highly elongated elements along solid surfaces using an advancing layer technique.

Mesh Movement

In the domains with a deformable ALE formulation, the movement of the interface nodes causes distortion of the elements connected to these nodes. Two techniques are used to minimize these effect: mesh movement and remeshing. Initially, a procedure akin to a mesh smoothing is used. The mesh is viewed as a network of elastic springs, see the Fig. 3, where each edge of the mesh corresponds to a spring whose stiffness is inversely proportional to the edge length. Then an elastic problem is solved, where the boundary conditions are given by the new positions of the nodes on the interfaces of the subdomains. A simple direct iteration is done, and typically, very few iterations are sufficient for convergence to the required tolerance. The location of the mesh's inner nodes are then determined by solving the static equilibrium equations which result from summing the forces in the spring system at each node. Hence the algorithm changes the positions of the interior nodes without altering the topology of the mesh. In this formulation the summation extends over all the nodes which surround node I and sufficient convergence is normally achieved in 10 iterations. By starting the iteration with the previous displacement at each node, the motion of the boundary is allowed to spread throughout the mesh. This procedure is therefore very fast and will be detailed in what follows.

For each edge connecting the nodes I and J, the stiffness is given by

A linear extrapolation is applied to pre-define the displacements for the node I, e , in x and y directions, as shown below:

The resulting system equations can be expressed in the following compact form:

where,

The new position of the node is obtained as:

This process is suitable for small displacements. When considering large displacements torsional springs are also necessary to avoid large deformation of the mesh (Farhat et al. 1997).

Another possibility is to prescribe the mesh velocity to be an analytic function based on the distance (r), of the domain point from the moving surface (Nomura, 1994, Azevedo, 1999 and Löhner, 2001), as sketched in the Fig. 4.a. This function f(r) assumes the value of unity for r = 0 and decays to zero as (r) increases. A linear function, as shown in Fig 4.b, is assumed, where given the distance r and the velocity point on the surface closest to domain point w(xG_C), the mesh velocity is given by

This procedure is extremely fast and efficient if the distance is obtained easily. However, if several moving bodies are present in the flow field, this procedure can be complex and insufficient, especially considering 3-D simulations. There are others alternatives, such as smoothing the velocity field, by solving a Laplace equation (Löhner, 2001). The quality of the elements inside the deformable ALE subdomains are verified at each time-step iteration. Such verification will define the need of reconstructing the mesh, i.e., to realize a remeshing.

Remeshing

There are cases, however, where the distortion of the elements is too severe, and the smoothing procedure breaks down. There are particular configurations, for instance, that are prone to element collapse. When this happens, either the whole mesh is re-generated (global remeshing) or the compromised elements are removed from the mesh and the mesh is re-generated in the remaining "holes" (local remeshing). The same mesh generation algorithm that was used to create the original mesh is used inside each "hole" (Lyra and Carvalho, 2000). The remeshing procedure, either local or global, is very attractive as the size of the mesh does not tend to grow unnecessarily. The local remeshing is even more attractive when just local effects are present in the problem, as the mesh is just re-built locally reducing the cpu time used by the mesh generation and inter-mesh interpolation steps. Furthermore, it reduces the error associated with the interpolation process.

The quality parameters of the computational elements refer to the change in area and the internal angles of the elements. An element is considered extremely deformed when:

1. The element has a area change which is greater than a predefined tolerance "A";

2. The element has an internal angle (q ) between q_ref< q < 180°– q_ref.

Once all distorted elements are determined by the previous procedure, a remeshing algorithm is applied if:

1. At least one element have an internal angle so that q < q_min or q > q_max;

2. More than prescribed percentage B of the elements are distorted.

In this work, A = 20%, q_ref = 30°, q_min = 2°, q_max = 170° and B = 2%.

The implementation of the "advancing front method" (Peraire et al., 1999) adopted in our mesh generator requires the knowledge of the mesh parameters (mesh spacing and stretching factor) desired for the mesh to be built, which is provided through the use of a background mesh that covers the whole domain to be discretised. Here, the background mesh is the mesh obtained after the displacement of the structure with the corresponding nodal movement procedure.

Since we did not incorporate an error analysis in our system yet, which would provide the "optimal" mesh for each stage of the analysis, some alternatives to compute the mesh parameters have been devised (Antunes et al., 2002). However, only isotropic meshes were adopted in the present study, so only this particular situation will be described here, as the other alternatives have been devised to deal with anisotropic meshes. Let's consider a node "p" in the mesh and the following definitions:

d_pe - height of the triangle "e" corresponding to node "p";

n_el - number of triangles surrounding node "p".

We loop through all the nodes of the background mesh, and

1. Identify all the elements around node 'p';

2. Compute for each element d_pe and for each node the average value

This strategy allows less mesh reconstructions, being computationally more efficient. An example of this procedure is shown in Fig. (6), where the circular cylinder suffers translation and rotation.

Applications

The scheme described above has already been utilized for the simulation of a variety of applications (Antunes et al., 2002, Antunes and Lyra, 2002 and Antunes, 2002), involving different level of interaction between an external fluid flow and a rigid circular cylinder, including the study of the flow around a fixed cylinder, the flow around a cylinder with an imposed periodic displacement, the free vibration of the cylinder in a stationary fluid and the coupled fluid-structure problem of a rigid cylinder supported by elastic springs free to interact with the surrounding fluid flow.

Vortex "Lock-in" Regime Under Forced Oscillating Cylinder

The numerical procedures described is used to simulate the external flow over a moving cylinder with a prescribed forced cross-flow oscillation, for several different Reynolds number. A detailed description of the computational domain with boundary and initial conditions can be seen in the Fig. 7.

The following boundary conditions are applied in this example: non-slip condition on the surface of the rigid body; uniform velocity field with u₁ = u⁰ e u₂ = 0 is prescribed in the face AB; for the faces AC and BD the boundary condition are u₂ = 0 e t1 = 0; on the boundary CD, the pressure value is prescribed, with p = 0 and traction free condition, e.g., t1 = t2 = 0. The initial condition is the velocity field with the components u₁ = u⁰ = u₂ = 0, which are specified all over the domain in the initial time t = t⁰, and a reference pressure p = 0. In this notation, the traction boundary condition is represented by " t_i", where the subscript i, represents the direction of application, and the time is represented by " tⁿ", where the superscript n represents the level time. The vertex of the domain are A(-5,-5), B(-5,5), C(10,-5) e D(10,5), and the circular cylinder have a unitary diameter and is centered on the origin of the orthogonal coordinates.

By computing the fluid flow considering a fixed cylinder we obtain the variation of the frequency of the vortex shedding at different Reynolds numbers. The computed frequencies of the vortex shedding, presented in the first column of Table 1, agree with the experimental data (Blevins, 1986).

The "lock-in" phenomenon (Blevins, 1986) occurs when the response frequency of the cylinder synchronises with the vortex shedding frequency. By imposing an harmonic oscillation to the cylinder and by analysing the modification on the vortex shedding frequency such phenomena can be studied. Taking the vortex shedding frequency obtained with the analysis of a fixed cylinder at Re=120 (0,1795 Hz) as the forced frequency imposed to the cylinder, varying the Reynolds number from 100 to 160 and the amplitude ratio (Ay/D) from 0% to 8% we built the Table 1 and Figure 8. These show the influence of the amplitude of the forced displacement on the vortex shedding frequency. The flow regime are "locked-in", i.e. the displacement and vortex shedding are in phase, for all values of Reynolds number and displacement amplitude ratio presented in the table with bold numbers. Whenever the vortex shedding frequency is different from the one obtained with a fixed cylinder we say that it is in a "perturbed" regime, and when the imposed movement does not change the vortex shedding frequency, the regime is unperturbed and out of "lock-in". Figure 8 is a graphic representation of Table 1 showing the discussed behaviour. Our results are in good agreement with that presented in the literature (Blevins, 1986; Correia, 2001). By increasing the displacement ratio to 2% we get a perturbed regime, and then a "lock-in" regime for all Reynolds number analysed, for 3,5% and 5%. If we increase the amplitude ratio to 8% the flow gets out of the "lock-in" regime for most Reynolds number studied except for Re=115.

The harmonic function to displace the cylinder adopted is

where, a is the amplitude displacement, f is the frequency, and t⁰ is the initial time. In this work t⁰ = 0, n = 1, f = 0,1795 Hz and a varying.

Time histories of the drag and lift coefficients are plotted in Figure 9(a) for Reynolds number of 160 and forced transverse oscillation with displacement amplitude ratio of 5%. For such Reynolds number and amplitude of oscillation the flow is initially perturbed by the oscillation of the cylinder and latter the "lock-in" phenomena is established. This can be further observed in Figure 9(b), which shows the transversal velocity (Uy) and the relative displacement amplitude, being in phase after time around 40.

A contour plot of the computed pressure can be seen in Figure 10, where the Von Karman vortex street behind the cylinder is characterized, for the flow at Re=160 and Ay/D=5%.

Suppressing Vortex Induced Vibration

In the literature (Hiejima et al., 1997 and references therein) several experimental and numerical results are reported in which acoustic excitation is applied to an external flow to increase the momentum transfer from the outside flow to the boundary layer and eliminating (or delaying) separation and suppressing (or reducing) vortex induced vibrations in different solid configurations. In this article, our computational system is used to perform an initial study on the behavior of the fluid-structure problem described by Hiejima et al. (1997), in which an idealization of the acoustic excitation is obtained through the application of a periodic velocity excitation on two points at the cylinder surface (see Fig. (11)). The angle between the stagnation point and the excitation point is f_a = 80º. The excitation velocity is given by

where, U_a and f_a refer to the periodic velocity excitation amplitude and frequency, respectively. The amplitude U_a is chosen to be equal to 10% of the free-stream velocity U, and the two excitation velocities are in phase.

Figure (12) shows the full description of the numerical model, including the: computational domain, boundary conditions, fluid properties and structural parameters of the solid-spring system. The free stream velocity is U = 0.0264 m/s and the Reynolds number based on the cylinder diameter is 2000. Initially, we performed some numerical simulations considering a fixed cylinder and the vortex shedding frequency obtained at Re = 2000 was f_S = 0.69 Hz. This value is higher than that reported by Hiejima et al. (1997), equals to 0.55Hz, but still in good agreement with the experimental curve presented by Blevins, 1986. The cylinder-spring parameters adopted here are that adopted by Hiejima et al. (1997), except the mass which is determined so that the natural frequency of the system is equal to the vortex shedding frequency of 0.69Hz. Therefore, for the values we obtained for the vortex shedding frequency (f_S = 0.69 Hz), this corresponds to the resonance frequency. These differences from the Hiejima et al. (1997) values have to be taken into account when comparing our results with theirs.

The domain was subdivided in subdomains where the different descriptions are utilized in order to facilitate the treatment of problems involving multiple physics as shown in Fig. 2. The finite element mesh adopted consists of a triangularization with 12076 elements and 6122 nodes. The circular cylinder is centered on the origin of the coordinates axes, and the portion of the ALE mesh is restricted to the circular region around the cylinder with diameter 1 m. The domain consists in a rectangle whose geometric control nodes coordinates in meters are: (-0.65; -0.65), (1.80; -0.65), (1.80; 0.65) e (-0.65; 0.65).

The effect of the periodic velocity excitation was investigated considering different ratio between the values of the excitation frequency (f_a) and the vortex shedding frequency (f_S), i.e. f_a / f_S = 3.78; 4.45 and 5.12. Accordingly to Hiejima et al. (1997), the value 4.45 is close to the experimental value near the transition wave frequency, which is an effective value of frequency for an acoustic excitation to change the flow around a stationary circular cylinder. With such value of excitation they were able to get a considerable increase on the vortex shedding frequency that was quite effective in reducing the vortex induced vibration amplitude, as the experimental results suggests. We picked up two other values around 4.45 in order to study the influence of the excitation frequency on our results. Considering a fixed cylinder and the different ratio (f_a/f_s) mentioned previously, the frequency of the velocity transversal to the flow in a point located inside the vortex shedding region behind the cylinder was studied. The frequency of the transversal velocity with and without the periodic excitation are plotted in Fig. 13. It can be observed that the periodic excitation modify this frequency and consequently the vortex shedding frequency.

The same analyses were performed considering the cylinder free to vibrate in the direction transversal to the flow. The numerical simulation set up consists of starting with a fixed cylinder and after the vortex shedding becomes periodic we allow the transverse movement, and after the vibration amplitude stabilizes on a constant value we start applying the periodic excitation. In all three cases the cylinder is set free when the time is around 80 seconds and the excitation starts when the time is around 180 seconds. In Fig. (14), ⁽¹⁵⁾ and ⁽¹⁶⁾ the displacement histories are plotted for f_a / f_S = 3.78; 4.45 and 5.12, respectively. In Fig. (14), f_a / f_S = 3.78, the transversal displacement amplitude reduces and we notice the presence of the pulse phenomena. In ^{Fig. (15)}, with f_a / f_S = 4.45, we have the most effective (for the three values analysed) reduction on the oscillation amplitude. The modification on amplitudes is according to the Bloor's curve, described in Hiejima et al. 1997. This results can be observed in Tab. 2.

For f_a / f_S = 5.12, ^{Fig. (16)}, there is no reduction but a small increase on the oscillatory amplitude, and the adopted excitation frequency has an adverse effect. For f_a / f_S = 4.45 the vortex shedding frequency differs more from the cylinder natural frequency (or resonance frequency) then the other ratio (see Tab. 2) and the results suggest that an even bigger variation on the vortex shedding frequency might reduce more or even suppress the vibration on the cylinder.

The table below shows the results for the different values of the excitation frequency.

It should be observed that the amplitude of the oscillations were small with the cylinder vibrating under the influence of the vortex formation and shedding behind the cylinder, and that the characteristic of the vortex induced vibrations were directly affected by the change on the frequency of such vortex formation and shedding. Qualitatively we were able to reproduce the expected results, but further investigation on the disagreement with Hiejima et al. 1997 on the vortex shedding frequency for the fixed cylinder must be pursued. Also, further investigation considering different ration between the values of the excitation frequency (f_a) and the vortex shedding value (f_S), and also considering different application points and amplitude of excitation must be pursued in order to gain a better insight on the behaviour of this application.

Conclusions

A general description of a computational system for fluid-structure interaction analysis was presented. The system incorporates several important tools (different mesh movement algorithms; possibility of decomposing the domain into several subdomains with different reference frame descriptions and of using global or local adaptive remeshing strategies) which renders it very flexible and capable of dealing with a large class of two dimensional applications. Some of these tools and strategies were not fully exploited in the present applications and must also be tested with more complex problems to stress their real importance. It was demonstrated through the numerical results and model problems presented that the ALE finite element formulation developed can be used to simulate correctly fluid-structure problems involving cylinder in cross flow, where the results for the important phenomenon of 'lock-in" are consistent with the literature. The results obtained through the study of using periodic acoustic excitation to suppress the vortex induced vibration on a circular cylinder are qualitatively consistent with the literature, however they are just preliminary and requires further investigation. The fluid dynamics solver using the adopted formulation has severe limitation in terms of the time step size leading to a large number of iterations and improvements on this feature or an alternative formulation without such a strong constraint must be investigated. Finally, other improvements are required in terms of efficiency to allow the analysis of large-scale problems within an acceptable time. This improvements would involve, for instance, the incorporation of an error estimator to control the adaptive procedure, a parallel or parallel/vector implementation and others elements of high performance computation. Some of this aspects are already under investigation.

Acknowledgement

The authors would like to thank the Brazilian Government, through the Agencies CAPES, CNPq and MCT-ANP for the financial support provided. The authors also wish to thanks Dr. Paulo A. B. de Sampaio for some useful discussions.

Paper accepted June, 2005.

Technical Editor: Aristeu da Silveira Neto.

• Antunes, A. R. E. and Lyra, P. R. M., 2002, "A Finite Element Petrov-Galerkin Formulation for the Solution of Incompressible Viscous Flow Involving Moving Boundaries", Proceedings of the II National Congress of Mechanical Engineering, João Pessoa, Brasil, 10p (in Portuguese and in CD Rom).
• Antunes, A. R. E., 2002, "A Flexible Computational System Using an "ALE" Finite Element Formulation for the Analysis of Fluid-Structure Interaction Problems". Master dissertation, Mechanical Engineering Department, Federal University of Pernambuco (UFPE), Recife/PE, Brasil, 83 p. (in Portuguese).
• Antunes, A. R. E., Carvalho, D. K. E., and Lyra, P. R. M., 2002, "Different Strategies for Mesh Movement and Remeshing of Unstructured Meshes Used When Dealing with Fluid-Structure Interaction Problem", Proceedings of the 9^th Brazilian Congress of Thermal Engineering and Science, Caxambu, Brasil, 11p (in Portuguese and in CD-ROM).
• Azevedo, R. L., 1999, "Análise de Problemas de Interação Fluido- Estrutura Usando o Método dos Elementos Finitos com um Acoplamento Monolítico", (DSc Thesis, Civil Engineering Departament/Federal University of Rio Grande do Sul), Porto Alegre/RS, Brasil, 116p. (in Portuguese).
• Blevins, R. D., 1986, "Flow-Induced Vibration", R. E. Krieger Publishing, Inc, Malabar/Florida, USA, 363 p.
• Blom, J. F. and Leyland, P., 1998, "Consistency Analysis of Fluid-Structure Interaction Algorithms", Proceedings of the European Cong. on Comp. Methods in App. Science and Eng. (ECCOMAS), Barcelona, Spain, 15 p. (in CD-ROM).
• Bonet, J. and Peraire, J., 1990, "An Alternating Digital Tree (ADT) Algorithm for 3-d Geometric Searching and Intersection Problems", Int. J. Num. Meth. Engng., Vol. 31, pp. 1-17.
• Brooks, A. N. and Hughes, T. J. R., 1982, "Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations", Computer Methods in Applied Mechanics and Engineering, Vol. 32, pp. 199-259.
• Childs, S. J., 1999, "The Energetic Implications of Using Deforming Reference Descriptions to Simulate the Motion of Incompressible Newtonian Fluids", Computer Methods in Applied Mechanics and Eng., Vol. 180, pp 219-238.
• Correia, A. C. A., "Computational Simulation of Vortex Formation at Cylinder, Focusing to the Study of Flow Induced Vibration", Mechanical Engng. Department, UFPE, MSc. Dissertation, 2001, 63p (in Portuguese).
• Farhat, C., Degand, C., Koobus and B., Lesoinne, M., 1997, "Torsional Springs for Two-Dimensional Dynamic Unstructured Fluid Meshes", Computer Methods in Applied Mechanics and Engineering, Vol. 163, pp. 231-245.
• Hiejima, S., Nomura, T., Kimura, K. & Fujino, Y., 1997, "Numerical Study on the Suppression of the Vortex-Induced Vibration of a Circular Cylinder by Acoustic Excitation", Journal of Wind Engineering and Industrial Aerodynamics, Vol. 67 & 68, pp 325-335.
• Hughes, T.J.R., Liu, W.K. and Zimmermann, T.K., 1981, "Lagrangian-Eulerian Formulation for Incompressible Viscous Flows", Computer Methods in Applied Mechanics and Engineering, Vol. 29, pp. 329-349.
• Hughes, T. J. R., 1987, "The Finite Element Method", 1Ş Edition, United States of America, Prentice-Hall, Inc., 803p.
• Löhner, R., 2001, "Applied CFD Techniques", 1Ş Edition, Great Britain, John Wiley & Sons, LTD, 366p.
• Lyra, P. R. M. and Carvalho, D. K. E., 2000, "A Flexible Unstructured Mesh Generator for Transient Anisotropic Remeshing", Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering (ECOMAS), Barcelona, Spain (in CD-ROM).
• Lyra, P. R. M., Antunes, A. R. E. and Willmersdorf, R. B., 2003, "A Flexible ALE Finite Element Procedure for the Analysis of Fluid-Structure Applications", Proceedings of the 17^th International Congress of Mechanical Ingeneering, São Paulo, Brasil, 10p (in CD-ROM).
• Nomura, T. and Hughes, T.J.R., 1992, "An Arbitrary Lagrangian-Eulerian Finite Element Method for Interaction of Fluid and a Rigid Body", Computer Methods in Applied Mechanics and Engineering, Vol. 95, pp. 115-138.
• Nomura, T., 1994, "ALE Finite Element Computations of Fluid-Structure Interaction Problems", Computer Methods in Applied Mechanics and Engineering, Vol. 112, pp. 291-308.
• Peraire, J.; Peiró, J. and Morgan, K., 1999, "Advancing Front Grid Generation", In: Handbook of Grid Generation. Florida, CRC Press, pp. 17.1-17.22.
• Sampaio, P.A.B., 1991, "A Petrov-Galerkin Formulation for the Incompressible Navier-Stokes Equations using Equal Order for Velocity and Pressure", Int. J. Numer. Meth. Eng., Vol. 31, pp. 1135-1149.
• Sampaio, P. A. B., Lyra, P. R. M., Morgan and K., Weatherill, N., 1993, "Petrov-Galerkin Solutions of Incompressible Navier-Stokes Equations in Primitive Variables with Adaptative Remeshing", Computer Methods and Applied Mechanics and Engineering, Vol. 106, pp. 143-178.

Publication Dates

• Publication in this collection
  05 Sept 2005
• Date of issue
  Sept 2005