1 Chapter I Introduction I.1 Motivation and Background Fluid-structure interaction (FSI) is a phenomenon generally encountered in nature, in which the motion of the structure is induced by the fluid forces, and the flow is also influenced by the structural response. In nature, this phenomenon can be observed in cardiovascular system, swimming modes of fish-like organisms, and flapping of insects. FSI is also widely utilized in a broad range of engineering fields such as aerospace engineering (aeroelasticity, for instance see Kamakoti and Shyy (2004)), biomedical engineering (analysis of cardiovascular system (Peskin, 1972 and 1977; Gilmanov et al., 2015)), and civil engineering (bridge and offshore platform design (Simiu and Scanlan, 1996)). It plays prominent roles to the design processes of bridges, buildings, aircrafts, where the interactions with the flows can exihibit catastrophic failures such as flutter, galloping, Limit-Cycle-Oscilation (LCO), etc. Even though FSI is an important issue for the design safety, this phenomenon can also be used for beneficial purposes, such as in the flutter energy harvesting system, in which an electic current is generated by extracting the energy stored in structure deformations resulted from its interaction with the flow. One particular class of FSI problems that has been increasingly popular over the past decades is FSI involving thin flexible structures thanks to the increasing trend of the use of thin flexible structures, such as for long-span bridges, tall buildings, and light aerial vehicles, etc. These slender structures are susceptible to FSI and their interactions with the flows are of interest since they can help prevent damage of important structures, understand biological locomotion in fluids, and also develop flutter-based power harvesting system. Such problems, however, are yet to be fully understood due to complex behaviours governing the interactions and limitations in measuring and modelling the systems. Several experiments have been conducted to study the interaction of flexible strcutures in fluid flows (Ghisalberti and Nepf, 2006; Py et al., 2005). However, due to the difficulties associated with some experimental methods (Okamoto et al., 2016) and the advances in computer technology, numerical simulation is becoming an attractive approach. Nevertheless, 2 the development of accurate numerical methods to study FSI is still a challenging task. The main difficulties, especially for FSI simulations of thin flexible structures, are: (a) the multidisciplinary nature where the two fields (i.e.fluid and structure) are coupled using a common interface whose position is not known a priori (Gilmanov et al., 2015); (b) the boundary condition enforcement on the thin body interface; (c) the inherently nonlinear nature involving large and complex deformations resulted from the flow interaction with thin and/or highly flexible structures; and (d) the challenges in developing stable and efficient FSI algorithms especially in problems with high fluid to solid mass ratios. There are two general approaches for FSI simulations: the monolithic approach and the partitioned approach. The monolithic approach (Hübner et al., 2004; Michler et al., 2004; Ryzhakov et al., 2010) employs a unified algorithm to solve simultaneously the fluid and structure equations since both systems are treated in the same mathematical framework. This approach is usually very accurate and stable, but requires more efforts to develop the code. On the other hand, the partitioned approach solves the fluid and the structure equations separately with their respective solvers. The FSI solution from this approach is obtained by coupling the flow and structure solvers to exchange the fluid and structure solutions. This approach allows the use of available flow and structural algorithms, thus reducing the code development time. Another well-known classification of the FSI numerical techniques is based upon the conformity of fluid meshes to the boundary interfaces: the conformal methods and non-conformal methods. The first methods treat the interface displacement and velocity as part of the solution. As a result, fluid meshes have to conform to the interface, thus re-meshing is needed in every time step. Arbitrary Lagrangian- Eulerian (ALE) methods (Hirt et al., 1974) are the example of these conformal methods. In contrast, the second methods treat the interface conditions as constraints imposed on the model equations. Hence, the fluid and structure equations can be solved independently with their respective meshes, and mesh- updating is not necessary. Immersed Boundary (IB) methods, which are still of 3 rising interest in computational FSI, are classified as the non-conformal methods. Several types of IB methods have been reviewed by Mittal and Iaccarino (2005) and Sotiropoulos and Yang (2014). Immersed boundary technique was initially developed by Peskin (Peskin, 1972 and 1977) and has since been improved and applied to simulate various types of flows (Choi et al., 2007; Gilmanov and Sotiropoulos, 2005; Kim et al., 2001; Ye et al., 1999), including turbulent flows (Fadlun et al., 2000; Ikeno and Kajishima, 2007; Tseng and Ferziger, 2003; Yang and Balaras, 2006). This method is particularly attractive since they are simple to implement in existing solvers and can handle complex deforming geometries. An example of the immersed methods is the volume penalization method. This method treats the bodies as external forces, and adds penalization forces to the governing equations. The penalization forces are defined such that kinematic boundary conditions on the solid surface are satisfied. This method can simply represent arbitrary solid geometries in fluids on structured mesh. The volume penalization method has been combined with vortex methods (Hejlesen et al., 2015; Kevlahan and Ghidaglia, 2001). Unlike the other numerical methods in Computational Fluid Dynamics (CFD) which use velocity-pressure formulation, vortex method solves the incompressible Navier-Stokes equations in velocity- vorticity formulation. This clearly simplifies the implementation since the pressure field can be independently computed from the solution procedure. This method has been used as a viable alternative to grid-based methods for its ability to preserve rotational flow features and to precisely and naturally simulate flows dominated by the advection phenomenon. Motivated by the fact that standard vortex methods encountered difficulties especially in handling boundary conditions for complex and/or moving geometries, the vortex method has been coupled with the volume penalization method. 4 In the last decade, the resulting method, or known as the Brinkman Vortex penalization method, has been applied to solve various types of flow simulations involving boundaries. Gallizio (2009) and Morency et al. (2012) used this method to simulate the flow over moving bodies. Their works were extended for the application of fluid-structure interactions by Coquerelle and Cottet (2008) to simulate the interaction of an incompressible flow with rigid bodies. Gazzola et al. further developed the penalization method combined to a projection approach and applied it to study obstacles’ motion induced by vortical flow-structures (Gazzola et al., 2012) and the flow dynamics of biologic swimmer (Gazzola et al., 2011, 2012, and 2014; Rees et al., 2013). The Brinkman Penalization Method (BPM) is relatively easy to code and to parallelize (Mimeau and Mortazavi, 2021). However, the major limitation of this method lies in its first-order accuracy. To overcome this weakness, Hejlesen et al. (2015) introduced a two-dimensional (2D) iterative Brinkman penalization (IBP) method capable to accurately enforce no-slip boundary conditions via iterations of the penalty sub-step. The overall computational cost of this method is of reasonable size since it allows the use of larger time steps than the ones usually required in a non-iterative Brinkman penalization. The method has been extended for the simulations of three-dimensional (3D) flows by Spietz et al. (2017). The vortex methods have also been developed at Institut Teknologi Bandung. The development was initially performed by Sepnov (2009) for 2D flow applications, following the works of Ploumhans and Winckelmans (2000). He extended the method for 3D unbounded flow simulations following the works of Ploumhans et al. (2002). The boundary condition was enforced using vortex sheet diffusion while the convection and diffusion process ware simulated using Fast Multipole Method (FMM) and Particle Strength Exchange (PSE), respectively. Widodo (2013) applied the existing 3D vortex method for 3D incompressible flows over a stationary and spinning bluff bodies at low Reynolds number (Re < 300). Dung (2015) further developed a 3D Vortex Particle Method (VPM) combined with the Brinkman penalization method to simulate the impulsively started flow past a 5 sphere and the flow dynamics around an anguiliform swimmer. Djutanta (2016) continued the work of Dung (2015) by employing the iterative version of Brinkman penalization method, following the works of Hejlesen et al. (2015). The efficiency of this 3D Vortex Particle Method was later improved by Canh (2019) by implementing linked-list scheme, thus reducing the costs in the velocity computation and diffusion steps. These works have paved the way for the development of multiresolution remeshed-Vortex Particle Method (MR-VPM) for incompressible flows around thin flexible bodies as will be discussed in the present work. Even though vortex methods have successfully been applied for bluff-bodies and biologic swimmers, there is no noticeable contributions of these methods or specifically the iterative Brinkman penalization VPM for FSI analysis of deformable geometry especially in the cases involving thin flexible structures. This partly due to the difficulty of vortex method to recover the forces distribution on the surface of complex deforming geometries, but also due to the requirement to have finite body volume when penalization technique is used. Therefore, the possibility of analysing FSI of thin flexible structures would allow the iterative Brinkman penalization VPM to study a new class of FSI problems such as the flow- induced flapping of a thin structure. FSI involving thin flexible structures also provides another challenge to the development of elastodynamics solver since it often exhibits geometric nonlinearity where the structure can undergo large deformation with finite rotation. The conventional FEM encounters some drawbacks when dealing with these geometrically nonlinear problems. The difficulties mainly lie on the existence of mesh distortion when the structures undergo large deformation (Durand et al., 2019; Luo, 2008; Wang et al., 2013), requiring adaptive remeshing, leading to high computational cost. One approach to analyzing large deformation structures, which is employed in the present study, is the corotational beam formulation. 6 I.2 Objectives and Methodologies In light of the previous discussion, the main objective of the present work is to develop a numerical algorithm for FSI simulations of flexible thin structures undergoing large displacements via a partitioned approach, where the flow is resolved using the iterative Brinkman penalization for Vortex Particle Method (IBP-VPM) and the structure is computed using a FEM-based elastodynamics solver. To support this, three primary research aims are identified as follows: ñ to extend the development of the IBP-VPM so that it can simulate flow around flexible and relatively thin structures, both for 2D and 3D cases; ñ to develop 2D and 3D structural solvers capable in simulating large displacements of highly flexible structures; ñ to couple the fluid and structural solver using both strong and loose coupling approaches to simulate FSI problems of various fluid-to-structure density ratios. To address these aims, the first task is to modify the VPM code developed at ITB, which is considered as a Lagrangian VPM since the computation of velocity field, diffusion and stretching terms, and advection of the particles are all evaluated on particles. The velocity field is computed by solving the Biot-Savart equation using Fast Multipole Method (FMM). The major drawback of the FMM lies on its large pre-factor resulted from the steps of constructing and traversing the tree data structure, thus reducing the efficiency of the FMM and making it less attractive (Hejlesen et al., 2013). As a result, the computational cost for the translation of multipole expansions becomes expensive especially for 3D simulations (Hejlesen et al., 2013). This also occured in our initial FSI simulations of 2D inverted flag using FMM-VPM. Constant tree data structure construction was required after advection step to make sure that all vortex particles including those located at far- most locations were evaluated. As a consequence, despite having utilizing the linked-list scheme, the simulations were still more expensive than those using remeshed-Vortex Particle Method (R-VPM) such as Vortex-In-Cell (VIC) method in which Fast Fourier Transform (FFT) is employed to resolve the particle velocity. 7 In constrast to the Lagrangian VPM, in the remeshed-VPM (R-VPM) the positions and the vorticity of the particles are updated in Lagrangian manner whereas the particle’s velocities, the stretching term, and the diffusion term, in which pure Lagrangian vortex methods experience difficulties, are computed using the Eulerian framework. This approach is more attractive, easier to implement, and can save the computational expense. A study to compare the performances of the Lagrangian and remeshed-VPMs in terms of accuracy and computational efficiency has been conducted using an impulsively started flow past a circular cylinder, as presented in Appendix A. The simulation concluded that the remeshed-VPM is more efficient and accurate than the Lagrangian VPM. In addition, the calculated forces are also less fluctuated than those obtained using the Lagrangian VPM. An accurate and stable output from the flow solver is very crucial for the partitioned FSI algorithm since it will affect the stability of the FSI simulations. These reasons have paved the way to adopt the remeshed-VPM as the fluid solver in the present work. To handle the immersed thin solid domain, the fluid solver is combined with the multiresolution scheme proposed in (Rasmussen et al., 2011). The scheme which was originally developed for 2D flow simulations is then extended for 3D applications. This multiresolution scheme enables simulating flow around thin structures by deviding the computational domain into two patches, where the patch within the body vicinity has a fine mesh size whereas the other patch uses a coarser mesh size. The finer mesh is required in the patch enclosed the solid boundary to make sure that the immersed body has a thickness of at least one or two cells to correctly enforce the boundary conditions via the penalization method. Finer mesh in the penalization region is also required in order to obtain accurate penalization forces, which is crucial for the present FSI simulations. The implementation of the smaller meshes would definitely increase the computational cost. To overcome this, a coarse mesh is applied for the rest of the computational domain. It should be noted that the present work is not intended to resolve the complex flow patterns resulted from the interaction between a flexible structure with an incompressible flow. 8 The subsequent task is to develop an elastodynamic solver that can accurately simulate long-term 2D and 3D dynamic analysis of flexible structures undergoing large displacements. In the FSI of flexible thin structures, large deformation with finite rotations can be experienced by the structures. Here, the Corotational Beam Formulation (CBF) is employed. This framework has been successfully used by many researchers to develop simple and efficient nonlinear beam and shell finite element methods of flexible structures (Alsafadie et al., 2011; Argyris et al., 1979; Battini, 2007 and 2008; Battini and Pacoste, 2002; Chhang et al., 2017a and 2017b; Crisfield et al., 1997; Crisfield and Cole, 1990; Galvanetto and Crisfield, 1996; Le et al., 2011, 2012, 2014a, and 2014b). For 3D dynamic analysis, however, there is no consistent and energy-momentum preserving algorithm applied to this formulation. This algorithm, which is capable to conserve the energy and momentum of the system, is beneficial for the cases of long-term simulations such as in FSI simulations of an inverted flag system. Inspired by the simple mid-point algorithm in (Crisfield et al., 1997) and the advantages of implementing 3D consistent corotational formulation as reported in (Le et al., 2014a), the present work develops a consistent and approximately energy-conserving 3D corotational formulation for nonlinear dynamic analysis of flexible beam structures. It is important to validate both the fluid and structural solvers before the FSI coupling step to evaluate the accuracy and efficiency of each solver.