A Quadratic Elasticity Formulation for Dynamic Interacting Structures in Flow

The deformation of slender elastic structures due to the motion of surrounding fluid is a common multiphysics problem encountered in many applications. In this work we detail the development of a numerical model capable of solving such strongly-coupled fluid-structure interaction problems, and analyse the dynamic behaviour of multiple interacting bodies under fluid loading. In most fluid-structure interaction problems the deformation of slender elastic bodies is significant and cannot be described by a purely linear analysis. We present a new formulation to model these larger displacements. By extending the standard modal analysis technique for linear structural analysis, the governing equations and boundary conditions are updated to account for non-linear terms and a new modal formulation with quadratic modes is derived. The quadratic modal approach is tested on standard benchmark problems of increasing complexity and compared with analytical and full non-linear numerical solutions. An analysis of the dynamic interactions between multiple finite plates in various configurations under fluid loading, as well as the effects of the spacing between the structures, is also conducted. Numerical results are compared with theoretical and experimental approaches. The inverted hydrodynamic drafting effect of elastic bodies in an in-line configuration can be confirmed from our numerical simulations.


Introduction
Computational mechanics is a growing discipline which uses computational methods to obtain approximate solutions to problems governed by the principles of mechanics. Fluidstructure interaction (FSI) constitutes a branch of computational mechanics in which there exists an intimate coupling between fluid and structural or solid domains; the behaviour of the system is influenced by the interaction of a moving fluid and a flexible solid structure. Examples of such FSI systems include wing flutter on aircraft, aeroelasticity of turbomachinery blades, flows in elastic pipes and blood vessels, heart valve dynamics, swimming of fish and in the processing of paper. In this work, we make use of a blend of mathematical and computational approaches to study strongly-coupled fluid-structure interaction problems involving long, thin structures.
2 Strongly-coupled fluid-structure interaction approaches Two computational fluid structure interaction (FSI) approaches are implemented in a partitioned manner to model FSI problems. In the first approach a finite volume method is developed and implemented in OpenFOAM by coupling fluid and structural solvers. The finite element method has primarily been used for modelling the mechanics of solids [1]. The finite volume method [2] has traditionally been more dominant in the field of fluid mechanics but has received increased attention for use in solid mechanics over the last two decades. Both schemes can be considered as methods of weighted residuals where they differ in the choice of weighting function [3]. In terms of computational FSI approaches, many recent studies have made use of a single discretisation scheme, either finite element [4][5][6] or finite volume [7][8][9], to solve the entire domain. This simplifies the treatment at the interface of the fluid and solid domains. The fluid is assumed to be viscous, incompressible and isothermal with the governing equations given by the continuity and Navier-Stokes equations in an arbitrary-Lagrangian-Eulerian (ALE) reference frame. The solid is assumed to be a homogeneous isotropic elastic solid undergoing large, non-linear deformation with the motion governed by Cauchy's first equation of motion. In this work, we make use of the finite volume method for discretisation of the entire domain. In the second approach a reducedorder modal approach for solving the structural equations is implemented and coupled with an incompressible finite volume fluid solver in OpenFOAM. Strong coupling of the fluid and solid domains is achieved by means of a fixed-point solver using a dynamic relaxation parameter based on Aitken's method. It is a simple yet highly robust and efficient approach.

Governing equations
The equations governing viscous incompressible isothermal fluid flow in an arbitrary-Lagrangian-Eulerian (ALE) reference frame [10] are given by the continuity and Navier-Stokes equations. To close the governing equations, a constitutive relation for stress is required. We assume a Newtonian fluid to relate the stress to strain.
The partial differential equations that describe a homogeneous isotropic elastic solid undergoing large, non-linear deformation are given by Cauchy's first equation of motion [11]. Assuming an isotropic hyperelastic St. Venant-Kirchhoff material model, the stress-strain relationship is given in terms of the second Piola-Kirchhoff stress. The first Piola-Kirchhoff stress is related to the second Piola-Kirchhoff stress by the deformation gradient, which relates quantities in the undeformed configuration to their counterparts in the deformed configuration. The constitutive stress-strain relationship is given in terms of the Green-Lagrange strain, which is related to the displacement field through the gradients of displacement.
The finite volume method of discretisation is used for the spatial discretisation of the equations above. This is done by placing them into weak form or integrating over an arbitrary control volume and applying the divergence theorem of Gauss.

Modal FSI analysis
Tangential to the full finite volume FSI approach, a reduced-order modal approach to solving the structural equations has been used and coupled with an incompressible finite volume fluid solver in OpenFOAM.
For a modal structural analysis, Lagrange's equations of motion in the general case are used, which provides an expression for the eigenvalues or natural frequencies and eigenvectors or mode shapes.

