127 Chapter V Fluid-Structure Coupling Methodology The two separate field solvers (fluid and structure) have been presented in the previous two chapters. This chapter now describes how the two solvers are coupled in order to perform fluid-structure interaction simulations. The first part presents a brief review of some numerical algorithms that have been applied for FSI simulations of thin flexible structures. The second part of this chapter details the approach to couple the moving/deforming boundary of the structure to the underlying fluid mesh. The third part briefly describes the transfer of information between the structural nodes and the underlying fluid mesh. The next part explains the FSI coupling scheme used in the present work. The subsequent part presents validation of the proposed FSI algorithm. Finally, the last part of this chapter summarizes this chapter and highlights some contributions of the present work. V.1 Introduction Most works concerning the FSI of thin flexible structures are around the flapping flag problem. This problem has gained a lot of interest because not only it can exhibit a rich set of nonlinear physics, but it also can be applied in areas such as energy harvesting and flow control. A flexible flag immersed in a fluid flow can exhibit a self-induced periodical flapping motion at certain high flow speed. Most studies of the FSI of thin flexible structures are around the conventional flag and the inverted flag configurations. A conventional flag has configuration where the leading edge is clamped and the trailing edge is free to oscillate, whereas in an inverted flag conguration the leading edge is free to move and the trailing edge is clamped. The conventional flag is stabilized by the structural stiffness and viscous drag, while the pressure difference tends to produce the opposite effect (Gurugubelli and Jaiman, 2015). This phenomenon has been implemented in some new surgical methods (Huang, 1995), nuclear plate assemblies (Guo and Paidoussis, 1999), flow control devices (Jaiman et al., 2004; Lucey, 1998) and small-scale energy harvesting systems (Allen and Smits, 2001; Tang et al., 2009). The flapping 128 dynamics of a conventional flexible flag has been extensively studied theoretically, experimentally, and numerically. an extensive literature review on this topic can be found in (Paidoussis, 2004; Shelley and Zhang, 2011). Zhang et al. (2000) is considered the first to initiate the recent rising of interest in this problem. By using a silk thread, Zhang et al. (2000) reported three distinct behavioural modes: a deflected mode, a regular flapping mode, and an irregular flapping mode. Following the work of Zhang et al. (2000), a large number of numerical studies have been conducted on a conventional flag configuration using various methods. For instance, Zhu and Peskin (2002 and 2003) employed the immersed boundary method to simulate the soap-film experiment of Zhang et al. (2000). This experiment was also simulated by Farnell et al. (2004) using a loosely coupled FSI algorithm. Connell and Yue (2007) employed a strongly coupled finite difference formulation and Liu et al. (2014) developed an algorithm based on the high-order finite-element method using arbitrary Eulerian–Lagrangian (ALE) coordinates to study the effects of the mass ratio and Reynolds number on the flapping dynamics. In contrast to the conventional flag configuration, the inverted flag system is more prone to the instabilities (lower critical velocities for flapping mode) and thus has significantly larger strain levels (Michelin and Doaré, 2013; Tang et al., 2009). Since this inverted flag configuration can harvest more electrical energy than its conventional counterpart with LE clamped, the flapping dynamics of an inverted flag in a uniform flow has become a subject of active research experimentally (Kim et al., 2013), analytically (Sader et al., 2016a, 2016b), and numerically (Gilmanov et al., 2015; Goza et al., 2018; Gurugubelli and Jaiman, 2015 and 2019; Nitti et al., 2020; Ryu et al., 2015; Shoele and Mittal, 2016; Tang et al., 2015). Since the experimental work of Kim et al. (2013), a large number of numerical studies have been performed on an inverted flag configuration using various FSI methods. Ryu et al. (2015) and Shoele and Mittal (2016) employed a two- dimensional (2D) immersed boundary method at Re = 250 and Re = 200, 129 respectively. Tang et al. (2015) used a three-dimensional (3D) immersed boundary lattice-Boltzmann method focusing also on low Reynolds number flows. Gilmanov et al. (2015) performed 3D simulation at high Reynolds number Re = 99,505 using curvilinear immersed boundary method (CURVIB) coupled with a rotation-free finite element (FE) model. Gurugubelli and Jaiman (2015) performed a two- dimensional (2D) numerical simulation at Re = 1000 using a second-order scheme combining the arbitrary Lagrangian–Eulerian (ALE) approach and the finite element procedure. They extended this 2D approach to 3D FSI formulation and performed simulation at Re = 30,000 (Gurugubelli and Jaiman, 2019). They used periodic boundary conditions in the spanwise direction to ignore the side vortices for simplicity and for computational efficiency. Goza et al. (2018) employed a 2D immersed boundary method to study the nonlinear stability of inverted flags. Chawdhury and Morgenthal (2019) developed a partitioned FSI solver using a 2D adaptive multiresolution VPM in combination with a 2D corotational beam method to simulate flapping dynamics of the inverted flag using varous flow speeds. They further developed a 3D partitioned FSI algorithm by combining a pseudo 3D VPM with linear shell model (Chawdhury and Morgenthal, 2021). Recently, Nitti et al. (2020) performed 3D simulations at low Reynolds number Re = 200 using immersed boundary method coupled with a finite element formulation in the Isogeometric Analysis (IGA). All of the aforementioned numerical studies on conventional and inverted flag configurations mostly employed a partitioned FSI approach. In addition, most of them are based on traditional velocity-pressure formulation in their fluid descriptions. The only vorticity-based FSI simulation was from Chawdhury and Morgenthal (2019), where a 2D adaptive multiresolution VPM was coupled with a 2D corotational beam method, very similar to the current work. However, there are two noticeable differences between the current 2D FSI solver and their algorithm. Firstly, they employed a non-consistent corotational beam formulation with the implicit Newmark time integration and artificial damping, whereas the proposed 2D structural solver used a consistent and energy-conserving corotational beam formulation. Secondly, they adopted a wavelets-based multiresolution method 130 called MRAG-I2D (Rossinelli et al., 2015), whereas the proposed fluid solver employed a simple multiresolution approach based on patches with varying resolution. To the best of author’s knowledge, the MRAG-I2D has not been extended to 3D applications. This may be the reason why they used a pseudo 3D VPM for their 3D FSI algorithm (Chawdhury and Morgenthal, 2021). This pseudo VPM employed multi-slice approximation in the spanwise direction, implying that their 3D VPM is actually constructed by a series of 2D VPM. Hence, neither the stretching term nor the cross-coupling of flow in between the slices was handled by their flow model. In light of the above discussion, it can be concluded that the current work may constribute to two aspects of the FSI simulations of thin flexible structures. Firstly, even though it is not the first 2D remeshed-VPM ever applied to this type of FSI problems, the proposed 2D coupled solver might become the first application of a 2D multiresolution remeshed-VPM combined with a consistent and energy- conserving 2D corotational beam formulation to FSI simulations of thin flexible structures. Secondly, the recently developed 3D FSI solver is definitely the first application of 3D multiresolution remeshed-VPM to the FSI simulations of this type of structures. V.2 FSI Coupling of an Immersed Body Due to the nature of partitioned solution scheme, there are two independent domains solved by two independent solvers. To couple the two solvers requires the information transfer at the shared interface. The information transferred include the structural displacements transferred from the solid domain to the fluid domain, and fluid forces transferred from the fluid domain to the solid domain. In the general FSI schemes, there are a number of challenges regarding how to transfer the information between the two independent domains. The fluid and solid domains are usually solved by different frameworks. The fluid domain is commonly solved using CFD solver based on finite volume method, which adopts Eulerian frameworks. The solid domain, however, is usually solved by Finite Element Method (FEM), which uses Lagrangian framework. Hence, the two domains are 131 descritized differently leading to information transfer between two non-matching meshes. This problem becomes one of the main difficulties in CFD, especially when dealing with complex geometries or moving boundaries. The immersed method (IM), the method that is adopted in Brinkman Penalization Method (BPM), solves this problem by obtaining the solution of the fluid equations separately from the boundary representation. The IM represents the effect of the boundary by adding a force in the governing fluid equations. In Brinkman penalization method, this force is called as a penalization force. The no-slip boundary condition on the surface of a body is enforced by introducing a source term (force) around the surface of the body, see (Angot et al., 1999). This penalization force can be considered as the imposed kinematic constraint such as no-slip, slip, stationary, moving, or flexible boundary conditions. In the present work, the approach is used to impose the no-slip boundary condition on both rigid and deforming boundaries. In the vortex method, based on the ideas of Lighthill (1963), the concept of vorticity generation at the solid surface is the key for satisfaction of boundary conditions in the vorticity equation. Satisfying no-slip boundary condition directly satisfies the no-through boundary conditions, see Koumoutsakos et al. (1994), since these boundary conditions are linked (linked boundary condition). This has been further employed by Kamemoto (2004), and Cooper and Barba (2009). The BPM employs two separate grids and different reference frames to represent the fluid and the boundary (see Figure V.1). Eulerian grid, which is Cartesian and fixed in space, is used to define the fluid. The boundary, on the other hand, follows the structural framework which is represented by a Lagrangian grid, hence the boundary is free to move over the top of the underlying fluid mesh. Since the grids are independent, the major challenge in the BP method is to ensure accurate and efficient data transfer from one to the other. Furthermore, since the fluid adopts the non-conforming mesh, this data transfer requires boundary/surface tracking and interpolation. The required boundary conditions are implemented into a distribution 132 of penalized velocities inside the body to impose the desired kinematics. The penalization forces are obtained from these penalized velocities. The fluid equations are solved in the whole computational domain, without any special treatment at the boundary interface (other than the imposed penalization field). Even though the structure and the boundary are defined in the Lagrangian frame, there is no requirement to transform the forces obtained from fluid solver (in Eulerian frame) into Lagrangian frame, thanks to the use of Corotational Beam Method (CBM) as the structural solver. Once the fluid forces are transfered to the structural solver, the CBM can directly use the fluid forces without any coordinate transformation as long as the Eulerian frame used in fluid solver and the global frame of reference used in the structural solver have the same orientation. The transformation from global to local lagrangian coordinate system and vice versa are inherently performed in the corotational formulation. Figure V.1 Conforming (left) and non-conforming (right) meshes One of the main benefits of the BPM or IM is that the solution of the fluid equations can be decoupled from the boundary representation. This offers a very simple mesh generation and can handle arbitrarily complex bodies undergoing large deformations as easily as simple rigid ones, without affecting the fluid mesh. This can greatly avoid problems with poor quality cells and costly remeshing procedures, which typically encounter in body-conforming methods. The IM, however, still have certain drawbacks. The interpolation used in velocity/displacement and force transfer smears the boundary over a finite region in space, resulting in diffused interfaces and reduction of accuracy at the boundary. 133 V.3 Interface Information Transfer The solid domain solved by Corotational Beam Method (CBM) is formulated at the mid-surface of the elements. A data mapping is required to transfer the fluid forces generated by penalization technique on the immersed body to the structural nodes. Moreover, to update the no-slip boundary conditions at the interface, the displacement and velocity of structural nodes are mapped through interpolation to the fluid meshes immersed inside the body. As described in Chapter III the flow solver used in the present work is based on mesh-less VPM.