Fluid-solid coupling
In this work we limit our analyses to incompressible flows and have made use of the incompressible pimpleDyMFoam flow solver. In its unmodified form, the coupling between the fluid and solid is done once at the beginning of every time-step when movement of the mesh occurs. The fluid flow provides a traction onto the structure. This results in a deformation of the structure which in turn affects the fluid flow. Coupling conditions for the traction, displacement and velocity are applied to ensure momentum conservation or force equilibrium at the interface, as well as to enforce the kinematic or geometric continuity and no-slip conditions respectively. However, the above conditions require an iterative procedure such that a strongly-coupled solution is obtained. FSI systems can be solved using a single or monolithic solution method, which is inherently strongly-coupled, or a partitioned solution method, which can be strongly-or weakly-coupled. Monolithic methods ensure stability and convergence of the solution, since all equations are discretised and solved simultaneously. However, this approach may suffer from ill conditioning and convergence is generally slow. Conversely, partitioned approaches allow for the use of two independent solution techniques for solid and fluid equations in isolation. Weakly-coupled partitioned methods may diverge or result in inaccurate solutions when applied to problems where there are strong interactions between solid and fluid domains. The FSI implementation described above is a weakly-coupled one: the modal equations are solved as a boundary condition only once at the beginning of each time-step. If applied to strongly-coupled problems, the solver results in inaccurate solutions or diverges.
A separate coupling algorithm or additional outer iterations between the fluid and solid is therefore required to achieve strong-coupling. The most popular strongly-coupled FSI algorithms either use fixed point iteration or interface Newton-Krylov methods. Both approaches have their own shortcomings: fixed point methods are generally slow to converge as they make use of Gauss-Seidel iterations and methods to speed up convergence are needed, whilst Newton-Raphson methods require the computation of Jacobians that may be difficult to calculate exactly. A strongly-coupled fixed-point solver with dynamic relaxation has been implemented in this work. It is the most basic of the above approaches yet is highly robust and efficient. The fixed point solver uses a relaxation parameter based on Aitken's method. The solver was modified to call the mesh movement algorithm after every iteration, thus allowing a re-calculation of interface forces, displacements and velocities with every iteration hence resulting in a fully-coupled scheme upon convergence.

Two-dimensional flapping beam
Both FSI approaches have been successfully validated and compared on benchmark cases and the results for a common FSI benchmark problem, proposed by Turek and Hron [12], of an elastic beam in the wake of a cylinder undergoing vortex-induced vibration, are included here. The properties of the fluid and solid are as described in [12]. A uniform constant fluid velocity, Reynolds number Re = 100, was applied at the inlet, while pressure was set to zero at the exit. Flow is assumed to be laminar. Snapshots of the beam deflection and velocity contours at various times are shown in Fig. 1, showing large oscillations of the beam in its second mode of vibration as vortices are shed periodically from either side of it. The vertical tip displacements of the beam using the modal FSI solver and the finite volume FSI method are compared with the results of Turek and Hron [12] in Fig. 2. There is a good correlation between these results. The finite volume FSI results, however, required a very fine mesh for the solid domain and consequently a long run-time. Upon further investigation of just a cantilever beam under a tip load it was found that the finite volume method for structural mechanics suffered from a phenomenon of shear locking. When the aspect ratio of the elements (ratio of element width to height) is large, a numerical shear strain component is produced in addition to a bending strain. This absorbs strain energy and results in a decrease in displacement or stiffening, as shown by the results in Fig. 3 (top). In addition, an evaluation of the accuracy of the finite volume scheme shows that the scheme is only first-order accurate, see Fig. 3 (bottom), which also explains the slow convergence in the result.

Stability analysis
The FSI system may become unstable due to the interaction between the fluid-induced pressure and structural rigidity, with flow velocity and structural length being critical parameters that influence the stability of the system. Three distinct theoretical approaches have traditionally been used to study the flapping stability of slender elastic structures immersed in axial flow. These include the traditional hydrodynamic stability theory, which can be viewed as a fluid-centric approach, the structural acoustic or wave-based approach and the aeroelastic theory or structure-centred approach. A number of detailed experimental studies have also been conducted including soap-film, water-tunnel and wind-tunnel experiments. If the flow velocity is below a certain critical value, the structure remains unmoved or in a so-called stretched-straight state. Above the critical flow velocity, the structure undergoes self-sustaining and regular flapping, while at a much greater velocity, irregular and chaotic flapping results. In terms of structural length, it has been shown that below a critical length the structure always remains in the stretched-straight state. Most theoretical studies have been restricted to linear analyses and two dimensions, whilst only a few studies have looked at the effect of multiple slender bodies. The use of computational methods may help to overcome these limitations. The modal FSI approach was validated by application to a two-dimensional, two degreeof-freedom airfoil. Theoretical methods predict a critical velocity of 51.5 m/s. The results from the transient analysis with an inlet velocity of 51.5m/s show constant amplitude flapping (Fig. 4) which corresponds well with theoretical predictions. A rigorous stability analysis of a two-dimensional finite plate under axial flow is conducted. The plate, which is initially at rest, is clamped at the leading edge and free at the trailing edge. A small initial perturbation is applied and the long-term transient behaviour of the plate is analysed. The bending rigidity of the structure stabilises the system, whilst an increase in aerodynamic pressure due to an increase in flow velocity destabilises the system and results in flutter. Following previous work [13][14][15], the critical stability curves can be defined in terms of the non-dimensional flag density and non-dimensional velocity. These critical stability curves are shown in Fig. 5. Also shown on this figure are various points that were studied in this analysis. To aid comparison, the flow is assumed to be inviscid in this case. The critical stability curves shown in Fig. 5 contain branches that correspond to a different mode structure, as also observed experimentally by Eloy et al. [14]. The first branch for M* < 1 corresponds with a flapping of the plate in its second mode of vibration, while the branch for 1 < M* < 8 corresponds with flapping in the third mode and M* > 8 corresponds with the fourth mode. An analysis of the flapping behaviour shows good agreement with the results of Eloy et al. [14] and Michelin et al. [13]. Snapshots of the flapping behaviour of the plate at M* = 0.002 and U* = 9 are shown in Fig. 6. The second mode of vibration dominates this flapping behaviour and can be observed by the single node, or position of zero displacement, along the plate.